From cd677685e9a3662ba44b151b3d4e189cafdeaaac Mon Sep 17 00:00:00 2001 From: Tom Byer <149726499+tjb-ltk@users.noreply.github.com> Date: Tue, 23 Sep 2025 10:38:11 -0700 Subject: [PATCH 1/6] change validation for region statistics for by resvol constraint --- .../wells/CompositionalMultiphaseWell.cpp | 39 +++++++++---------- .../wells/CompositionalMultiphaseWell.hpp | 6 +-- .../fluidFlow/wells/SinglePhaseWell.cpp | 34 ++++++++-------- .../fluidFlow/wells/SinglePhaseWell.hpp | 8 ++-- .../fluidFlow/wells/WellSolverBase.cpp | 7 ++-- .../fluidFlow/wells/WellSolverBase.hpp | 3 +- 6 files changed, 50 insertions(+), 47 deletions(-) diff --git a/src/coreComponents/physicsSolvers/fluidFlow/wells/CompositionalMultiphaseWell.cpp b/src/coreComponents/physicsSolvers/fluidFlow/wells/CompositionalMultiphaseWell.cpp index 14537f0d999..483ac3c5d3c 100644 --- a/src/coreComponents/physicsSolvers/fluidFlow/wells/CompositionalMultiphaseWell.cpp +++ b/src/coreComponents/physicsSolvers/fluidFlow/wells/CompositionalMultiphaseWell.cpp @@ -453,17 +453,7 @@ void CompositionalMultiphaseWell::validateWellConstraints( real64 const & time_n GEOS_FMT( "{}: Region {} is not a target of the reservoir solver and cannot be used for referenceReservoirRegion in WellControl {}.", getDataContext(), regionName, wellControls.getName() ) ); - ElementRegionBase const & region = elemManager.getRegion( wellControls.referenceReservoirRegion()); - // Check if regions statistics are being computed - GEOS_ERROR_IF( !region.hasWrapper( CompositionalMultiphaseStatistics::catalogName()), - GEOS_FMT( "{}: No region average quantities computed. WellControl {} referenceReservoirRegion field requires CompositionalMultiphaseStatistics to be configured for region {} ", - getDataContext(), wellControls.getName(), regionName )); - - CompositionalMultiphaseStatistics::RegionStatistics const & stats = region.getReference< CompositionalMultiphaseStatistics::RegionStatistics >( - CompositionalMultiphaseStatistics::regionStatisticsName() ); - wellControls.setRegionAveragePressure( stats.averagePressure ); - wellControls.setRegionAverageTemperature( stats.averageTemperature ); } } string const & fluidName = subRegion.getReference< string >( viewKeyStruct::fluidNamesString()); @@ -682,7 +672,7 @@ void CompositionalMultiphaseWell::updateBHPForConstraint( WellElementSubRegion & } -void CompositionalMultiphaseWell::updateVolRatesForConstraint( WellElementSubRegion & subRegion ) +void CompositionalMultiphaseWell::updateVolRatesForConstraint( ElementRegionManager const & elemManager, WellElementSubRegion & subRegion ) { GEOS_MARK_FUNCTION; @@ -737,8 +727,16 @@ void CompositionalMultiphaseWell::updateVolRatesForConstraint( WellElementSubReg } else { + ElementRegionBase const & region = elemManager.getRegion( wellControls.referenceReservoirRegion()); + CompositionalMultiphaseStatistics::RegionStatistics const & stats = region.getReference< CompositionalMultiphaseStatistics::RegionStatistics >( + CompositionalMultiphaseStatistics::regionStatisticsName() ); + wellControls.setRegionAveragePressure( stats.averagePressure ); + wellControls.setRegionAverageTemperature( stats.averageTemperature ); + GEOS_ERROR_IF( stats.averagePressure <= 0.0, + GEOS_FMT( "{}: No region average quantities computed. WellControl {} referenceReservoirRegion field requires CompositionalMultiphaseStatistics to be configured for region {} ", + getDataContext(), wellControls.getName(), wellControls.referenceReservoirRegion() )); // If flashPressure is not set by region the value is defaulted to -1 and indicates to use top segment conditions - flashPressure = wellControls.getRegionAveragePressure(); + flashPressure = stats.averagePressure; if( flashPressure < 0.0 ) { // region name not set, use segment conditions @@ -1004,13 +1002,14 @@ void CompositionalMultiphaseWell::updateState( DomainPartition & domain ) MeshLevel & mesh, string_array const & regionNames ) { - mesh.getElemManager().forElementSubRegions< WellElementSubRegion >( regionNames, [&]( localIndex const, - WellElementSubRegion & subRegion ) + ElementRegionManager & elemManager = mesh.getElemManager(); + elemManager.forElementSubRegions< WellElementSubRegion >( regionNames, [&]( localIndex const, + WellElementSubRegion & subRegion ) { WellControls & wellControls = getWellControls( subRegion ); if( wellControls.getWellStatus() == WellControls::Status::OPEN ) { - real64 const maxRegionPhaseVolFrac = updateSubRegionState( subRegion ); + real64 const maxRegionPhaseVolFrac = updateSubRegionState( elemManager, subRegion ); maxPhaseVolFrac = LvArray::math::max( maxRegionPhaseVolFrac, maxPhaseVolFrac ); } } ); @@ -1023,14 +1022,14 @@ void CompositionalMultiphaseWell::updateState( DomainPartition & domain ) } -real64 CompositionalMultiphaseWell::updateSubRegionState( WellElementSubRegion & subRegion ) +real64 CompositionalMultiphaseWell::updateSubRegionState( ElementRegionManager const & elemManager, WellElementSubRegion & subRegion ) { // update properties updateGlobalComponentFraction( subRegion ); // update volumetric rates for the well constraints // note: this must be called before updateFluidModel - updateVolRatesForConstraint( subRegion ); + updateVolRatesForConstraint( elemManager, subRegion ); // update densities, phase fractions, phase volume fractions @@ -1148,7 +1147,7 @@ void CompositionalMultiphaseWell::initializeWells( DomainPartition & domain, rea wellElemCompDens ); // 5) Recompute the pressure-dependent properties - updateSubRegionState( subRegion ); + updateSubRegionState( elemManager, subRegion ); // 6) Estimate the well rates // TODO: initialize rates using perforation rates @@ -1954,7 +1953,7 @@ void CompositionalMultiphaseWell::resetStateToBeginningOfStep( DomainPartition & if( wellControls.isWellOpen( ) ) { - updateSubRegionState( subRegion ); + updateSubRegionState( elemManager, subRegion ); } } ); } ); @@ -2148,7 +2147,7 @@ void CompositionalMultiphaseWell::implicitStepSetup( real64 const & time_n, validateWellConstraints( time_n, dt, subRegion, elemManager ); - updateSubRegionState( subRegion ); + updateSubRegionState( elemManager, subRegion ); } } ) ; diff --git a/src/coreComponents/physicsSolvers/fluidFlow/wells/CompositionalMultiphaseWell.hpp b/src/coreComponents/physicsSolvers/fluidFlow/wells/CompositionalMultiphaseWell.hpp index fd874ec5a0e..819beafa0d8 100644 --- a/src/coreComponents/physicsSolvers/fluidFlow/wells/CompositionalMultiphaseWell.hpp +++ b/src/coreComponents/physicsSolvers/fluidFlow/wells/CompositionalMultiphaseWell.hpp @@ -140,10 +140,10 @@ class CompositionalMultiphaseWell : public WellSolverBase /** * @brief Recompute the volumetric rates that are used in the well constraints + * @param elemManager the well region manager containing the well * @param subRegion the well subregion containing all the primary and dependent fields - * @param targetIndex the targetIndex of the subRegion */ - void updateVolRatesForConstraint( WellElementSubRegion & subRegion ); + void updateVolRatesForConstraint( ElementRegionManager const & elemManager, WellElementSubRegion & subRegion ); /** * @brief Recompute the current BHP pressure @@ -185,7 +185,7 @@ class CompositionalMultiphaseWell : public WellSolverBase */ virtual void updateState( DomainPartition & domain ) override; - virtual real64 updateSubRegionState( WellElementSubRegion & subRegion ) override; + virtual real64 updateSubRegionState( ElementRegionManager const & elemManager, WellElementSubRegion & subRegion ) override; virtual string wellElementDofName() const override { return viewKeyStruct::dofFieldString(); } diff --git a/src/coreComponents/physicsSolvers/fluidFlow/wells/SinglePhaseWell.cpp b/src/coreComponents/physicsSolvers/fluidFlow/wells/SinglePhaseWell.cpp index 4d49a5062e4..9db9bb4ce88 100644 --- a/src/coreComponents/physicsSolvers/fluidFlow/wells/SinglePhaseWell.cpp +++ b/src/coreComponents/physicsSolvers/fluidFlow/wells/SinglePhaseWell.cpp @@ -164,16 +164,6 @@ void SinglePhaseWell::validateWellConstraints( real64 const & time_n, GEOS_FMT( "{}: Region {} is not a target of the reservoir solver and cannot be used for referenceReservoirRegion in WellControl {}.", getDataContext(), regionName, wellControls.getName() ) ); - ElementRegionBase const & region = elemManager.getRegion( wellControls.referenceReservoirRegion() ); - - // Check if regions statistics are being computed - GEOS_ERROR_IF( !region.hasWrapper( SinglePhaseStatistics::catalogName()), - GEOS_FMT( "{}: No region average quantities computed. WellControl {} referenceReservoirRegion field requires SinglePhaseStatistics to be configured for region {} ", - getDataContext(), wellControls.getName(), regionName )); - - SinglePhaseStatistics::RegionStatistics const & stats = region.getReference< SinglePhaseStatistics::RegionStatistics >( SinglePhaseStatistics::regionStatisticsName() ); - wellControls.setRegionAveragePressure( stats.averagePressure ); - wellControls.setRegionAverageTemperature( stats.averageTemperature ); } } WellControls::Control currentControl = wellControls.getControl(); @@ -262,7 +252,7 @@ void SinglePhaseWell::updateBHPForConstraint( WellElementSubRegion & subRegion ) } -void SinglePhaseWell::updateVolRateForConstraint( WellElementSubRegion & subRegion ) +void SinglePhaseWell::updateVolRateForConstraint( ElementRegionManager const & elemManager, WellElementSubRegion & subRegion ) { GEOS_MARK_FUNCTION; @@ -303,8 +293,18 @@ void SinglePhaseWell::updateVolRateForConstraint( WellElementSubRegion & subRegi } else { + ElementRegionBase const & region = elemManager.getRegion( wellControls.referenceReservoirRegion() ); + + // Check if regions statistics are being computed + SinglePhaseStatistics::RegionStatistics const & stats = region.getReference< SinglePhaseStatistics::RegionStatistics >( SinglePhaseStatistics::regionStatisticsName() ); + GEOS_ERROR_IF( stats.averagePressure <= 0.0, + GEOS_FMT( "{}: No region average quantities computed. WellControl {} referenceReservoirRegion field requires SinglePhaseStatistics to be configured for region {} ", + getDataContext(), wellControls.getName(), wellControls.referenceReservoirRegion() )); + wellControls.setRegionAveragePressure( stats.averagePressure ); + wellControls.setRegionAverageTemperature( stats.averageTemperature ); + // use region conditions - flashPressure = wellControls.getRegionAveragePressure(); + flashPressure = stats.averagePressure; if( flashPressure < 0.0 ) { // use segment conditions @@ -401,11 +401,11 @@ void SinglePhaseWell::updateFluidModel( WellElementSubRegion & subRegion ) const } ); } -real64 SinglePhaseWell::updateSubRegionState( WellElementSubRegion & subRegion ) +real64 SinglePhaseWell::updateSubRegionState( ElementRegionManager const & elemManager, WellElementSubRegion & subRegion ) { // update volumetric rates for the well constraints // Warning! This must be called before updating the fluid model - updateVolRateForConstraint( subRegion ); + updateVolRateForConstraint( elemManager, subRegion ); // update density in the well elements updateFluidModel( subRegion ); @@ -491,7 +491,7 @@ void SinglePhaseWell::initializeWells( DomainPartition & domain, real64 const & // 4) Recompute the pressure-dependent properties // Note: I am leaving that here because I would like to use the perforationRates (computed in UpdateState) // to better initialize the rates - updateSubRegionState( subRegion ); + updateSubRegionState( elemManager, subRegion ); string const & fluidName = subRegion.getReference< string >( viewKeyStruct::fluidNamesString() ); SingleFluidBase & fluid = subRegion.getConstitutiveModel< SingleFluidBase >( fluidName ); @@ -1152,7 +1152,7 @@ void SinglePhaseWell::resetStateToBeginningOfStep( DomainPartition & domain ) subRegion.getField< well::connectionRate_n >(); connRate.setValues< parallelDevicePolicy<> >( connRate_n ); - updateSubRegionState( subRegion ); + updateSubRegionState( elemManager, subRegion ); } ); } ); } @@ -1194,7 +1194,7 @@ void SinglePhaseWell::implicitStepSetup( real64 const & time, validateWellConstraints( time, dt, subRegion, elemManager ); - updateSubRegionState( subRegion ); + updateSubRegionState( elemManager, subRegion ); } ); } ); } diff --git a/src/coreComponents/physicsSolvers/fluidFlow/wells/SinglePhaseWell.hpp b/src/coreComponents/physicsSolvers/fluidFlow/wells/SinglePhaseWell.hpp index 7bfeb996aa4..ac73a4a3c76 100644 --- a/src/coreComponents/physicsSolvers/fluidFlow/wells/SinglePhaseWell.hpp +++ b/src/coreComponents/physicsSolvers/fluidFlow/wells/SinglePhaseWell.hpp @@ -140,9 +140,10 @@ class SinglePhaseWell : public WellSolverBase /** * @brief Recompute the volumetric rate that are used in the well constraints + * @param elemManager the well region manager * @param subRegion the well subregion containing all the primary and dependent fields */ - virtual void updateVolRateForConstraint( WellElementSubRegion & subRegion ); + virtual void updateVolRateForConstraint( ElementRegionManager const & elemManager, WellElementSubRegion & subRegion ); /** * @brief Recompute the BHP pressure that is used in the well constraints @@ -164,10 +165,11 @@ class SinglePhaseWell : public WellSolverBase real64 const & dt, DomainPartition & domain ) override; /** - * @brief Recompute all dependent quantities from primary variables (including constitutive models) on the well + * @brief Recompute all dependent quantities from primary variables (including constitutive models) + * @param elemManager the elemManager containing the well * @param subRegion the well subRegion containing the well elements and their associated fields */ - virtual real64 updateSubRegionState( WellElementSubRegion & subRegion ) override; + virtual real64 updateSubRegionState( ElementRegionManager const & elemManager, WellElementSubRegion & subRegion ) override; /** * @brief function to assemble the linear system matrix and rhs diff --git a/src/coreComponents/physicsSolvers/fluidFlow/wells/WellSolverBase.cpp b/src/coreComponents/physicsSolvers/fluidFlow/wells/WellSolverBase.cpp index e470984d7cd..f5136ffb70d 100644 --- a/src/coreComponents/physicsSolvers/fluidFlow/wells/WellSolverBase.cpp +++ b/src/coreComponents/physicsSolvers/fluidFlow/wells/WellSolverBase.cpp @@ -305,9 +305,10 @@ void WellSolverBase::updateState( DomainPartition & domain ) MeshLevel & mesh, string_array const & regionNames ) { - mesh.getElemManager().forElementSubRegions< WellElementSubRegion >( regionNames, [&]( localIndex const, - WellElementSubRegion & subRegion ) - { updateSubRegionState( subRegion ); } ); + ElementRegionManager & elemManager = mesh.getElemManager(); + elemManager.forElementSubRegions< WellElementSubRegion >( regionNames, [&]( localIndex const, + WellElementSubRegion & subRegion ) + { updateSubRegionState( elemManager, subRegion ); } ); } ); } diff --git a/src/coreComponents/physicsSolvers/fluidFlow/wells/WellSolverBase.hpp b/src/coreComponents/physicsSolvers/fluidFlow/wells/WellSolverBase.hpp index a8159378b50..d7634be8b1c 100644 --- a/src/coreComponents/physicsSolvers/fluidFlow/wells/WellSolverBase.hpp +++ b/src/coreComponents/physicsSolvers/fluidFlow/wells/WellSolverBase.hpp @@ -250,9 +250,10 @@ class WellSolverBase : public PhysicsSolverBase /** * @brief Recompute all dependent quantities from primary variables (including constitutive models) + * @param elemManager the elemManager containing the well * @param subRegion the well subRegion containing the well elements and their associated fields */ - virtual real64 updateSubRegionState( WellElementSubRegion & subRegion ) = 0; + virtual real64 updateSubRegionState( ElementRegionManager const & elemManager, WellElementSubRegion & subRegion ) = 0; /** * @brief Recompute the perforation rates for all the wells From 5ced2b7b9cb84cc2231c5a981740cc954a7ad1f5 Mon Sep 17 00:00:00 2001 From: Tom Byer <149726499+tjb-ltk@users.noreply.github.com> Date: Tue, 23 Sep 2025 10:53:06 -0700 Subject: [PATCH 2/6] rel mode compile fix --- .../fluidFlow/wells/CompositionalMultiphaseWell.cpp | 5 ++--- .../fluidFlow/wells/CompositionalMultiphaseWell.hpp | 3 +-- .../physicsSolvers/fluidFlow/wells/SinglePhaseWell.cpp | 5 ++--- .../physicsSolvers/fluidFlow/wells/SinglePhaseWell.hpp | 4 +--- .../physicsSolvers/fluidFlow/wells/WellSolverBase.cpp | 2 +- .../physicsSolvers/fluidFlow/wells/WellSolverBase.hpp | 3 +-- 6 files changed, 8 insertions(+), 14 deletions(-) diff --git a/src/coreComponents/physicsSolvers/fluidFlow/wells/CompositionalMultiphaseWell.cpp b/src/coreComponents/physicsSolvers/fluidFlow/wells/CompositionalMultiphaseWell.cpp index 483ac3c5d3c..fad2d34b505 100644 --- a/src/coreComponents/physicsSolvers/fluidFlow/wells/CompositionalMultiphaseWell.cpp +++ b/src/coreComponents/physicsSolvers/fluidFlow/wells/CompositionalMultiphaseWell.cpp @@ -427,8 +427,7 @@ void CompositionalMultiphaseWell::validateInjectionStreams( WellElementSubRegion void CompositionalMultiphaseWell::validateWellConstraints( real64 const & time_n, real64 const & GEOS_UNUSED_PARAM( dt ), - WellElementSubRegion const & subRegion, - ElementRegionManager const & elemManager ) + WellElementSubRegion const & subRegion ) { WellControls & wellControls = getWellControls( subRegion ); if( !wellControls.useSurfaceConditions() ) @@ -2145,7 +2144,7 @@ void CompositionalMultiphaseWell::implicitStepSetup( real64 const & time_n, MultiFluidBase const & fluid = getConstitutiveModel< MultiFluidBase >( subRegion, fluidName ); fluid.saveConvergedState(); - validateWellConstraints( time_n, dt, subRegion, elemManager ); + validateWellConstraints( time_n, dt, subRegion ); updateSubRegionState( elemManager, subRegion ); } diff --git a/src/coreComponents/physicsSolvers/fluidFlow/wells/CompositionalMultiphaseWell.hpp b/src/coreComponents/physicsSolvers/fluidFlow/wells/CompositionalMultiphaseWell.hpp index 819beafa0d8..38311177552 100644 --- a/src/coreComponents/physicsSolvers/fluidFlow/wells/CompositionalMultiphaseWell.hpp +++ b/src/coreComponents/physicsSolvers/fluidFlow/wells/CompositionalMultiphaseWell.hpp @@ -345,8 +345,7 @@ class CompositionalMultiphaseWell : public WellSolverBase */ virtual void validateWellConstraints( real64 const & time_n, real64 const & dt, - WellElementSubRegion const & subRegion, - ElementRegionManager const & elemManager ) override; + WellElementSubRegion const & subRegion ) override; /** * @brief Create well separator diff --git a/src/coreComponents/physicsSolvers/fluidFlow/wells/SinglePhaseWell.cpp b/src/coreComponents/physicsSolvers/fluidFlow/wells/SinglePhaseWell.cpp index 9db9bb4ce88..943c81677ce 100644 --- a/src/coreComponents/physicsSolvers/fluidFlow/wells/SinglePhaseWell.cpp +++ b/src/coreComponents/physicsSolvers/fluidFlow/wells/SinglePhaseWell.cpp @@ -140,8 +140,7 @@ string SinglePhaseWell::resElementDofName() const void SinglePhaseWell::validateWellConstraints( real64 const & time_n, real64 const & GEOS_UNUSED_PARAM( dt ), - WellElementSubRegion const & subRegion, - ElementRegionManager const & elemManager ) + WellElementSubRegion const & subRegion ) { WellControls & wellControls = getWellControls( subRegion ); if( !wellControls.useSurfaceConditions() ) @@ -1192,7 +1191,7 @@ void SinglePhaseWell::implicitStepSetup( real64 const & time, getConstitutiveModel< SingleFluidBase >( subRegion, subRegion.getReference< string >( viewKeyStruct::fluidNamesString() ) ); fluid.saveConvergedState(); - validateWellConstraints( time, dt, subRegion, elemManager ); + validateWellConstraints( time, dt, subRegion ); updateSubRegionState( elemManager, subRegion ); } ); diff --git a/src/coreComponents/physicsSolvers/fluidFlow/wells/SinglePhaseWell.hpp b/src/coreComponents/physicsSolvers/fluidFlow/wells/SinglePhaseWell.hpp index ac73a4a3c76..8c07e223fb7 100644 --- a/src/coreComponents/physicsSolvers/fluidFlow/wells/SinglePhaseWell.hpp +++ b/src/coreComponents/physicsSolvers/fluidFlow/wells/SinglePhaseWell.hpp @@ -292,12 +292,10 @@ class SinglePhaseWell : public WellSolverBase * @param time_n the time at the beginning of the time step * @param dt the time step dt * @param subRegion the well subRegion - * @param elemManager the element manager */ virtual void validateWellConstraints( real64 const & time_n, real64 const & dt, - WellElementSubRegion const & subRegion, - ElementRegionManager const & elemManager ) override; + WellElementSubRegion const & subRegion ) override; }; diff --git a/src/coreComponents/physicsSolvers/fluidFlow/wells/WellSolverBase.cpp b/src/coreComponents/physicsSolvers/fluidFlow/wells/WellSolverBase.cpp index f5136ffb70d..bba04073fbd 100644 --- a/src/coreComponents/physicsSolvers/fluidFlow/wells/WellSolverBase.cpp +++ b/src/coreComponents/physicsSolvers/fluidFlow/wells/WellSolverBase.cpp @@ -159,7 +159,7 @@ void WellSolverBase::initializePostSubGroups() [&]( localIndex const, WellElementSubRegion & subRegion ) { - validateWellConstraints( 0, 0, subRegion, elemManager ); + validateWellConstraints( 0, 0, subRegion ); // validate perforation status table PerforationData & perforationData = *subRegion.getPerforationData(); diff --git a/src/coreComponents/physicsSolvers/fluidFlow/wells/WellSolverBase.hpp b/src/coreComponents/physicsSolvers/fluidFlow/wells/WellSolverBase.hpp index d7634be8b1c..04fe58112b4 100644 --- a/src/coreComponents/physicsSolvers/fluidFlow/wells/WellSolverBase.hpp +++ b/src/coreComponents/physicsSolvers/fluidFlow/wells/WellSolverBase.hpp @@ -321,8 +321,7 @@ class WellSolverBase : public PhysicsSolverBase */ virtual void validateWellConstraints( real64 const & time_n, real64 const & dt, - WellElementSubRegion const & subRegion, - ElementRegionManager const & elemManager ) = 0; + WellElementSubRegion const & subRegion ) = 0; virtual void printRates( real64 const & time_n, real64 const & dt, From b2cbc13cec17710960dc323cb06b0880517831aa Mon Sep 17 00:00:00 2001 From: Tom Byer <149726499+tjb-ltk@users.noreply.github.com> Date: Tue, 23 Sep 2025 13:45:49 -0700 Subject: [PATCH 3/6] check if reference region defined --- .../wells/CompositionalMultiphaseWell.cpp | 21 +++++++++------- .../fluidFlow/wells/SinglePhaseWell.cpp | 24 ++++++++++--------- 2 files changed, 25 insertions(+), 20 deletions(-) diff --git a/src/coreComponents/physicsSolvers/fluidFlow/wells/CompositionalMultiphaseWell.cpp b/src/coreComponents/physicsSolvers/fluidFlow/wells/CompositionalMultiphaseWell.cpp index fad2d34b505..1493279b08a 100644 --- a/src/coreComponents/physicsSolvers/fluidFlow/wells/CompositionalMultiphaseWell.cpp +++ b/src/coreComponents/physicsSolvers/fluidFlow/wells/CompositionalMultiphaseWell.cpp @@ -726,16 +726,19 @@ void CompositionalMultiphaseWell::updateVolRatesForConstraint( ElementRegionMana } else { - ElementRegionBase const & region = elemManager.getRegion( wellControls.referenceReservoirRegion()); - CompositionalMultiphaseStatistics::RegionStatistics const & stats = region.getReference< CompositionalMultiphaseStatistics::RegionStatistics >( - CompositionalMultiphaseStatistics::regionStatisticsName() ); - wellControls.setRegionAveragePressure( stats.averagePressure ); - wellControls.setRegionAverageTemperature( stats.averageTemperature ); - GEOS_ERROR_IF( stats.averagePressure <= 0.0, - GEOS_FMT( "{}: No region average quantities computed. WellControl {} referenceReservoirRegion field requires CompositionalMultiphaseStatistics to be configured for region {} ", - getDataContext(), wellControls.getName(), wellControls.referenceReservoirRegion() )); + if( wellControls.referenceReservoirRegion() != "" ) + { + ElementRegionBase const & region = elemManager.getRegion( wellControls.referenceReservoirRegion()); + CompositionalMultiphaseStatistics::RegionStatistics const & stats = region.getReference< CompositionalMultiphaseStatistics::RegionStatistics >( + CompositionalMultiphaseStatistics::regionStatisticsName() ); + wellControls.setRegionAveragePressure( stats.averagePressure ); + wellControls.setRegionAverageTemperature( stats.averageTemperature ); + GEOS_ERROR_IF( stats.averagePressure <= 0.0, + GEOS_FMT( "{}: No region average quantities computed. WellControl {} referenceReservoirRegion field requires CompositionalMultiphaseStatistics to be configured for region {} ", + getDataContext(), wellControls.getName(), wellControls.referenceReservoirRegion() )); + } // If flashPressure is not set by region the value is defaulted to -1 and indicates to use top segment conditions - flashPressure = stats.averagePressure; + flashPressure = wellControls.getRegionAveragePressure(); if( flashPressure < 0.0 ) { // region name not set, use segment conditions diff --git a/src/coreComponents/physicsSolvers/fluidFlow/wells/SinglePhaseWell.cpp b/src/coreComponents/physicsSolvers/fluidFlow/wells/SinglePhaseWell.cpp index 943c81677ce..60f09a7f1a7 100644 --- a/src/coreComponents/physicsSolvers/fluidFlow/wells/SinglePhaseWell.cpp +++ b/src/coreComponents/physicsSolvers/fluidFlow/wells/SinglePhaseWell.cpp @@ -292,18 +292,20 @@ void SinglePhaseWell::updateVolRateForConstraint( ElementRegionManager const & e } else { - ElementRegionBase const & region = elemManager.getRegion( wellControls.referenceReservoirRegion() ); - - // Check if regions statistics are being computed - SinglePhaseStatistics::RegionStatistics const & stats = region.getReference< SinglePhaseStatistics::RegionStatistics >( SinglePhaseStatistics::regionStatisticsName() ); - GEOS_ERROR_IF( stats.averagePressure <= 0.0, - GEOS_FMT( "{}: No region average quantities computed. WellControl {} referenceReservoirRegion field requires SinglePhaseStatistics to be configured for region {} ", - getDataContext(), wellControls.getName(), wellControls.referenceReservoirRegion() )); - wellControls.setRegionAveragePressure( stats.averagePressure ); - wellControls.setRegionAverageTemperature( stats.averageTemperature ); - + if( wellControls.referenceReservoirRegion() != "" ) + { + ElementRegionBase const & region = elemManager.getRegion( wellControls.referenceReservoirRegion() ); + + // Check if regions statistics are being computed + SinglePhaseStatistics::RegionStatistics const & stats = region.getReference< SinglePhaseStatistics::RegionStatistics >( SinglePhaseStatistics::regionStatisticsName() ); + GEOS_ERROR_IF( stats.averagePressure <= 0.0, + GEOS_FMT( "{}: No region average quantities computed. WellControl {} referenceReservoirRegion field requires SinglePhaseStatistics to be configured for region {} ", + getDataContext(), wellControls.getName(), wellControls.referenceReservoirRegion() )); + wellControls.setRegionAveragePressure( stats.averagePressure ); + wellControls.setRegionAverageTemperature( stats.averageTemperature ); + } // use region conditions - flashPressure = stats.averagePressure; + flashPressure = wellControls.getRegionAveragePressure(); if( flashPressure < 0.0 ) { // use segment conditions From 7771fbe69877b0a7e6eb3cd06c9cebc0d02089e0 Mon Sep 17 00:00:00 2001 From: MelReyCG Date: Tue, 2 Dec 2025 10:53:54 +0100 Subject: [PATCH 4/6] =?UTF-8?q?=F0=9F=9A=A7=20wip=201?= MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit --- src/CompositionalMultiphaseStatist.plantuml | 86 +++++ ...ew2CompositionalMultiphaseStatist.plantuml | 83 +++++ ...ew3CompositionalMultiphaseStatist.plantuml | 115 +++++++ ...ew4CompositionalMultiphaseStatist.plantuml | 115 +++++++ ...NewCompositionalMultiphaseStatist.plantuml | 119 +++++++ .../physicsSolvers/fluidFlow/CMakeLists.txt | 6 +- ...itionalMultiphaseStatisticsAggregator.cpp} | 220 +++--------- ...sitionalMultiphaseStatisticsAggregator.hpp | 188 ++++++++++ .../CompositionalMultiphaseStatisticsTask.cpp | 320 ++++++++++++++++++ ...CompositionalMultiphaseStatisticsTask.hpp} | 114 +++---- 10 files changed, 1119 insertions(+), 247 deletions(-) create mode 100644 src/CompositionalMultiphaseStatist.plantuml create mode 100644 src/New2CompositionalMultiphaseStatist.plantuml create mode 100644 src/New3CompositionalMultiphaseStatist.plantuml create mode 100644 src/New4CompositionalMultiphaseStatist.plantuml create mode 100644 src/NewCompositionalMultiphaseStatist.plantuml rename src/coreComponents/physicsSolvers/fluidFlow/{CompositionalMultiphaseStatistics.cpp => CompositionalMultiphaseStatisticsAggregator.cpp} (69%) create mode 100644 src/coreComponents/physicsSolvers/fluidFlow/CompositionalMultiphaseStatisticsAggregator.hpp create mode 100644 src/coreComponents/physicsSolvers/fluidFlow/CompositionalMultiphaseStatisticsTask.cpp rename src/coreComponents/physicsSolvers/fluidFlow/{CompositionalMultiphaseStatistics.hpp => CompositionalMultiphaseStatisticsTask.hpp} (56%) diff --git a/src/CompositionalMultiphaseStatist.plantuml b/src/CompositionalMultiphaseStatist.plantuml new file mode 100644 index 00000000000..4133442d469 --- /dev/null +++ b/src/CompositionalMultiphaseStatist.plantuml @@ -0,0 +1,86 @@ + + +class Group #text:00000060 +class ExecutableGroup #text:00000060 +class TaskBase #text:00000060 +class FieldStatisticsBase #text:00000060 +class ObjectManagerBase #text:00000060 +class ElementSubRegionBase #text:00000060 +class ElementSubRegionBase #text:00000060 +class WellSolverBase #text:00000060 +class PhysicsSolverBase #text:00000060 + +class FieldStatisticsBase< SOLVER_T > #text:00000060 { + string m_solverName (user input : "solverName") + string m_outputDir (user input : "outputDir") + bool m_writeCSV (user input : "writeCSV") +} + +Group <|-- ExecutableGroup +ExecutableGroup <|-- TaskBase +TaskBase <|-- FieldStatisticsBase + +Group <|-- ObjectManagerBase +ObjectManagerBase <|-- ElementSubRegionBase +ElementSubRegionBase <|-- CellElementRegion + +WellSolverBase <|-- CompositionalMultiphaseWell +PhysicsSolverBase <|-- WellSolverBase +ExecutableGroup <|-- PhysicsSolverBase + +class CompositionalMultiphaseStatistics { + m_computeCFLNumbers (user input: "computeCFLNumbers") + m_computeRegionStatistics (user input: "computeRegionStatistics") + m_relpermThreshold (user input: "relpermThreshold") + -- + void postInputInitialization() override + void registerDataOnMesh( Group & meshBodies ) override + void execute( time_n, dt, [...], domain ) + void computeRegionStatistics( time, meshLevel, regionNames ) + void computeCFLNumbers( time, dt, domainPartition ) +} +FieldStatisticsBase <|-- CompositionalMultiphaseStatistics : SOLVER_T = CompositionalMultiphaseBase + + +class RegionStatistics { + real64 averagePressure + real64 minPressure + real64 maxPressure + + real64 minDeltaPressure + real64 maxDeltaPressure + + real64 averageTemperature + real64 minTemperature + real64 maxTemperature + + real64 totalPoreVolume + real64 totalUncompactedPoreVolume + array1d phasePoreVolume + + array1d phaseMass + array1d trappedPhaseMass + array1d immobilePhaseMass + array2d componentMass +} +note left of RegionStatistics + Also exists in single-phase version with: + real64 averagePressure + real64 minPressure + real64 maxPressure + + real64 minDeltaPressure + real64 maxDeltaPressure + + real64 averageTemperature + real64 minTemperature + real64 maxTemperature + + real64 totalPoreVolume + real64 totalUncompactedPoreVolume + + real64 totalMass +end note +CompositionalMultiphaseStatistics ....> RegionStatistics : "instanciate\n&\ncomputes\n&\nreads data from" +CompositionalMultiphaseWell ....> RegionStatistics : "reads data\nfrom" +CellElementRegion *--l--> RegionStatistics : ""regionStatistics"\nWrapper" diff --git a/src/New2CompositionalMultiphaseStatist.plantuml b/src/New2CompositionalMultiphaseStatist.plantuml new file mode 100644 index 00000000000..f4449422be1 --- /dev/null +++ b/src/New2CompositionalMultiphaseStatist.plantuml @@ -0,0 +1,83 @@ + + +class Group #text:00000060 +class ExecutableGroup #text:00000060 +class TaskBase #text:00000060 +class FieldStatisticsBase #text:00000060 +class ObjectManagerBase #text:00000060 +class ElementSubRegionBase #text:00000060 +class ElementSubRegionBase #text:00000060 +class WellSolverBase #text:00000060 +class PhysicsSolverBase #text:00000060 + +class FieldStatisticsBase< SOLVER_T > #text:00000060 { + string m_solverName (user input : "solverName") + string m_outputDir (user input : "outputDir") + bool m_writeCSV (user input : "writeCSV") +} + +Group <|-- ExecutableGroup +ExecutableGroup <|-- TaskBase +TaskBase <|-- FieldStatisticsBase + +Group <|-- ObjectManagerBase +ObjectManagerBase <|-- ElementSubRegionBase +ElementSubRegionBase <|-- CellElementRegion + +WellSolverBase <|-- CompositionalMultiphaseWell +PhysicsSolverBase <|-- WellSolverBase +ExecutableGroup <|-- PhysicsSolverBase + +class CompositionalMultiphaseStatisticsTask #text:Green { + bool m_computeCFLNumbers (wrapper: "computeCFLNumbers") + bool m_computeRegionStatistics (wrapper: "computeRegionStatistics") + float m_relpermThreshold (wrapper: "relpermThreshold") + -- + void postInputInitialization() override + void execute( time_n, dt, [...], domain ) +} +note bottom of CompositionalMultiphaseStatisticsTask + User defined executable task for + statistics computation & output +end note +FieldStatisticsBase <|-- CompositionalMultiphaseStatisticsTask : SOLVER_T = CompositionalMultiphaseBase + +class CompositionalMultiphaseStatisticsAggregator #text:Green { + void registerDataOnMesh( Group & meshBodies ) override + void computeRegionStatistics( time, meshLevel, regionNames ) + void computeCFLNumbers( time, dt, domainPartition, m_relpermThreshold ) + void forRegionStatstics( functor ) +} +note bottom of CompositionalMultiphaseStatisticsAggregator + Generated sub-group for instanciating + RegionStatistics over the regions & + doing statistics computation. +end note +CompositionalMultiphaseStatisticsTask *---> CompositionalMultiphaseStatisticsAggregator : "contains one" +CompositionalMultiphaseWell *---> CompositionalMultiphaseStatisticsAggregator : "contains one" + +class RegionStatistics { + real64 averagePressure + real64 minPressure + real64 maxPressure + + real64 minDeltaPressure + real64 maxDeltaPressure + + real64 averageTemperature + real64 minTemperature + real64 maxTemperature + + real64 totalPoreVolume + real64 totalUncompactedPoreVolume + array1d phasePoreVolume + + array1d phaseMass + array1d trappedPhaseMass + array1d immobilePhaseMass + array2d componentMass +} +CompositionalMultiphaseStatisticsAggregator ...> RegionStatistics : "instanciate\n&\ncomputes" +CompositionalMultiphaseStatisticsTask ..> RegionStatistics : "reads data\nfrom" +CompositionalMultiphaseWell ...> RegionStatistics : "reads data\nfrom" +CellElementRegion *--l--> RegionStatistics : ""regionStatistics"\nWrapper" diff --git a/src/New3CompositionalMultiphaseStatist.plantuml b/src/New3CompositionalMultiphaseStatist.plantuml new file mode 100644 index 00000000000..e41f9e36ca1 --- /dev/null +++ b/src/New3CompositionalMultiphaseStatist.plantuml @@ -0,0 +1,115 @@ + + +class Group #text:00000060 +class ExecutableGroup #text:00000060 +class TaskBase #text:00000060 +class FieldStatisticsBase #text:00000060 +class ObjectManagerBase #text:00000060 +class ElementSubRegionBase #text:00000060 +class ElementSubRegionBase #text:00000060 +class WellSolverBase #text:00000060 +class PhysicsSolverBase #text:00000060 + +class FieldStatisticsBase< SOLVER_T > #text:00000060 { + string m_solverName (user input : "solverName") + string m_outputDir (user input : "outputDir") + bool m_writeCSV (user input : "writeCSV") +} + +Group <|-- ExecutableGroup +ExecutableGroup <|-- TaskBase +TaskBase <|-- FieldStatisticsBase + +Group <|-- ObjectManagerBase +ObjectManagerBase <|-- ElementSubRegionBase +ElementSubRegionBase <|-- CellElementRegion + +WellSolverBase <|-- CompositionalMultiphaseWell +PhysicsSolverBase <|-- WellSolverBase +ExecutableGroup <|-- PhysicsSolverBase + +class CompositionalMultiphaseStatisticsTask #text:Green { + bool m_computeCFLNumbers (wrapper: "computeCFLNumbers") + bool m_computeRegionStatistics (wrapper: "computeRegionStatistics") + float m_relpermThreshold (wrapper: "relpermThreshold") + -- + void postInputInitialization() override + void execute( time_n, dt, [...], domain ) +} +note bottom of CompositionalMultiphaseStatisticsTask + User defined executable task for + statistics computation & output +end note +FieldStatisticsBase <|-- CompositionalMultiphaseStatisticsTask : SOLVER_T = CompositionalMultiphaseBase + +class CompositionalMultiphaseStatisticsAggregator #text:Green { + void registerDataOnMesh( Group & meshBodies ) override + void computeRegionStatistics( time, meshLevel, regionNames ) + void computeCFLNumbers( time, dt, domainPartition, m_relpermThreshold ) + void forRegionStatstics( functor ) +} +note bottom of CompositionalMultiphaseStatisticsAggregator + Generated sub-group for instanciating + RegionStatistics over the regions & + doing statistics computation. +end note +CompositionalMultiphaseStatisticsTask *---> CompositionalMultiphaseStatisticsAggregator : "contains one" +CompositionalMultiphaseWell *---> CompositionalMultiphaseStatisticsAggregator : "contains one" + +class RegionStatistics #text:Green { +} +note bottom of RegionStatistics + May be removed in favor of direct statistics +end note +CompositionalMultiphaseStatisticsAggregator ...> RegionStatistics : "instanciate\n&\ncomputes" +CompositionalMultiphaseStatisticsTask ..> RegionStatistics : "reads data\nfrom" +CompositionalMultiphaseWell ...> RegionStatistics : "reads data\nfrom" +CellElementRegion *-l--> RegionStatistics : "'regionStatistics' wrapper" + +class PressureStatistics #text:Green { + real64 averagePressure + real64 minPressure + real64 maxPressure + real64 minDeltaPressure + real64 maxDeltaPressure +} +RegionStatistics *--> PressureStatistics : "generated wrapper\n"pressure"" +PressureStatistics [PressureStatistics] --> RegionStatisticsType : "implements" + +class TemperatureStatistics #text:Green { + real64 averageTemperature + real64 minTemperature + real64 maxTemperature +} +RegionStatistics *--> TemperatureStatistics : "generated wrapper\n"temperature"" +TemperatureStatistics [TemperatureStatistics] --> RegionStatisticsType : "implements" + +class PoreVolumeStatistics #text:Green { + real64 totalPoreVolume + real64 totalUncompactedPoreVolume + array1d phasePoreVolume +} +RegionStatistics *--> PoreVolumeStatistics : "generated wrapper\n"poreVolume"" +PoreVolumeStatistics [PoreVolumeStatistics] --> RegionStatisticsType : "implements" + +class MassStatistics #text:Green { + array1d phaseMass + array1d trappedPhaseMass + array1d immobilePhaseMass + array2d componentMass +} +RegionStatistics *--> MassStatistics : "generated wrapper\n"mass"" +MassStatistics [MassStatistics] --> RegionStatisticsType : "implements" + +interface RegionStatisticsType < STATS_TYPE > #text:Green { + void init() + void setKernelViews( RegionStatisticsKernelViews< STATS_TYPE > & views ) + HOST_DEVICE void collectElementStats( RegionStatisticsKernelViews< STATS_TYPE > const & views ) + void kernelReduce() + void mpiReduce() +} + +class STATS_TYPE::KernelViews < STATS_TYPE > #text:Green { + ( structure of views necessary for collectElementStats() ) +} +RegionStatisticsType [STATS_TYPE] ..> STATS_TYPE::KernelViews \ No newline at end of file diff --git a/src/New4CompositionalMultiphaseStatist.plantuml b/src/New4CompositionalMultiphaseStatist.plantuml new file mode 100644 index 00000000000..8c64115b7a4 --- /dev/null +++ b/src/New4CompositionalMultiphaseStatist.plantuml @@ -0,0 +1,115 @@ + + +class Group #text:00000060 +class ExecutableGroup #text:00000060 +class TaskBase #text:00000060 +class FieldStatisticsBase #text:00000060 +class ObjectManagerBase #text:00000060 +class ElementSubRegionBase #text:00000060 +class ElementSubRegionBase #text:00000060 +class WellSolverBase #text:00000060 +class PhysicsSolverBase #text:00000060 + +class FieldStatisticsBase< SOLVER_T > #text:00000060 { + string m_solverName (user input : "solverName") + string m_outputDir (user input : "outputDir") + bool m_writeCSV (user input : "writeCSV") +} + +Group <|-- ExecutableGroup +ExecutableGroup <|-- TaskBase +TaskBase <|-- FieldStatisticsBase + +Group <|-- ObjectManagerBase +ObjectManagerBase <|-- ElementSubRegionBase +ElementSubRegionBase <|-- CellElementRegion + +WellSolverBase <|-- CompositionalMultiphaseWell +PhysicsSolverBase <|-- WellSolverBase +ExecutableGroup <|-- PhysicsSolverBase + +class CompositionalMultiphaseStatisticsTask #text:Green { + bool m_computeCFLNumbers (wrapper: "computeCFLNumbers") + bool m_computeRegionStatistics (wrapper: "computeRegionStatistics") + float m_relpermThreshold (wrapper: "relpermThreshold") + -- + void postInputInitialization() override + void execute( time_n, dt, [...], domain ) +} +note bottom of CompositionalMultiphaseStatisticsTask + User defined executable task for + statistics computation & output +end note +FieldStatisticsBase <|-- CompositionalMultiphaseStatisticsTask : SOLVER_T = CompositionalMultiphaseBase + +class CompositionalMultiphaseStatisticsAggregator #text:Green { + void registerDataOnMesh( Group & meshBodies ) override + void computeRegionStatistics( time, meshLevel, regionNames ) + void computeCFLNumbers( time, dt, domainPartition, m_relpermThreshold ) + void forRegionStatstics( functor ) +} +note bottom of CompositionalMultiphaseStatisticsAggregator + Generated sub-group for instanciating + RegionStatistics over the regions & + doing statistics computation. +end note +CompositionalMultiphaseStatisticsTask *---> CompositionalMultiphaseStatisticsAggregator : "contains one" +CompositionalMultiphaseWell *---> CompositionalMultiphaseStatisticsAggregator : "contains one" + +class RegionStatisticsData #text:Green { +} +note left of RegionStatisticsData + May be removed in favor of direct statistics +end note +CompositionalMultiphaseStatisticsAggregator ...> RegionStatisticsData : "instanciate\n&\ncomputes" +CompositionalMultiphaseStatisticsTask ..> RegionStatisticsData : "reads data\nfrom" +CompositionalMultiphaseWell ...> RegionStatisticsData : "reads data\nfrom" +CellElementRegion *-l--> RegionStatisticsData : "'regionStatistics' wrapper" + +class PressureStatistics #text:Green { + real64 averagePressure + real64 minPressure + real64 maxPressure + real64 minDeltaPressure + real64 maxDeltaPressure +} +RegionStatisticsData *--> PressureStatistics : "generated wrapper\n"pressure"" +PressureStatistics [PressureStatistics] --> RegionStatisticsType : "implements" + +class TemperatureStatistics #text:Green { + real64 averageTemperature + real64 minTemperature + real64 maxTemperature +} +RegionStatisticsData *--> TemperatureStatistics : "generated wrapper\n"temperature"" +TemperatureStatistics [TemperatureStatistics] --> RegionStatisticsType : "implements" + +class PoreVolumeStatistics #text:Green { + real64 totalPoreVolume + real64 totalUncompactedPoreVolume + array1d phasePoreVolume +} +RegionStatisticsData *--> PoreVolumeStatistics : "generated wrapper\n"poreVolume"" +PoreVolumeStatistics [PoreVolumeStatistics] --> RegionStatisticsType : "implements" + +class MassStatistics #text:Green { + array1d phaseMass + array1d trappedPhaseMass + array1d immobilePhaseMass + array2d componentMass +} +RegionStatisticsData *--> MassStatistics : "generated wrapper\n"mass"" +MassStatistics [MassStatistics] --> RegionStatisticsType : "implements" + +interface RegionStatisticsType < STATS_TYPE > #text:Green { + void init() + void setKernelViews( RegionStatisticsKernelViews< STATS_TYPE > & views ) + HOST_DEVICE void collectElementStats( RegionStatisticsKernelViews< STATS_TYPE > const & views ) + void kernelReduce() + void mpiReduce() +} + +class STATS_TYPE::KernelViews < STATS_TYPE > #text:Green { + ( structure of views necessary for collectElementStats() ) +} +RegionStatisticsType [STATS_TYPE] ..> STATS_TYPE::KernelViews \ No newline at end of file diff --git a/src/NewCompositionalMultiphaseStatist.plantuml b/src/NewCompositionalMultiphaseStatist.plantuml new file mode 100644 index 00000000000..782ff4b0263 --- /dev/null +++ b/src/NewCompositionalMultiphaseStatist.plantuml @@ -0,0 +1,119 @@ + + +class Group #text:00000060 +class ExecutableGroup #text:00000060 +class TaskBase #text:00000060 +class FieldStatisticsBase #text:00000060 +class ObjectManagerBase #text:00000060 +class ElementSubRegionBase #text:00000060 +class ElementSubRegionBase #text:00000060 +class WellSolverBase #text:00000060 +class PhysicsSolverBase #text:00000060 + +class FieldStatisticsBase< SOLVER_T > #text:00000060 { + string m_solverName (user input : "solverName") + string m_outputDir (user input : "outputDir") + bool m_writeCSV (user input : "writeCSV") +} + +Group <|-- ExecutableGroup +ExecutableGroup <|-- TaskBase +TaskBase <|-- FieldStatisticsBase + +Group <|-- ObjectManagerBase +ObjectManagerBase <|-- ElementSubRegionBase +ElementSubRegionBase <|-- CellElementRegion + +WellSolverBase <|-- CompositionalMultiphaseWell +WellSolverBase <|-- SinglePhaseWell +PhysicsSolverBase <|-- WellSolverBase +ExecutableGroup <|-- PhysicsSolverBase + +class CompositionalMultiphaseStatistics { + m_computeCFLNumbers (user input: "computeCFLNumbers") + m_computeRegionStatistics (user input: "computeRegionStatistics") + m_relpermThreshold (user input: "relpermThreshold") + -- + void postInputInitialization() override + void registerDataOnMesh( Group & meshBodies ) override + + void computeRegionStatistics( time, meshLevel, regionNames ) + void computeCFLNumbers( time, dt, domainPartition ) +} +FieldStatisticsBase <|-- CompositionalMultiphaseStatistics : SOLVER_T = CompositionalMultiphaseBase + + +class RegionStatistics { + real64 dataTime + ---- + void computeRegionStatistics< STATS_TYPES... >() +} +RegionStatistics <..l.. CompositionalMultiphaseStatistics +CellElementRegion *--> RegionStatistics : "Generated Group\nregionStatistics" + +note right of RegionStatistics::dataTime + Added to keep track if stats has + been computed for the current time. + To be aware that stats are computer + after solver computation, "dataTime" + must be //timeStep + dt// (taking + into account the //dt// of cuts). +end note + + +class PressureStatistics #text:Green { + real64 averagePressure + real64 minPressure + real64 maxPressure + real64 minDeltaPressure + real64 maxDeltaPressure +} +RegionStatistics *--> PressureStatistics : "Lazily Generated Wrapper\npressure" +PressureStatistics [PressureStatistics] --> RegionStatisticsType : "implements" + +class TemperatureStatistics #text:Green { + real64 averageTemperature + real64 minTemperature + real64 maxTemperature +} +RegionStatistics *--> TemperatureStatistics : "Lazily Generated Wrapper\ntemperature" +TemperatureStatistics [TemperatureStatistics] --> RegionStatisticsType : "implements" + +class PoreVolumeStatistics #text:Green { + real64 totalPoreVolume + real64 totalUncompactedPoreVolume + array1d phasePoreVolume +} +RegionStatistics *--> PoreVolumeStatistics : "Lazily Generated Wrapper\nporeVolume" +PoreVolumeStatistics [PoreVolumeStatistics] --> RegionStatisticsType : "implements" + +class MassStatistics #text:Green { + array1d phaseMass + array1d trappedPhaseMass + array1d immobilePhaseMass + array2d componentMass +} +RegionStatistics *--> MassStatistics : "Lazily Generated Wrapper\nmass" +MassStatistics [MassStatistics] --> RegionStatisticsType : "implements" + +class CFLNumbers #text:Green { + array1d phaseMass + array1d trappedPhaseMass + array1d immobilePhaseMass + array2d componentMass +} +RegionStatistics *--> CFLNumbers : "Lazily Generated Wrapper\ncflNumbers" +CFLNumbers [CFLNumbers] --> RegionStatisticsType : "implements" + +interface RegionStatisticsType < STATS_TYPE > #text:Green { + void init() + void setKernelViews( RegionStatisticsKernelViews< STATS_TYPE > & views ) + HOST_DEVICE void collectElementStats( RegionStatisticsKernelViews< STATS_TYPE > const & views ) + void kernelReduce() + void mpiReduce() +} + +class STATS_TYPE::KernelViews < STATS_TYPE > #text:Green { + ( structure of views necessary for collectElementStats() ) +} +RegionStatisticsType [STATS_TYPE] ..> STATS_TYPE::KernelViews \ No newline at end of file diff --git a/src/coreComponents/physicsSolvers/fluidFlow/CMakeLists.txt b/src/coreComponents/physicsSolvers/fluidFlow/CMakeLists.txt index 28ffa57e587..ebfeb403a4f 100644 --- a/src/coreComponents/physicsSolvers/fluidFlow/CMakeLists.txt +++ b/src/coreComponents/physicsSolvers/fluidFlow/CMakeLists.txt @@ -22,7 +22,8 @@ set( fluidFlowSolvers_headers FlowSolverBaseFields.hpp CompositionalMultiphaseBase.hpp CompositionalMultiphaseBaseFields.hpp - CompositionalMultiphaseStatistics.hpp + CompositionalMultiphaseStatisticsAggregator.hpp + CompositionalMultiphaseStatisticsTask.hpp CompositionalMultiphaseFVM.hpp CompositionalMultiphaseHybridFVM.hpp CompositionalMultiphaseUtilities.hpp @@ -129,7 +130,8 @@ set( fluidFlowSolvers_headers set( fluidFlowSolvers_sources CompositionalMultiphaseBase.cpp CompositionalMultiphaseFVM.cpp - CompositionalMultiphaseStatistics.cpp + CompositionalMultiphaseStatisticsAggregator.cpp + CompositionalMultiphaseStatisticsTask.cpp CompositionalMultiphaseHybridFVM.cpp ImmiscibleMultiphaseFlow.cpp ReactiveCompositionalMultiphaseOBL.cpp diff --git a/src/coreComponents/physicsSolvers/fluidFlow/CompositionalMultiphaseStatistics.cpp b/src/coreComponents/physicsSolvers/fluidFlow/CompositionalMultiphaseStatisticsAggregator.cpp similarity index 69% rename from src/coreComponents/physicsSolvers/fluidFlow/CompositionalMultiphaseStatistics.cpp rename to src/coreComponents/physicsSolvers/fluidFlow/CompositionalMultiphaseStatisticsAggregator.cpp index ed8d8e5eb80..13d4fc78b9c 100644 --- a/src/coreComponents/physicsSolvers/fluidFlow/CompositionalMultiphaseStatistics.cpp +++ b/src/coreComponents/physicsSolvers/fluidFlow/CompositionalMultiphaseStatisticsAggregator.cpp @@ -17,78 +17,36 @@ * @file CompositionalMultiphaseStatistics.cpp */ -#include "CompositionalMultiphaseStatistics.hpp" +#include "CompositionalMultiphaseStatisticsAggregator.hpp" -#include "mesh/DomainPartition.hpp" -#include "constitutive/fluid/multifluid/MultiFluidBase.hpp" -#include "constitutive/relativePermeability/RelativePermeabilityBase.hpp" -#include "constitutive/solid/CoupledSolidBase.hpp" -#include "physicsSolvers/LogLevelsInfo.hpp" -#include "physicsSolvers/fluidFlow/LogLevelsInfo.hpp" +#include "common/logger/Logger.hpp" #include "physicsSolvers/fluidFlow/CompositionalMultiphaseBase.hpp" -#include "physicsSolvers/fluidFlow/CompositionalMultiphaseBaseFields.hpp" -#include "physicsSolvers/fluidFlow/CompositionalMultiphaseHybridFVM.hpp" -#include "physicsSolvers/fluidFlow/FlowSolverBaseFields.hpp" #include "physicsSolvers/fluidFlow/kernels/compositional/StatisticsKernel.hpp" -#include "common/format/table/TableData.hpp" -#include "common/format/table/TableFormatter.hpp" -#include "common/format/table/TableLayout.hpp" +#include "physicsSolvers/LogLevelsInfo.hpp" +#include "constitutive/relativePermeability/RelativePermeabilityBase.hpp" +#include namespace geos { -using namespace constitutive; -using namespace fields; -using namespace dataRepository; - -CompositionalMultiphaseStatistics::CompositionalMultiphaseStatistics( const string & name, - Group * const parent ): - Base( name, parent ), - m_computeCFLNumbers( 0 ), - m_computeRegionStatistics( 1 ) +namespace compositionalMultiphaseStatistics { - registerWrapper( viewKeyStruct::computeCFLNumbersString(), &m_computeCFLNumbers ). - setApplyDefaultValue( 0 ). - setInputFlag( InputFlags::OPTIONAL ). - setDescription( "Flag to decide whether CFL numbers are computed or not" ); - - registerWrapper( viewKeyStruct::computeRegionStatisticsString(), &m_computeRegionStatistics ). - setApplyDefaultValue( 1 ). - setInputFlag( InputFlags::OPTIONAL ). - setDescription( "Flag to decide whether region statistics are computed or not" ); - - registerWrapper( viewKeyStruct::relpermThresholdString(), &m_relpermThreshold ). - setApplyDefaultValue( 1e-6 ). - setInputFlag( InputFlags::OPTIONAL ). - setDescription( "Flag to decide whether a phase is considered mobile (when the relperm is above the threshold) or immobile (when the relperm is below the threshold) in metric 2" ); - - addLogLevel< logInfo::CFL >(); - addLogLevel< logInfo::Statistics >(); -} -void CompositionalMultiphaseStatistics::postInputInitialization() -{ - Base::postInputInitialization(); +using namespace constitutive; +// using namespace fields; +using namespace dataRepository; - if( dynamicCast< CompositionalMultiphaseHybridFVM * >( m_solver ) && m_computeCFLNumbers != 0 ) - { - GEOS_THROW( GEOS_FMT( "{} {}: the option to compute CFL numbers is incompatible with CompositionalMultiphaseHybridFVM", - catalogName(), getDataContext() ), - InputError ); - } -} +StatsAggregator::StatsAggregator(): + m_params(), + m_cflStats(), + m_isRegionStatsEnabled( false ), + m_isCFLNumberEnabled( false ) +{} -void CompositionalMultiphaseStatistics::registerDataOnMesh( Group & meshBodies ) +void StatsAggregator::enableRegionStatistics( dataRepository::Group & meshBodies ) { - // the fields have to be registered in "registerDataOnMesh" (and not later) - // otherwise they cannot be targeted by TimeHistory - - // for now, this guard is needed to avoid breaking the xml schema generation - if( m_solver == nullptr ) - { - return; - } + GEOS_ASSERT_NE( m_solver, nullptr ); m_solver->forDiscretizationOnMeshTargets( meshBodies, [&] ( string const &, MeshLevel & mesh, @@ -99,126 +57,36 @@ void CompositionalMultiphaseStatistics::registerDataOnMesh( Group & meshBodies ) integer const numPhases = m_solver->numFluidPhases(); integer const numComps = m_solver->numFluidComponents(); - // if we have to report region statistics, we have to register them first here - if( m_computeRegionStatistics ) + for( size_t i = 0; i < regionNames.size(); ++i ) { - for( size_t i = 0; i < regionNames.size(); ++i ) - { - ElementRegionBase & region = elemManager.getRegion( regionNames[i] ); - - region.registerWrapper< RegionStatistics >( viewKeyStruct::regionStatisticsString() ). - setRestartFlags( RestartFlags::NO_WRITE ); - region.excludeWrappersFromPacking( { viewKeyStruct::regionStatisticsString() } ); - RegionStatistics & stats = region.getReference< RegionStatistics >( viewKeyStruct::regionStatisticsString() ); - - stats.phasePoreVolume.resizeDimension< 0 >( numPhases ); - stats.phaseMass.resizeDimension< 0 >( numPhases ); - stats.trappedPhaseMass.resizeDimension< 0 >( numPhases ); - stats.immobilePhaseMass.resizeDimension< 0 >( numPhases ); - stats.componentMass.resizeDimension< 0, 1 >( numPhases, numComps ); - - if( m_writeCSV > 0 && MpiWrapper::commRank() == 0 ) - { - auto addStatsValue = []( std::ostringstream & pstatsLayout, TableLayout & ptableLayout, - string const & description, string_view punit, - integer pnumPhases, integer pnumComps = 0 ) - { - for( int ip = 0; ip < pnumPhases; ++ip ) - { - if( pnumComps == 0 ) - { - pstatsLayout << description << " (phase " << ip << ") [" << punit << "]"; - } - else - { - for( int ic = 0; ic < pnumComps; ++ic ) - { - pstatsLayout << description << " (component " << ic << " / phase " << ip << ") [" << punit << "]"; - if( ic == 0 ) - { - pstatsLayout << ","; - } - } - } - if( ip == 0 ) - { - pstatsLayout << ","; - } - } - - ptableLayout.addColumn( pstatsLayout.str()); - pstatsLayout.str( "" ); - }; - - string_view massUnit = units::getSymbol( m_solver->getMassUnit() ); - - TableLayout tableLayout( { - TableLayout::Column().setName( "Time [s]" ), - TableLayout::Column().setName( "Min pressure [Pa]" ), - TableLayout::Column().setName( "Average pressure [Pa]" ), - TableLayout::Column().setName( "Max pressure [Pa]" ), - TableLayout::Column().setName( "Min delta pressure [Pa]" ), - TableLayout::Column().setName( "Max delta pressure [Pa]" ), - TableLayout::Column().setName( "Min temperature [Pa]" ), - TableLayout::Column().setName( "Average temperature [Pa]" ), - TableLayout::Column().setName( "Max temperature [Pa]" ), - TableLayout::Column().setName( "Total dynamic pore volume [rm^3]" ), - } ); - - std::ostringstream statsLayout; - addStatsValue( statsLayout, tableLayout, "Phase dynamic pore volume", "rm^3", numPhases ); - addStatsValue( statsLayout, tableLayout, "Phase mass", massUnit, numPhases ); - addStatsValue( statsLayout, tableLayout, "Trapped phase mass (metric 1)", massUnit, numPhases ); - addStatsValue( statsLayout, tableLayout, "Non-trapped phase mass (metric 1)", massUnit, numPhases ); - addStatsValue( statsLayout, tableLayout, "Immobile phase mass (metric 2)", massUnit, numPhases ); - addStatsValue( statsLayout, tableLayout, "Mobile phase mass (metric 2)", massUnit, numPhases ); - addStatsValue( statsLayout, tableLayout, "Component mass", massUnit, numPhases, numComps ); - - std::ofstream outputFile( m_outputDir + "/" + regionNames[i] + ".csv" ); - TableCSVFormatter csvFormatter( tableLayout ); - outputFile << csvFormatter.headerToString(); - } - } + ElementRegionBase & region = elemManager.getRegion( regionNames[i] ); + + region.registerWrapper< RegionStatistics >( viewKeyStruct::regionStatisticsString() ). + setRestartFlags( RestartFlags::NO_WRITE ); + region.excludeWrappersFromPacking( { viewKeyStruct::regionStatisticsString() } ); + RegionStatistics & stats = region.getReference< RegionStatistics >( viewKeyStruct::regionStatisticsString() ); + + stats.phasePoreVolume.resizeDimension< 0 >( numPhases ); + stats.phaseMass.resizeDimension< 0 >( numPhases ); + stats.trappedPhaseMass.resizeDimension< 0 >( numPhases ); + stats.immobilePhaseMass.resizeDimension< 0 >( numPhases ); + stats.componentMass.resizeDimension< 0, 1 >( numPhases, numComps ); } } ); - - // if we have to compute CFL numbers later, we need to register additional variables - if( m_computeCFLNumbers ) - { - m_solver->registerDataForCFL( meshBodies ); - } + m_isRegionStatsEnabled = true; } -bool CompositionalMultiphaseStatistics::execute( real64 const time_n, - real64 const dt, - integer const GEOS_UNUSED_PARAM( cycleNumber ), - integer const GEOS_UNUSED_PARAM( eventCounter ), - real64 const GEOS_UNUSED_PARAM( eventProgress ), - DomainPartition & domain ) +void StatsAggregator::enableCFLStatistics( dataRepository::Group & meshBodies ) { - m_solver->forDiscretizationOnMeshTargets( domain.getMeshBodies(), [&] ( string const &, - MeshLevel & mesh, - string_array const & regionNames ) - { - if( m_computeRegionStatistics ) - { - // current time is time_n + dt - computeRegionStatistics( time_n + dt, mesh, regionNames ); - } - } ); + GEOS_ASSERT_NE( m_solver, nullptr ); - if( m_computeCFLNumbers ) - { - // current time is time_n + dt - computeCFLNumbers( time_n + dt, dt, domain ); - } - - return false; + m_solver->registerDataForCFL( meshBodies ); + m_isCFLNumberEnabled = true; } -void CompositionalMultiphaseStatistics::computeRegionStatistics( real64 const time, - MeshLevel & mesh, - string_array const & regionNames ) const +void StatsAggregator::computeDiscretizationStatistics( real64 const time, + MeshLevel & mesh, + string_array const & regionNames ) const { GEOS_MARK_FUNCTION; @@ -260,11 +128,11 @@ void CompositionalMultiphaseStatistics::computeRegionStatistics( real64 const ti arrayView1d< integer const > const elemGhostRank = subRegion.ghostRank(); arrayView1d< real64 const > const volume = subRegion.getElementVolume(); - arrayView1d< real64 const > const pres = subRegion.getField< flow::pressure >(); - arrayView1d< real64 const > const temp = subRegion.getField< flow::temperature >(); + arrayView1d< real64 const > const pres = subRegion.getField< fields::flow::pressure >(); + arrayView1d< real64 const > const temp = subRegion.getField< fields::flow::temperature >(); arrayView2d< real64 const, compflow::USD_PHASE > const phaseVolFrac = - subRegion.getField< flow::phaseVolumeFraction >(); - arrayView1d< real64 const > const deltaPres = subRegion.getField< flow::deltaPressure >(); + subRegion.getField< fields::flow::phaseVolumeFraction >(); + arrayView1d< real64 const > const deltaPres = subRegion.getField< fields::flow::deltaPressure >(); Group const & constitutiveModels = subRegion.getGroup( ElementSubRegionBase::groupKeyStruct::constitutiveModelsString() ); @@ -306,7 +174,7 @@ void CompositionalMultiphaseStatistics::computeRegionStatistics( real64 const ti launch< parallelDevicePolicy<> >( subRegion.size(), numComps, numPhases, - m_relpermThreshold, + m_params.m_relpermThreshold, elemGhostRank, volume, pres, @@ -553,4 +421,6 @@ REGISTER_CATALOG_ENTRY( TaskBase, CompositionalMultiphaseStatistics, string const &, dataRepository::Group * const ) +} /* namespace compositionalMultiphaseStatistics */ + } /* namespace geos */ diff --git a/src/coreComponents/physicsSolvers/fluidFlow/CompositionalMultiphaseStatisticsAggregator.hpp b/src/coreComponents/physicsSolvers/fluidFlow/CompositionalMultiphaseStatisticsAggregator.hpp new file mode 100644 index 00000000000..419a7e6d35f --- /dev/null +++ b/src/coreComponents/physicsSolvers/fluidFlow/CompositionalMultiphaseStatisticsAggregator.hpp @@ -0,0 +1,188 @@ +/* + * ------------------------------------------------------------------------------------------------------------ + * SPDX-License-Identifier: LGPL-2.1-only + * + * Copyright (c) 2016-2024 Lawrence Livermore National Security LLC + * Copyright (c) 2018-2024 TotalEnergies + * Copyright (c) 2018-2024 The Board of Trustees of the Leland Stanford Junior University + * Copyright (c) 2023-2024 Chevron + * Copyright (c) 2019- GEOS/GEOSX Contributors + * All rights reserved + * + * See top level LICENSE, COPYRIGHT, CONTRIBUTORS, NOTICE, and ACKNOWLEDGEMENTS files for details. + * ------------------------------------------------------------------------------------------------------------ + */ + +/** + * @file CompositionalMultiphaseStatistics.hpp + */ + +#ifndef SRC_CORECOMPONENTS_PHYSICSSOLVERS_FLUIDFLOW_COMPOSITIONALMULTIPHASESTATISTICSAGGREGATOR_HPP_ +#define SRC_CORECOMPONENTS_PHYSICSSOLVERS_FLUIDFLOW_COMPOSITIONALMULTIPHASESTATISTICSAGGREGATOR_HPP_ + +#include "dataRepository/Group.hpp" +#include "mesh/DomainPartition.hpp" +#include "mesh/MeshLevel.hpp" + +namespace geos +{ + +class CompositionalMultiphaseBase; + +namespace compositionalMultiphaseStatistics +{ + +struct AggregatorParameters +{ + /// Threshold to decide whether a phase is considered "mobile" or not + real64 m_relpermThreshold; + + // TODO: add other params like views and stuff +}; + +struct RegionStatistics +{ + /// Time of statistics computation + real64 time; + + /// average region pressure + real64 averagePressure; + /// minimum region pressure + real64 minPressure; + /// maximum region pressure + real64 maxPressure; + + /// minimum region delta pressure + real64 minDeltaPressure; + /// maximum region delta pressure + real64 maxDeltaPressure; + + /// average region temperature + real64 averageTemperature; + /// minimum region temperature + real64 minTemperature; + /// maximum region temperature + real64 maxTemperature; + + /// total region pore volume + real64 totalPoreVolume; + /// total region uncompacted pore volume + real64 totalUncompactedPoreVolume; + /// phase region phase pore volume + array1d< real64 > phasePoreVolume; + + /// region phase mass (trapped and non-trapped, immobile and mobile) + array1d< real64 > phaseMass; + /// trapped region phase mass + array1d< real64 > trappedPhaseMass; + /// immobile region phase mass + array1d< real64 > immobilePhaseMass; + /// region component mass + array2d< real64 > componentMass; + + // TODO: -> split to struct PressureStats...MassStats: + // - VKS for struct name ("pressureStats"..."massStats") + // - current RegionStatistics struct bits +}; + +struct CFLStatistics +{ + /// Time of statistics computation + real64 time; + + /// Maximum Courant Friedrichs Lewy number in the grid for each phase + real64 maxPhaseCFL; + + /// Maximum Courant-Friedrichs-Lewy number in the grid for each component + real64 maxCompCFL; +}; + +class StatsAggregator +{ +public: + + struct viewKeyStruct + { + /// String for the region statistics + constexpr static char const * regionStatisticsString() { return "regionStatistics"; } + }; + + using RegionFunctor = std::function< void (string, RegionStatistics const &) >; + + StatsAggregator(); + + /** + * @brief Set the reference flow solver object to retrieve: + - the simulated regions, + - fields for statistics computation. + * @param solver The reference flow solver + */ + void setFlowSolver( CompositionalMultiphaseBase & solver ) + { m_solver = &solver; } + + // void forDiscretizations( DomainPartition const &, + // std::function< void(MeshLevel const &, + // string_array const & regionNames) > functor ) const; + + // void forRegionStatistics( DomainPartition const &, + // MeshLevel const & mesh, + // RegionFunctor functor ) const; + + /** + * @brief Register the results structs & wrappers so they will be targeted by TimeHistory output + * @note Must be called in "registerDataOnMesh" initialization phase + * @param solver The flow solver + * @param meshBodies The Group containing the MeshBody objects + */ + void enableRegionStatistics( dataRepository::Group & meshBodies ); + + /** + * @brief Register the results structs & wrappers so they will be targeted by TimeHistory output + * @note Must be called in "registerDataOnMesh" initialization phase + * @param solver The flow solver + * @param meshBodies The Group containing the MeshBody objects + */ + void enableCFLStatistics( dataRepository::Group & meshBodies ); + + /** + * @brief Compute some statistics on a given mesh discretization (average field pressure, etc) + * @param[in] time current time + * @param[in] mesh the mesh level object + * @param[in] regionNames the array of target region names + */ + void computeDiscretizationStatistics( real64 const time, + MeshLevel & mesh, + string_array const & regionNames ) const; + + /** + * @brief Compute CFL numbers + * @param[in] time current time + * @param[in] dt the time step size + * @param[in] domain the domain partition + */ + void computeCFLNumbers( real64 const time, + real64 const dt, + DomainPartition & domain ) const; + + CFLStatistics const & getCFLStatistics() const + { return m_cflStats; } + +private: + + CompositionalMultiphaseBase * m_solver; + + AggregatorParameters m_params; + + CFLStatistics m_cflStats; + + bool m_isRegionStatsEnabled; + + bool m_isCFLNumberEnabled; + +}; + +} /* namespace compositionalMultiphaseStatistics */ + +} /* namespace geos */ + +#endif /* SRC_CORECOMPONENTS_PHYSICSSOLVERS_FLUIDFLOW_COMPOSITIONALMULTIPHASESTATISTICSAGGREGATOR_HPP_ */ diff --git a/src/coreComponents/physicsSolvers/fluidFlow/CompositionalMultiphaseStatisticsTask.cpp b/src/coreComponents/physicsSolvers/fluidFlow/CompositionalMultiphaseStatisticsTask.cpp new file mode 100644 index 00000000000..b3d6a6b59b8 --- /dev/null +++ b/src/coreComponents/physicsSolvers/fluidFlow/CompositionalMultiphaseStatisticsTask.cpp @@ -0,0 +1,320 @@ +/* + * ------------------------------------------------------------------------------------------------------------ + * SPDX-License-Identifier: LGPL-2.1-only + * + * Copyright (c) 2016-2024 Lawrence Livermore National Security LLC + * Copyright (c) 2018-2024 TotalEnergies + * Copyright (c) 2018-2024 The Board of Trustees of the Leland Stanford Junior University + * Copyright (c) 2023-2024 Chevron + * Copyright (c) 2019- GEOS/GEOSX Contributors + * All rights reserved + * + * See top level LICENSE, COPYRIGHT, CONTRIBUTORS, NOTICE, and ACKNOWLEDGEMENTS files for details. + * ------------------------------------------------------------------------------------------------------------ + */ + +/** + * @file CompositionalMultiphaseStatistics.cpp + */ + +#include "CompositionalMultiphaseStatisticsTask.hpp" + +#include "common/DataTypes.hpp" +#include "common/StdContainerWrappers.hpp" +#include "common/format/Format.hpp" +#include "mesh/DomainPartition.hpp" +#include "constitutive/fluid/multifluid/MultiFluidBase.hpp" +#include "constitutive/relativePermeability/RelativePermeabilityBase.hpp" +#include "constitutive/solid/CoupledSolidBase.hpp" +#include "physicsSolvers/LogLevelsInfo.hpp" +#include "physicsSolvers/fluidFlow/CompositionalMultiphaseBase.hpp" +#include "physicsSolvers/fluidFlow/CompositionalMultiphaseBaseFields.hpp" +#include "physicsSolvers/fluidFlow/CompositionalMultiphaseHybridFVM.hpp" +#include "physicsSolvers/fluidFlow/FlowSolverBaseFields.hpp" +#include "common/format/table/TableData.hpp" +#include "common/format/table/TableFormatter.hpp" +#include "common/format/table/TableLayout.hpp" +#include + + +namespace geos +{ + +using namespace constitutive; +using namespace fields; +using namespace dataRepository; + +namespace compositionalMultiphaseStatistics +{ + +StatsTask::StatsTask( const string & name, Group * const parent ): + Base( name, parent ), + m_aggregator(), + m_computeCFLNumbers( 0 ), + m_computeRegionStatistics( 1 ) +{ + registerWrapper( viewKeyStruct::computeCFLNumbersString(), &m_computeCFLNumbers ). + setApplyDefaultValue( 0 ). + setInputFlag( InputFlags::OPTIONAL ). + setDescription( "Flag to decide whether CFL numbers are computed or not" ); + + registerWrapper( viewKeyStruct::computeRegionStatisticsString(), &m_computeRegionStatistics ). + setApplyDefaultValue( 1 ). + setInputFlag( InputFlags::OPTIONAL ). + setDescription( "Flag to decide whether region statistics are computed or not" ); + + registerWrapper( viewKeyStruct::relpermThresholdString(), &m_relpermThreshold ). + setApplyDefaultValue( 1e-6 ). + setInputFlag( InputFlags::OPTIONAL ). + setDescription( "Flag to decide whether a phase is considered mobile (when the relperm is above the threshold) or immobile (when the relperm is below the threshold) in metric 2" ); + + addLogLevel< logInfo::CFL >(); + addLogLevel< logInfo::Statistics >(); +} + +void StatsTask::postInputInitialization() +{ + Base::postInputInitialization(); + + m_aggregator.setFlowSolver( *m_solver ); + + if( dynamicCast< CompositionalMultiphaseHybridFVM * >( m_solver ) && m_computeCFLNumbers != 0 ) + { + GEOS_THROW( GEOS_FMT( "{} {}: the option to compute CFL numbers is incompatible with CompositionalMultiphaseHybridFVM", + catalogName(), getDataContext() ), + InputError ); + } +} + +void StatsTask::registerDataOnMesh( Group & meshBodies ) +{ + // for now, this guard is needed to avoid breaking the xml schema generation + if( m_solver == nullptr ) + { + return; + } + + prepareFluidMetaData(); + + if( m_computeRegionStatistics ) + m_aggregator.enableRegionStatistics( meshBodies ); + + // if we have to compute CFL numbers later, we need to register additional variables + if( m_computeCFLNumbers ) + m_solver->registerDataForCFL( meshBodies ); + + m_solver->forDiscretizationOnMeshTargets( meshBodies, [&] ( string const &, + MeshLevel & mesh, + string_array const & regionNames ) + { + prepareLogTableLayouts( mesh.getName() ); + prepareCsvTableLayouts( mesh.getName() ); + } ); +} + +void StatsTask::prepareFluidMetaData() +{ + ConstitutiveManager const & constitutiveManager = this->getGroupByPath< ConstitutiveManager >( "/Problem/domain/Constitutive" ); + MultiFluidBase const & fluid = constitutiveManager.getGroup< MultiFluidBase >( m_solver->referenceFluidModelName() ); + + m_fluid.m_numPhases = fluid.numFluidPhases(); + m_fluid.m_numComps = fluid.numFluidComponents(); + + m_fluid.m_phaseNames = fluid.phaseNames(); + m_fluid.m_compNames = fluid.componentNames(); + + m_fluid.m_phaseCompNames.resize( + m_fluid.m_numPhases, + stdVector< string >( m_fluid.m_numComps, {} ) ); + + for( int ip = 0; ip < m_fluid.m_numPhases; ++ip ) + for( int ic = 0; ic < m_fluid.m_numComps; ++ic ) + m_fluid.m_phaseCompNames[ip][ic] = GEOS_FMT( "{}/{}", m_fluid.m_phaseNames[ip], m_fluid.m_compNames[ic] ); +} + +void StatsTask::prepareLogTableLayouts( string_view meshName ) +{ + // only output from rank 0 + if( MpiWrapper::commRank() != 0 ) + return; + + TableLayout const tableLayout = TableLayout() + .setTitle( GEOS_FMT( "{}: mesh {}", getName(), meshName ) ); + + m_csvFormatters.emplace( meshName, std::make_unique< TableTextFormatter >( tableLayout ) ); +} + +void StatsTask::prepareCsvTableLayouts( string_view meshName ) +{ + // only output from rank 0 + if( MpiWrapper::commRank() != 0 || !m_writeCSV ) + return; + + integer const numPhases = m_solver->numFluidPhases(); + integer const numComps = m_solver->numFluidComponents(); + + auto addPhaseColumns = []( TableLayout & ptableLayout, string const & description, string_view punit, + integer pnumPhases ) + { + for( int ip = 0; ip < pnumPhases; ++ip ) + ptableLayout.addColumn( GEOS_FMT( "{} (phase {}) [{}]", + description, ip, punit ) ); + }; + auto addPhaseCompColumns = []( TableLayout & ptableLayout, string const & description, string_view punit, + integer pnumPhases, integer pnumComps ) + { + for( int ip = 0; ip < pnumPhases; ++ip ) + for( int ic = 0; ic < pnumComps; ++ic ) + ptableLayout.addColumn( GEOS_FMT( "{} (component {} / phase {}) [{}]", + description, ic, ip, punit ) ); + }; + + string_view massUnit = units::getSymbol( m_solver->getMassUnit() ); + + TableLayout tableLayout( { + TableLayout::Column().setName( "Time [s]" ), + TableLayout::Column().setName( "Region" ), // TODO : mention this change in PR description + TableLayout::Column().setName( "Min pressure [Pa]" ), + TableLayout::Column().setName( "Average pressure [Pa]" ), + TableLayout::Column().setName( "Max pressure [Pa]" ), + TableLayout::Column().setName( "Min delta pressure [Pa]" ), + TableLayout::Column().setName( "Max delta pressure [Pa]" ), + TableLayout::Column().setName( "Min temperature [Pa]" ), + TableLayout::Column().setName( "Average temperature [Pa]" ), + TableLayout::Column().setName( "Max temperature [Pa]" ), + TableLayout::Column().setName( "Total dynamic pore volume [rm^3]" ), + } ); + addPhaseColumns( tableLayout, "Phase dynamic pore volume", "rm^3", numPhases ); + addPhaseColumns( tableLayout, "Phase mass", massUnit, numPhases ); + addPhaseColumns( tableLayout, "Trapped phase mass (metric 1)", massUnit, numPhases ); + addPhaseColumns( tableLayout, "Non-trapped phase mass (metric 1)", massUnit, numPhases ); + addPhaseColumns( tableLayout, "Immobile phase mass (metric 2)", massUnit, numPhases ); + addPhaseColumns( tableLayout, "Mobile phase mass (metric 2)", massUnit, numPhases ); + addPhaseCompColumns( tableLayout, "Component mass", massUnit, numPhases, numComps ); + + auto csvFormatter = std::make_unique< TableCSVFormatter >( tableLayout ); + m_csvFormatters.emplace( meshName, std::move( csvFormatter ) ); + + // output CSV header + std::ofstream outputFile( GEOS_FMT( "{}/{}.csv", m_outputDir, meshName )); + outputFile << csvFormatter->headerToString(); + GEOS_LOG( GEOS_FMT( "table {} : {}", meshName, csvFormatter->headerToString() ) ); // TODO : remove this log +} + +bool StatsTask::execute( real64 const time_n, + real64 const dt, + integer const GEOS_UNUSED_PARAM( cycleNumber ), + integer const GEOS_UNUSED_PARAM( eventCounter ), + real64 const GEOS_UNUSED_PARAM( eventProgress ), + DomainPartition & domain ) +{ + // current time is time_n + dt. TODO: verify implication of events ordering in 'time_n+dt' validity + real64 statsTime = time_n + dt; + + m_solver->forDiscretizationOnMeshTargets( domain.getMeshBodies(), [&] ( string const &, + MeshLevel & mesh, + string_array const & regionNames ) + { + if( m_computeRegionStatistics ) + { + m_aggregator.computeRegionsStatistics( statsTime, mesh, regionNames ); + outputLogStats( statsTime, mesh, regionNames ); + outputCsvStats( statsTime, mesh, regionNames ); + } + } ); + + if( m_computeCFLNumbers ) + m_aggregator.computeCFLNumbers( statsTime, dt, domain ); + + return false; +} + +void StatsTask::outputLogStats( real64 const statsTime, + MeshLevel & mesh, + string_array const & regionNames ) +{ + auto const formatterIter = m_logFormatters.find( mesh.getName()); + if( formatterIter==m_logFormatters.end()) + return; + + TableFormatter const & formatter = *formatterIter->second; + TableData tableData; + static constexpr auto merge = CellType::MergeNext; + + tableData.addRow( "Statistics time", merge, merge, statsTime ); + forRegionStatistics( []( string_view regionName, RegionStatistics const & stats ) + { + string_view massUnit = units::getSymbol( m_solver->getMassUnit() ); + string_view pressureUnit = units::getSymbol( units::Pressure ); + string_view tempUnit = units::getSymbol( units::Temperature ); + string_view resVolUnit = "rm3"; + + tableData.addRow( merge, merge, merge, "" ); + + tableData.addSeparator(); + tableData.addRow( merge, merge, merge, GEOS_FMT( "Region '{}'", regionName ) ); + tableData.addRow( "statistics", "min", "average", "max" ); + tableData.addSeparator(); + + tableData.addRow( GEOS_FMT( "Pressure [{}]", pressureUnit ), + stats.minPressure, stats.averagePressure, stats.maxPressure ); + tableData.addRow( GEOS_FMT( "Delta pressure [{}]", pressureUnit ), + stats.minDeltaPressure, "/", stats.maxDeltaPressure ); + tableData.addRow( GEOS_FMT( "Temperature [{}]", tempUnit ), + stats.minTemperature, stats.averageTemperature, stats.maxTemperature ); + + tableData.addSeparator(); + + tableData.addRow( GEOS_FMT( "Total dynamic pore volume [{}]", resVolUnit ), CellType::MergeNext, CellType::MergeNext, stats.totalPoreVolume ); + tableData.addRow( GEOS_FMT( "Phase dynamic pore volume [{}]", resVolUnit ), + stringutilities::joinLambda( phaseNames, "\n", []( auto data ) { return data[0]; } ), + CellType::MergeNext, + stringutilities::joinLambda( stats.phasePoreVolume, "\n", []( auto data ) { return data[0]; } ) ); + + tableData.addSeparator(); + + tableData.addRow( GEOS_FMT( "Phase mass [{}]", massUnit ), + stringutilities::joinLambda( phaseNames, "\n", []( auto data ) { return data[0]; } ), + CellType::MergeNext, + stringutilities::joinLambda( stats.phaseMass, "\n", []( auto data ) { return data[0]; } ) ); + + tableData.addSeparator(); + + tableData.addRow( GEOS_FMT( "Trapped phase mass (metric 1) [{}]", massUnit ), + stringutilities::joinLambda( phaseNames, "\n", []( auto value ) { return value[0]; } ), + CellType::MergeNext, + stringutilities::joinLambda( stats.trappedPhaseMass, "\n", []( auto value ) { return value[0]; } ) ); + tableData.addRow( GEOS_FMT( "Non-trapped phase mass (metric 1) [{}]", massUnit ), + stringutilities::joinLambda( phaseNames, "\n", []( auto value ) { return value[0]; } ), + CellType::MergeNext, + stringutilities::joinLambda( nonTrappedPhaseMass, "\n", []( auto value ) { return value[0]; } ) ); + + tableData.addSeparator(); + + tableData.addRow( GEOS_FMT( "Immobile phase mass (metric 2) [{}]", massUnit ), + stringutilities::joinLambda( phaseNames, "\n", []( auto value ) { return value[0]; } ), + CellType::MergeNext, + stringutilities::joinLambda( stats.immobilePhaseMass, "\n", []( auto value ) { return value[0]; } ) ); + tableData.addRow( GEOS_FMT( "Mobile phase mass (metric 2) [{}]", massUnit ), + stringutilities::joinLambda( phaseNames, "\n", []( auto value ) { return value[0]; } ), + CellType::MergeNext, + stringutilities::joinLambda( mobilePhaseMass, "\n", []( auto value ) { return value[0]; } ) ); + + tableData.addSeparator(); + + tableData.addRow( GEOS_FMT( "Component mass [{}]", massUnit ), + stringutilities::join( phaseCompName, '\n' ), + CellType::MergeNext, + stringutilities::join( massValues, '\n' ) ); + + tableData.addSeparator(); + } ); +} + +} /* namespace compositionalMultiphaseStatistics */ + +REGISTER_CATALOG_ENTRY( TaskBase, + compositionalMultiphaseStatistics::StatsTask, + string const &, dataRepository::Group * const ) + +} /* namespace geos */ diff --git a/src/coreComponents/physicsSolvers/fluidFlow/CompositionalMultiphaseStatistics.hpp b/src/coreComponents/physicsSolvers/fluidFlow/CompositionalMultiphaseStatisticsTask.hpp similarity index 56% rename from src/coreComponents/physicsSolvers/fluidFlow/CompositionalMultiphaseStatistics.hpp rename to src/coreComponents/physicsSolvers/fluidFlow/CompositionalMultiphaseStatisticsTask.hpp index dbe275b7b6f..78fb07eac3e 100644 --- a/src/coreComponents/physicsSolvers/fluidFlow/CompositionalMultiphaseStatistics.hpp +++ b/src/coreComponents/physicsSolvers/fluidFlow/CompositionalMultiphaseStatisticsTask.hpp @@ -17,22 +17,29 @@ * @file CompositionalMultiphaseStatistics.hpp */ -#ifndef SRC_CORECOMPONENTS_PHYSICSSOLVERS_FLUIDFLOW_COMPOSITIONALMULTIPHASESTATISTICS_HPP_ -#define SRC_CORECOMPONENTS_PHYSICSSOLVERS_FLUIDFLOW_COMPOSITIONALMULTIPHASESTATISTICS_HPP_ +#ifndef SRC_CORECOMPONENTS_PHYSICSSOLVERS_FLUIDFLOW_COMPOSITIONALMULTIPHASESTATISTICSTASK_HPP_ +#define SRC_CORECOMPONENTS_PHYSICSSOLVERS_FLUIDFLOW_COMPOSITIONALMULTIPHASESTATISTICSTASK_HPP_ +#include "common/DataTypes.hpp" +#include "common/format/table/TableFormatter.hpp" +#include "common/format/table/TableLayout.hpp" #include "physicsSolvers/FieldStatisticsBase.hpp" +#include "physicsSolvers/fluidFlow/CompositionalMultiphaseStatisticsAggregator.hpp" namespace geos { class CompositionalMultiphaseBase; +namespace compositionalMultiphaseStatistics +{ + /** - * @class CompositionalMultiphaseStatistics + * @class CompositionalMultiphaseStatistics::Task * * Task class allowing for the computation of aggregate statistics in compositional multiphase simulations */ -class CompositionalMultiphaseStatistics : public FieldStatisticsBase< CompositionalMultiphaseBase > +class StatsTask : public FieldStatisticsBase< CompositionalMultiphaseBase > { public: @@ -41,16 +48,11 @@ class CompositionalMultiphaseStatistics : public FieldStatisticsBase< Compositio * @param[in] name the name of the task coming from the xml * @param[in] parent the parent group of the task */ - CompositionalMultiphaseStatistics( const string & name, - Group * const parent ); + StatsTask( const string & name, dataRepository::Group * const parent ); /// Accessor for the catalog name static string catalogName() { return "CompositionalMultiphaseStatistics"; } - /// Accessor for the region statistics catalog name - static string regionStatisticsName() { return "regionStatistics"; } - - /** * @defgroup Tasks Interface Functions * @@ -67,43 +69,6 @@ class CompositionalMultiphaseStatistics : public FieldStatisticsBase< Compositio /**@}*/ - struct RegionStatistics - { - /// average region pressure - real64 averagePressure; - /// minimum region pressure - real64 minPressure; - /// maximum region pressure - real64 maxPressure; - - /// minimum region delta pressure - real64 minDeltaPressure; - /// maximum region delta pressure - real64 maxDeltaPressure; - - /// average region temperature - real64 averageTemperature; - /// minimum region temperature - real64 minTemperature; - /// maximum region temperature - real64 maxTemperature; - - /// total region pore volume - real64 totalPoreVolume; - /// total region uncompacted pore volume - real64 totalUncompactedPoreVolume; - /// phase region phase pore volume - array1d< real64 > phasePoreVolume; - - /// region phase mass (trapped and non-trapped, immobile and mobile) - array1d< real64 > phaseMass; - /// trapped region phase mass - array1d< real64 > trappedPhaseMass; - /// immobile region phase mass - array1d< real64 > immobilePhaseMass; - /// region component mass - array2d< real64 > componentMass; - }; private: using Base = FieldStatisticsBase< CompositionalMultiphaseBase >; @@ -117,37 +82,36 @@ class CompositionalMultiphaseStatistics : public FieldStatisticsBase< Compositio constexpr static char const * computeCFLNumbersString() { return "computeCFLNumbers"; } /// String for the flag deciding the computation of the region statistics constexpr static char const * computeRegionStatisticsString() { return "computeRegionStatistics"; } - /// String for the region statistics - constexpr static char const * regionStatisticsString() { return "regionStatistics"; } /// String for the relperm threshold constexpr static char const * relpermThresholdString() { return "relpermThreshold"; } }; + void postInputInitialization() override; - /** - * @brief Compute some statistics on the reservoir (average field pressure, etc) - * @param[in] time current time - * @param[in] mesh the mesh level object - * @param[in] regionNames the array of target region names - */ - void computeRegionStatistics( real64 const time, - MeshLevel & mesh, - string_array const & regionNames ) const; + void registerDataOnMesh( Group & meshBodies ) override; - /** - * @brief Compute CFL numbers - * @param[in] time current time - * @param[in] dt the time step size - * @param[in] domain the domain partition - */ - void computeCFLNumbers( real64 const time, - real64 const dt, - DomainPartition & domain ) const; + void prepareFluidMetaData(); - void postInputInitialization() override; + void prepareLogTableLayouts( string_view tableName ); - void registerDataOnMesh( Group & meshBodies ) override; + void prepareCsvTableLayouts( string_view tableName ); + + void outputLogStats( real64 statsTime, + MeshLevel & mesh, + string_array const & regionNames ); + + void outputCsvStats( real64 statsTime, + MeshLevel & mesh, + string_array const & regionNames ); + + /// For each discretization (MeshLevel name), table formatter for log output. + stdMap< string, std::unique_ptr< TableFormatter > > m_logFormatters; + + /// For each discretization (MeshLevel name), table formatter for csv output. + stdMap< string, std::unique_ptr< TableFormatter > > m_csvFormatters; + + StatsAggregator m_aggregator; /// Flag to decide whether CFL numbers are computed or not integer m_computeCFLNumbers; @@ -158,9 +122,19 @@ class CompositionalMultiphaseStatistics : public FieldStatisticsBase< Compositio /// Threshold to decide whether a phase is considered "mobile" or not real64 m_relpermThreshold; + struct FluidMetaData + { + integer m_numPhases; + integer m_numComps; + stdVector< string > m_phaseNames; + stdVector< string > m_compNames; + stdVector< stdVector< string > > m_phaseCompNames; + } m_fluid; + }; +} /* namespace compositionalMultiphaseStatistics */ } /* namespace geos */ -#endif /* SRC_CORECOMPONENTS_PHYSICSSOLVERS_FLUIDFLOW_COMPOSITIONALMULTIPHASESTATISTICS_HPP_ */ +#endif /* SRC_CORECOMPONENTS_PHYSICSSOLVERS_FLUIDFLOW_COMPOSITIONALMULTIPHASESTATISTICSTASK_HPP_ */ From 37d75e5a193dd5a6d16ac6f44d186ede1c903777 Mon Sep 17 00:00:00 2001 From: MelReyCG Date: Tue, 2 Dec 2025 10:54:12 +0100 Subject: [PATCH 5/6] =?UTF-8?q?=F0=9F=9A=A7wip=202?= MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit --- ...sitionalMultiphaseStatisticsAggregator.cpp | 32 +++++++-------- ...sitionalMultiphaseStatisticsAggregator.hpp | 34 ++++++++++----- .../CompositionalMultiphaseStatisticsTask.cpp | 41 ++++++++++--------- 3 files changed, 62 insertions(+), 45 deletions(-) diff --git a/src/coreComponents/physicsSolvers/fluidFlow/CompositionalMultiphaseStatisticsAggregator.cpp b/src/coreComponents/physicsSolvers/fluidFlow/CompositionalMultiphaseStatisticsAggregator.cpp index 13d4fc78b9c..e459eb9b10d 100644 --- a/src/coreComponents/physicsSolvers/fluidFlow/CompositionalMultiphaseStatisticsAggregator.cpp +++ b/src/coreComponents/physicsSolvers/fluidFlow/CompositionalMultiphaseStatisticsAggregator.cpp @@ -46,7 +46,10 @@ StatsAggregator::StatsAggregator(): void StatsAggregator::enableRegionStatistics( dataRepository::Group & meshBodies ) { - GEOS_ASSERT_NE( m_solver, nullptr ); + GEOS_ERROR_IF_EQ_MSG( m_solver, nullptr, "Flow solver must be set." ); + + integer const numPhases = m_solver->numFluidPhases(); + integer const numComps = m_solver->numFluidComponents(); m_solver->forDiscretizationOnMeshTargets( meshBodies, [&] ( string const &, MeshLevel & mesh, @@ -54,9 +57,6 @@ void StatsAggregator::enableRegionStatistics( dataRepository::Group & meshBodies { ElementRegionManager & elemManager = mesh.getElemManager(); - integer const numPhases = m_solver->numFluidPhases(); - integer const numComps = m_solver->numFluidComponents(); - for( size_t i = 0; i < regionNames.size(); ++i ) { ElementRegionBase & region = elemManager.getRegion( regionNames[i] ); @@ -69,7 +69,9 @@ void StatsAggregator::enableRegionStatistics( dataRepository::Group & meshBodies stats.phasePoreVolume.resizeDimension< 0 >( numPhases ); stats.phaseMass.resizeDimension< 0 >( numPhases ); stats.trappedPhaseMass.resizeDimension< 0 >( numPhases ); + stats.nonTrappedPhaseMass.resizeDimension< 0 >( numPhases ); stats.immobilePhaseMass.resizeDimension< 0 >( numPhases ); + stats.mobilePhaseMass.resizeDimension< 0 >( numPhases ); stats.componentMass.resizeDimension< 0, 1 >( numPhases, numComps ); } } ); @@ -78,15 +80,15 @@ void StatsAggregator::enableRegionStatistics( dataRepository::Group & meshBodies void StatsAggregator::enableCFLStatistics( dataRepository::Group & meshBodies ) { - GEOS_ASSERT_NE( m_solver, nullptr ); + GEOS_ERROR_IF_EQ_MSG( m_solver, nullptr, "Flow solver must be set." ); m_solver->registerDataForCFL( meshBodies ); m_isCFLNumberEnabled = true; } -void StatsAggregator::computeDiscretizationStatistics( real64 const time, - MeshLevel & mesh, - string_array const & regionNames ) const +bool StatsAggregator::computeRegionsStatistics( real64 const time, + MeshLevel & mesh, + string_array const & regionNames ) const { GEOS_MARK_FUNCTION; @@ -117,7 +119,9 @@ void StatsAggregator::computeDiscretizationStatistics( real64 const time, stats.phaseMass.setValues< serialPolicy >( 0.0 ); stats.trappedPhaseMass.setValues< serialPolicy >( 0.0 ); + stats.nonTrappedPhaseMass.setValues< serialPolicy >( 0.0 ); stats.immobilePhaseMass.setValues< serialPolicy >( 0.0 ); + stats.mobilePhaseMass.setValues< serialPolicy >( 0.0 ); stats.componentMass.setValues< serialPolicy >( 0.0 ); } @@ -289,18 +293,14 @@ void StatsAggregator::computeDiscretizationStatistics( real64 const time, stats.averagePressure = 0.0; stats.averageTemperature = 0.0; GEOS_LOG_LEVEL_RANK_0( logInfo::Statistics, - GEOS_FMT( "{}, {}: Cannot compute average pressure because region pore volume is zero.", - getName(), regionNames[i] ) ); + GEOS_FMT( "Cannot compute average pressure because region pore volume is zero in region '{}'.", + regionNames[i] ) ); } - - // helpers to report statistics - array1d< real64 > nonTrappedPhaseMass( numPhases ); - array1d< real64 > mobilePhaseMass( numPhases ); for( integer ip = 0; ip < numPhases; ++ip ) { - nonTrappedPhaseMass[ip] = stats.phaseMass[ip] - stats.trappedPhaseMass[ip]; - mobilePhaseMass[ip] = stats.phaseMass[ip] - stats.immobilePhaseMass[ip]; + stats.nonTrappedPhaseMass[ip] = stats.phaseMass[ip] - stats.trappedPhaseMass[ip]; + stats.mobilePhaseMass[ip] = stats.phaseMass[ip] - stats.immobilePhaseMass[ip]; } string_view massUnit = units::getSymbol( m_solver->getMassUnit() ); diff --git a/src/coreComponents/physicsSolvers/fluidFlow/CompositionalMultiphaseStatisticsAggregator.hpp b/src/coreComponents/physicsSolvers/fluidFlow/CompositionalMultiphaseStatisticsAggregator.hpp index 419a7e6d35f..3735f197efe 100644 --- a/src/coreComponents/physicsSolvers/fluidFlow/CompositionalMultiphaseStatisticsAggregator.hpp +++ b/src/coreComponents/physicsSolvers/fluidFlow/CompositionalMultiphaseStatisticsAggregator.hpp @@ -20,6 +20,7 @@ #ifndef SRC_CORECOMPONENTS_PHYSICSSOLVERS_FLUIDFLOW_COMPOSITIONALMULTIPHASESTATISTICSAGGREGATOR_HPP_ #define SRC_CORECOMPONENTS_PHYSICSSOLVERS_FLUIDFLOW_COMPOSITIONALMULTIPHASESTATISTICSAGGREGATOR_HPP_ +#include "common/StdContainerWrappers.hpp" #include "dataRepository/Group.hpp" #include "mesh/DomainPartition.hpp" #include "mesh/MeshLevel.hpp" @@ -75,8 +76,12 @@ struct RegionStatistics array1d< real64 > phaseMass; /// trapped region phase mass array1d< real64 > trappedPhaseMass; + /// non-trapped region phase mass + array1d< real64 > nonTrappedPhaseMass; /// immobile region phase mass array1d< real64 > immobilePhaseMass; + /// mobile region phase mass + array1d< real64 > mobilePhaseMass; /// region component mass array2d< real64 > componentMass; @@ -107,7 +112,7 @@ class StatsAggregator constexpr static char const * regionStatisticsString() { return "regionStatistics"; } }; - using RegionFunctor = std::function< void (string, RegionStatistics const &) >; + using RegionFunctor = std::function< void (string_view, RegionStatistics const &) >; StatsAggregator(); @@ -124,14 +129,14 @@ class StatsAggregator // std::function< void(MeshLevel const &, // string_array const & regionNames) > functor ) const; - // void forRegionStatistics( DomainPartition const &, - // MeshLevel const & mesh, - // RegionFunctor functor ) const; + void forRegionStatistics( DomainPartition const &, + MeshLevel const & mesh, + RegionFunctor functor ) const; /** - * @brief Register the results structs & wrappers so they will be targeted by TimeHistory output + * @brief Enable the computation of region statistics, initialize data structure to collect them. + * Register the resulting data wrappers so they will be targeted by TimeHistory output * @note Must be called in "registerDataOnMesh" initialization phase - * @param solver The flow solver * @param meshBodies The Group containing the MeshBody objects */ void enableRegionStatistics( dataRepository::Group & meshBodies ); @@ -139,20 +144,21 @@ class StatsAggregator /** * @brief Register the results structs & wrappers so they will be targeted by TimeHistory output * @note Must be called in "registerDataOnMesh" initialization phase - * @param solver The flow solver * @param meshBodies The Group containing the MeshBody objects */ void enableCFLStatistics( dataRepository::Group & meshBodies ); /** * @brief Compute some statistics on a given mesh discretization (average field pressure, etc) + * Results are reduced on rank 0, and broadcasted over all ranks. * @param[in] time current time * @param[in] mesh the mesh level object * @param[in] regionNames the array of target region names + * return false if */ - void computeDiscretizationStatistics( real64 const time, - MeshLevel & mesh, - string_array const & regionNames ) const; + bool computeRegionsStatistics( real64 const time, + MeshLevel & mesh, + string_array const & regionNames ) const; /** * @brief Compute CFL numbers @@ -167,6 +173,12 @@ class StatsAggregator CFLStatistics const & getCFLStatistics() const { return m_cflStats; } + CompositionalMultiphaseBase const * getSolver() const + { return m_solver; } + + stdVector< string_view > const & getIssues() const + { return m_issues; } + private: CompositionalMultiphaseBase * m_solver; @@ -175,6 +187,8 @@ class StatsAggregator CFLStatistics m_cflStats; + stdVector< string_view > m_issues; + bool m_isRegionStatsEnabled; bool m_isCFLNumberEnabled; diff --git a/src/coreComponents/physicsSolvers/fluidFlow/CompositionalMultiphaseStatisticsTask.cpp b/src/coreComponents/physicsSolvers/fluidFlow/CompositionalMultiphaseStatisticsTask.cpp index b3d6a6b59b8..608e57fe9b9 100644 --- a/src/coreComponents/physicsSolvers/fluidFlow/CompositionalMultiphaseStatisticsTask.cpp +++ b/src/coreComponents/physicsSolvers/fluidFlow/CompositionalMultiphaseStatisticsTask.cpp @@ -27,6 +27,7 @@ #include "constitutive/relativePermeability/RelativePermeabilityBase.hpp" #include "constitutive/solid/CoupledSolidBase.hpp" #include "physicsSolvers/LogLevelsInfo.hpp" +#include "physicsSolvers/fluidFlow/LogLevelsInfo.hpp" #include "physicsSolvers/fluidFlow/CompositionalMultiphaseBase.hpp" #include "physicsSolvers/fluidFlow/CompositionalMultiphaseBaseFields.hpp" #include "physicsSolvers/fluidFlow/CompositionalMultiphaseHybridFVM.hpp" @@ -233,7 +234,7 @@ void StatsTask::outputLogStats( real64 const statsTime, MeshLevel & mesh, string_array const & regionNames ) { - auto const formatterIter = m_logFormatters.find( mesh.getName()); + auto const formatterIter = m_logFormatters.find( mesh.getName() ); if( formatterIter==m_logFormatters.end()) return; @@ -241,13 +242,15 @@ void StatsTask::outputLogStats( real64 const statsTime, TableData tableData; static constexpr auto merge = CellType::MergeNext; + string_view massUnit = units::getSymbol( m_aggregator.getSolver()->getMassUnit() ); + string_view pressureUnit = units::getSymbol( units::Pressure ); + string_view tempUnit = units::getSymbol( units::Temperature ); + string_view resVolUnit = units::getSymbol( units::ReservoirVolume ); + tableData.addRow( "Statistics time", merge, merge, statsTime ); - forRegionStatistics( []( string_view regionName, RegionStatistics const & stats ) + m_aggregator.forRegionStatistics( MeshLevel & mesh, + [&]( string_view regionName, RegionStatistics const & stats ) { - string_view massUnit = units::getSymbol( m_solver->getMassUnit() ); - string_view pressureUnit = units::getSymbol( units::Pressure ); - string_view tempUnit = units::getSymbol( units::Temperature ); - string_view resVolUnit = "rm3"; tableData.addRow( merge, merge, merge, "" ); @@ -267,54 +270,54 @@ void StatsTask::outputLogStats( real64 const statsTime, tableData.addRow( GEOS_FMT( "Total dynamic pore volume [{}]", resVolUnit ), CellType::MergeNext, CellType::MergeNext, stats.totalPoreVolume ); tableData.addRow( GEOS_FMT( "Phase dynamic pore volume [{}]", resVolUnit ), - stringutilities::joinLambda( phaseNames, "\n", []( auto data ) { return data[0]; } ), + stringutilities::joinLambda( m_fluid.m_phaseNames, "\n", []( auto data ) { return data[0]; } ), CellType::MergeNext, stringutilities::joinLambda( stats.phasePoreVolume, "\n", []( auto data ) { return data[0]; } ) ); tableData.addSeparator(); tableData.addRow( GEOS_FMT( "Phase mass [{}]", massUnit ), - stringutilities::joinLambda( phaseNames, "\n", []( auto data ) { return data[0]; } ), + stringutilities::joinLambda( m_fluid.m_phaseNames, "\n", []( auto data ) { return data[0]; } ), CellType::MergeNext, stringutilities::joinLambda( stats.phaseMass, "\n", []( auto data ) { return data[0]; } ) ); tableData.addSeparator(); tableData.addRow( GEOS_FMT( "Trapped phase mass (metric 1) [{}]", massUnit ), - stringutilities::joinLambda( phaseNames, "\n", []( auto value ) { return value[0]; } ), + stringutilities::joinLambda( m_fluid.m_phaseNames, "\n", []( auto value ) { return value[0]; } ), CellType::MergeNext, stringutilities::joinLambda( stats.trappedPhaseMass, "\n", []( auto value ) { return value[0]; } ) ); tableData.addRow( GEOS_FMT( "Non-trapped phase mass (metric 1) [{}]", massUnit ), - stringutilities::joinLambda( phaseNames, "\n", []( auto value ) { return value[0]; } ), + stringutilities::joinLambda( m_fluid.m_phaseNames, "\n", []( auto value ) { return value[0]; } ), CellType::MergeNext, - stringutilities::joinLambda( nonTrappedPhaseMass, "\n", []( auto value ) { return value[0]; } ) ); + stringutilities::joinLambda( stats.nonTrappedPhaseMass, "\n", []( auto value ) { return value[0]; } ) ); tableData.addSeparator(); tableData.addRow( GEOS_FMT( "Immobile phase mass (metric 2) [{}]", massUnit ), - stringutilities::joinLambda( phaseNames, "\n", []( auto value ) { return value[0]; } ), + stringutilities::joinLambda( m_fluid.m_phaseNames, "\n", []( auto value ) { return value[0]; } ), CellType::MergeNext, stringutilities::joinLambda( stats.immobilePhaseMass, "\n", []( auto value ) { return value[0]; } ) ); tableData.addRow( GEOS_FMT( "Mobile phase mass (metric 2) [{}]", massUnit ), - stringutilities::joinLambda( phaseNames, "\n", []( auto value ) { return value[0]; } ), + stringutilities::joinLambda( m_fluid.m_phaseNames, "\n", []( auto value ) { return value[0]; } ), CellType::MergeNext, - stringutilities::joinLambda( mobilePhaseMass, "\n", []( auto value ) { return value[0]; } ) ); + stringutilities::joinLambda( stats.mobilePhaseMass, "\n", []( auto value ) { return value[0]; } ) ); tableData.addSeparator(); tableData.addRow( GEOS_FMT( "Component mass [{}]", massUnit ), - stringutilities::join( phaseCompName, '\n' ), + stringutilities::join( m_fluid.m_phaseCompNames, '\n' ), CellType::MergeNext, - stringutilities::join( massValues, '\n' ) ); + stringutilities::join( stats.componentMass, '\n' ) ); tableData.addSeparator(); } ); } -} /* namespace compositionalMultiphaseStatistics */ - REGISTER_CATALOG_ENTRY( TaskBase, - compositionalMultiphaseStatistics::StatsTask, + StatsTask, string const &, dataRepository::Group * const ) +} /* namespace compositionalMultiphaseStatistics */ + } /* namespace geos */ From 239ff243eabac0a3f04abba2757e811d983b81da Mon Sep 17 00:00:00 2001 From: MelReyCG Date: Tue, 2 Dec 2025 11:47:13 +0100 Subject: [PATCH 6/6] =?UTF-8?q?=F0=9F=9A=A7=20applying=20merge=20change=20?= =?UTF-8?q?to=20task=20component?= MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit --- .../CompositionalMultiphaseStatisticsTask.cpp | 26 +++++++++---------- 1 file changed, 13 insertions(+), 13 deletions(-) diff --git a/src/coreComponents/physicsSolvers/fluidFlow/CompositionalMultiphaseStatisticsTask.cpp b/src/coreComponents/physicsSolvers/fluidFlow/CompositionalMultiphaseStatisticsTask.cpp index 608e57fe9b9..bf340da9478 100644 --- a/src/coreComponents/physicsSolvers/fluidFlow/CompositionalMultiphaseStatisticsTask.cpp +++ b/src/coreComponents/physicsSolvers/fluidFlow/CompositionalMultiphaseStatisticsTask.cpp @@ -173,19 +173,19 @@ void StatsTask::prepareCsvTableLayouts( string_view meshName ) string_view massUnit = units::getSymbol( m_solver->getMassUnit() ); TableLayout tableLayout( { - TableLayout::Column().setName( "Time [s]" ), - TableLayout::Column().setName( "Region" ), // TODO : mention this change in PR description - TableLayout::Column().setName( "Min pressure [Pa]" ), - TableLayout::Column().setName( "Average pressure [Pa]" ), - TableLayout::Column().setName( "Max pressure [Pa]" ), - TableLayout::Column().setName( "Min delta pressure [Pa]" ), - TableLayout::Column().setName( "Max delta pressure [Pa]" ), - TableLayout::Column().setName( "Min temperature [Pa]" ), - TableLayout::Column().setName( "Average temperature [Pa]" ), - TableLayout::Column().setName( "Max temperature [Pa]" ), - TableLayout::Column().setName( "Total dynamic pore volume [rm^3]" ), - } ); - addPhaseColumns( tableLayout, "Phase dynamic pore volume", "rm^3", numPhases ); + TableLayout::Column().setName( GEOS_FMT( "Time [{}]", units::getSymbol( units::Unit::Time ))), + TableLayout::Column().setName( "Region" ), // TODO : mention this change in PR description + TableLayout::Column().setName( GEOS_FMT( "Min pressure [{}]", units::getSymbol( units::Unit::Pressure ))), + TableLayout::Column().setName( GEOS_FMT( "Average pressure [{}]", units::getSymbol( units::Unit::Pressure )) ), + TableLayout::Column().setName( GEOS_FMT( "Max pressure [{}]", units::getSymbol( units::Unit::Pressure ) ) ), + TableLayout::Column().setName( GEOS_FMT( "Min delta pressure [{}]", units::getSymbol( units::Unit::Pressure ))), + TableLayout::Column().setName( GEOS_FMT( "Max delta pressure [{}]", units::getSymbol( units::Unit::Pressure ))), + TableLayout::Column().setName( GEOS_FMT( "Min temperature [{}]", units::getSymbol( units::Unit::Temperature ) )), + TableLayout::Column().setName( GEOS_FMT( "Average temperature [{}]", units::getSymbol( units::Unit::Temperature ) )), + TableLayout::Column().setName( GEOS_FMT( "Max temperature [{}]", units::getSymbol( units::Unit::Temperature ) )), + TableLayout::Column().setName( GEOS_FMT( "Total dynamic pore volume [{}]", units::getSymbol( units::Unit::ReservoirVolume ) )), + } ); + addPhaseColumns( tableLayout, "Phase dynamic pore volume", units::getSymbol( units::Unit::ReservoirVolume ) ), numPhases ); addPhaseColumns( tableLayout, "Phase mass", massUnit, numPhases ); addPhaseColumns( tableLayout, "Trapped phase mass (metric 1)", massUnit, numPhases ); addPhaseColumns( tableLayout, "Non-trapped phase mass (metric 1)", massUnit, numPhases );