diff --git a/src/coreComponents/common/format/table/TableFormatter.hpp b/src/coreComponents/common/format/table/TableFormatter.hpp index c2035bc2a1c..e73927a9b2d 100644 --- a/src/coreComponents/common/format/table/TableFormatter.hpp +++ b/src/coreComponents/common/format/table/TableFormatter.hpp @@ -57,6 +57,12 @@ class TableFormatter TableErrorListing & getErrorsList() const { return *m_errors; } + /** + * @return The prepared table layout + */ + PreparedTableLayout const & getLayout() const + { return m_tableLayout; } + protected: /// Layout for a table diff --git a/src/coreComponents/common/format/table/TableTypes.hpp b/src/coreComponents/common/format/table/TableTypes.hpp index 59826ddde6a..11e2f86ef65 100644 --- a/src/coreComponents/common/format/table/TableTypes.hpp +++ b/src/coreComponents/common/format/table/TableTypes.hpp @@ -82,7 +82,7 @@ class TableErrorListing * @brief Append a vector of string to the errors vector. * @param errors A vector of string to append */ - void appendErrors( std::vector< string > & errors ) + void appendErrors( std::vector< string > const & errors ) { m_errorList.insert( m_errorList.end(), errors.begin(), errors.end() );} /** diff --git a/src/coreComponents/dataRepository/DataContext.hpp b/src/coreComponents/dataRepository/DataContext.hpp index 8a43381b438..f34ca3c8687 100644 --- a/src/coreComponents/dataRepository/DataContext.hpp +++ b/src/coreComponents/dataRepository/DataContext.hpp @@ -78,7 +78,7 @@ class DataContext /** * @return Get the target object name */ - string getTargetName() const + string const & getTargetName() const { return m_targetName; } /** * @brief Insert contextual information in the provided stream. @@ -179,13 +179,13 @@ class DataFileContext final : public DataContext /** * @return the type name in the source file (XML node tag name / attribute name). */ - string getTypeName() const + string const & getTypeName() const { return m_typeName; } /** * @return the source file path where the target object has been declared. */ - string getFilePath() const + string const & getFilePath() const { return m_filePath; } /** diff --git a/src/coreComponents/mesh/MeshBody.hpp b/src/coreComponents/mesh/MeshBody.hpp index c0f30715c0c..a8f8cbe70af 100644 --- a/src/coreComponents/mesh/MeshBody.hpp +++ b/src/coreComponents/mesh/MeshBody.hpp @@ -96,8 +96,7 @@ class MeshBody : public dataRepository::Group * @param[in] level The lookup key of the MeshLevel * @return const reference to the MeshLevel */ - template< typename T, std::enable_if_t< std::is_same< T, string >::value || - std::is_same< T, const char * >::value, bool > = false > + template< typename T > MeshLevel & getMeshLevel( T const & level ) const { return m_meshLevels.getGroup< MeshLevel >( level ); } @@ -108,8 +107,7 @@ class MeshBody : public dataRepository::Group * @param[in] level The lookup key of the MeshLevel * @return Reference to the MeshLevel */ - template< typename T, std::enable_if_t< std::is_same< T, string >::value || - std::is_same< T, const char * >::value, bool > = false > + template< typename T > MeshLevel & getMeshLevel( T const & level ) { return m_meshLevels.getGroup< MeshLevel >( level ); } diff --git a/src/coreComponents/physicsSolvers/CMakeLists.txt b/src/coreComponents/physicsSolvers/CMakeLists.txt index 246fc46be5f..ffce065fbbf 100644 --- a/src/coreComponents/physicsSolvers/CMakeLists.txt +++ b/src/coreComponents/physicsSolvers/CMakeLists.txt @@ -29,6 +29,8 @@ set( physicsSolversBase_headers PhysicsSolverBase.hpp PhysicsSolverBaseKernels.hpp SolverStatistics.hpp + StatisticsAggregatorBase.hpp + StatisticsAggregatorBaseHelpers.hpp FieldStatisticsBase.hpp ) # # Specify solver sources diff --git a/src/coreComponents/physicsSolvers/FieldStatisticsBase.hpp b/src/coreComponents/physicsSolvers/FieldStatisticsBase.hpp index 13390fe314c..1548ed4b6e6 100644 --- a/src/coreComponents/physicsSolvers/FieldStatisticsBase.hpp +++ b/src/coreComponents/physicsSolvers/FieldStatisticsBase.hpp @@ -50,8 +50,7 @@ class FieldStatisticsBase : public TaskBase m_outputDir( joinPath( OutputBase::getOutputDirectory(), name ) ) { - string const key = SOLVER::coupledSolverAttributePrefix() + "SolverName"; - registerWrapper( key, &m_solverName ). + registerWrapper( getSolverWrapperKey(), &m_solverName ). setRTTypeName( rtTypes::CustomTypes::groupNameRef ). setInputFlag( dataRepository::InputFlags::REQUIRED ). setDescription( "Name of the " + SOLVER::coupledSolverAttributePrefix() + " solver" ); @@ -80,6 +79,20 @@ class FieldStatisticsBase : public TaskBase protected: + struct viewKeyStruct + { + static constexpr char const * writeCSVFlagString() { return "writeCSV"; } + }; + + /// Pointer to the physics solver + SOLVER * m_solver; + + // Output directory + string const m_outputDir; + + // Flag to enable writing CSV output + integer m_writeCSV; + void postInputInitialization() override { Group & problemManager = this->getGroupByPath( "/Problem" ); @@ -103,19 +116,8 @@ class FieldStatisticsBase : public TaskBase } } - struct viewKeyStruct - { - static constexpr char const * writeCSVFlagString() { return "writeCSV"; } - }; - - /// Pointer to the physics solver - SOLVER * m_solver; - - // Output directory - string const m_outputDir; - - // Flag to enable writing CSV output - integer m_writeCSV; + string getSolverWrapperKey() const + { return SOLVER::coupledSolverAttributePrefix() + "SolverName"; } private: diff --git a/src/coreComponents/physicsSolvers/StatisticsAggregatorBase.hpp b/src/coreComponents/physicsSolvers/StatisticsAggregatorBase.hpp new file mode 100644 index 00000000000..07392c10402 --- /dev/null +++ b/src/coreComponents/physicsSolvers/StatisticsAggregatorBase.hpp @@ -0,0 +1,268 @@ +/* + * ------------------------------------------------------------------------------------------------------------ + * 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 CompositionalMultiphaseStatisticsAggregator.hpp + */ + +#ifndef SRC_CORECOMPONENTS_PHYSICSSOLVERS_STATISTICSAGGREGATOR_HPP_ +#define SRC_CORECOMPONENTS_PHYSICSSOLVERS_STATISTICSAGGREGATOR_HPP_ + +#include "common/DataTypes.hpp" +#include "dataRepository/DataContext.hpp" +#include "dataRepository/Group.hpp" +#include "mesh/CellElementRegion.hpp" +#include "mesh/MeshBody.hpp" +#include "mesh/MeshLevel.hpp" + +namespace geos +{ + +/** + * @brief Output data group to contain the result of a given stat aggregator on the dataRepository. + * Attributes are public since the class is a POD. + * @todo repair 1D HDF5 outputs to enable stats HDF5 outputs + */ +class RegionStatisticsBase : public dataRepository::Group +{ +public: + + /// Time of statistics computation + real64 m_time; + + /** + * @brief Construct a new Region Statistics Base object + * @param targetName name of the data-repository object that is targeted by the statistics + * (mesh level / region / sub-region). + * @param parent the instance parent in data-repository + */ + RegionStatisticsBase( string const & targetName, dataRepository::Group * const parent ); + + /** + * @return the name of the data-repository object that is targeted by the statistics + * (mesh level / region / sub-region). + */ + string_view getTargetName() const + { return getName(); } + +}; + +template< typename T > +struct StatsAggregatorTraits; + +/** + * @brief Reponsible of computing physical statistics over the grid, registering the result in the + * data repository, but not storing / outputing it by itself. It does not have mutable state + * except the encountered issues. + * @todo repair 1D HDF5 outputs to enable stats HDF5 outputs + * @tparam Impl the derived type of the statistics aggregator which contains all necessary implementations (CRTP) + */ +template< typename Impl > +class StatsAggregatorBase +{ +public: + + using SolverType = typename StatsAggregatorTraits< Impl >::SolverType; + + using StatsGroupType = typename StatsAggregatorTraits< Impl >::StatsGroupType; + + /** + * @brief Standard function signature for any functor that applies on statistics group instances (StatsGroupType) + * - param 0: OwnerType &, the group instance containing the data for which we want to aggregate the statistics (MeshLevel, + * CellElementRegion...) + * - param 1: StatsAggregateGroupType &, the statistics aggregate Group where to store the data + * @tparam OwnerType the concrete type of the OwnerType param + */ + template< typename OwnerType > + using RegionStatsFunc = std::function< void ( OwnerType &, + StatsGroupType & ) >; + + /** + * @brief TODO Document + */ + using RegionStatsRegisterFunc = std::function< StatsGroupType & ( dataRepository::Group &, + string const & ) >; + + /** + * @brief the associated view keys + */ + struct ViewKeys + { + /// String for the discretization statistics group + constexpr static char const * statisticsString() { return "statistics"; } + /// String for the region statistics group + constexpr static char const * regionsStatisticsString() { return "regionsStatistics"; } + }; + + /** + * @brief Construct a new Stats Aggregator object + * @param ownerName the unique name of the entity requesting the statistics. + * An error is thrown if not unique in this context. + */ + StatsAggregatorBase( dataRepository::DataContext const & ownerDataContext ); + + /** + * @brief Enable the computation of any statistics, initialize data structure to collect them. + * Register the resulting data wrappers so they will be targeted by TimeHistory output + * @param solver flow solver object to retrieve: + - the simulated regions, + - fields for statistics computation. + * @param meshBodies The Group containing the MeshBody objects + */ + void initStatisticsAggregation( dataRepository::Group & meshBodies, + SolverType & solver ); + + /** + * @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 or after the "registerDataOnMesh" initialization phase + * @param registerStatsFunc The functor which register each statistics group whithin the regions hierarchy + */ + void enableRegionStatisticsAggregation( RegionStatsRegisterFunc && registerStatsFunc ); + + void forRegionStatistics( RegionStatsFunc< MeshLevel > const & functor ) const; + + void forRegionStatistics( MeshLevel & mesh, + StatsGroupType & meshRegionsStatistics, + RegionStatsFunc< CellElementRegion > const & functor ) const; + + void forRegionStatistics( CellElementRegion & region, + StatsGroupType & regionStatistics, + RegionStatsFunc< CellElementSubRegion > const & functor ) const; + + /** + * @param[in] timeRequest The time for which we want to know if the statistics are computed. + * @param[in] stats the statistics data structure we want to know if it has been computed + * @return true if the statistics have been computed. + */ + bool isComputed( real64 const timeRequest, StatsGroupType const & stats ); + + /** + * @brief set the statistics as dirty, ensuring isComputed() will be false until the next computation. + */ + void setDirty(); + + /** + * @brief Compute statistics on the mesh discretizations (average field pressure, etc) + * Results are reduced on rank 0, and broadcasted over all ranks. + * @param[in] timeRequest The time for which we want to compute the statistics. + * @return false if there was a problem that prevented the statistics to be computed correctly. + */ + bool computeRegionsStatistics( real64 const timeRequest ); + + /** + * @return the name of the entity that needs the statistics. + */ + string const & getOwnerName() const + { return m_ownerDataContext.getTargetName(); } + + /** + * @return The encountered issues during the last computing method call. + */ + stdVector< string > const & getWarnings() const + { return m_warnings; } + + dataRepository::Group & getInstanceStatisticsGroup( MeshLevel & mesh ) const; + + StatsGroupType & getRegionsStatistics( MeshLevel & mesh ) const; + + /** + * @brief TODO + * @throw InputError if no statistics data is found for the given region name. + * @param mesh TODO + * @param regionNname TODO + * @return TODO + */ + StatsGroupType & getRegionStatistics( MeshLevel & mesh, string_view regionName ) const; + +protected: + + struct StatsState + { + bool m_isEnabled = false; + bool m_isDirty = false; + }; + + struct DiscretizationGroupPath + { + localIndex meshBody; + localIndex meshLevel; + string_array regionNames; + }; + + /// @see getOwnerName() + dataRepository::DataContext const & m_ownerDataContext; + + dataRepository::Group * m_meshBodies = nullptr; + + std::vector< DiscretizationGroupPath > m_discretizationsPaths; + + /// @see getWarnings() + stdVector< string > m_warnings; + + /// The current state of the region statistics + StatsState m_regionStatsState; + + MeshLevel & getMeshLevel( DiscretizationGroupPath const & path ) const; + + /** + * @brief Initialize all statistics values to aggregable default values, + * before any computation / reduction for the current timestep. + * @param stats the statistics instance + * @param time start time of the current timestep (s) + * @note Must be implemented for each type that implements this template (CRTP). + */ + void initStats( StatsGroupType & stats, real64 time ) const + { static_cast< Impl const * >(this)->initStats( stats, time ); } + + /** + * @brief Compute the rank-local stats for the given sub-region and store the results in the given stats group. + * @param subRegion + * @param subRegionStats the stats group instance for the subregion + * @note Must be implemented for each type that implements this template (CRTP). + */ + void computeSubRegionRankStats( CellElementSubRegion & subRegion, StatsGroupType & subRegionStats ) const + { static_cast< Impl const * >(this)->computeSubRegionRankStats( subRegion, subRegionStats ); } + + /** + * @brief Aggregate all instance statistics with those of another instance on the current rank. + * @param stats the statistics instance + * @param other the other instance to aggregate with. + * @note Must be implemented for each type that implements this template (CRTP). + */ + void aggregateStats( StatsGroupType & stats, StatsGroupType const & other ) const + { static_cast< Impl const * >(this)->aggregateStats( stats, other ); } + + /** + * @brief Aggregate all instance statistics with those of other ranks. + * @param stats the statistics instance + * @note Must be implemented for each type that implements this template (CRTP). + */ + void mpiAggregateStats( StatsGroupType & stats ) const + { static_cast< Impl const * >(this)->mpiAggregateStats( stats ); } + + /** + * @brief Do the final computations for the statistics. Must be called after computations & aggregations. + * @param stats the statistics instance + * @note Must be implemented for each type that implements this template (CRTP). + */ + void postAggregateStats( StatsGroupType & stats ) + { static_cast< Impl * >(this)->postAggregateStats( stats ); } + +}; + +} /* namespace geos */ + +#endif /* SRC_CORECOMPONENTS_PHYSICSSOLVERS_STATISTICSAGGREGATOR_HPP_ */ diff --git a/src/coreComponents/physicsSolvers/StatisticsAggregatorBaseHelpers.hpp b/src/coreComponents/physicsSolvers/StatisticsAggregatorBaseHelpers.hpp new file mode 100644 index 00000000000..7f70116c8e4 --- /dev/null +++ b/src/coreComponents/physicsSolvers/StatisticsAggregatorBaseHelpers.hpp @@ -0,0 +1,291 @@ +/* + * ------------------------------------------------------------------------------------------------------------ + * 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. + * ------------------------------------------------------------------------------------------------------------ + */ + +#include "mesh/MeshBody.hpp" +#include "physicsSolvers/StatisticsAggregatorBase.hpp" + +#include "LvArray/src/math.hpp" +#include "common/DataTypes.hpp" +#include "common/format/Format.hpp" +#include "common/format/StringUtilities.hpp" +#include "common/logger/Logger.hpp" +#include "dataRepository/DataContext.hpp" +#include "dataRepository/Group.hpp" +#include "mesh/CellElementRegion.hpp" +#include "mesh/CellElementSubRegion.hpp" +#include "mesh/ElementRegionManager.hpp" +#include "mesh/MeshLevel.hpp" + +namespace geos +{ + +inline RegionStatisticsBase::RegionStatisticsBase( string const & targetName, + dataRepository::Group * const parent ): + dataRepository::Group( targetName, parent ), + m_time( std::numeric_limits< double >::lowest() ) +{} + +template< typename Impl > +StatsAggregatorBase< Impl >::StatsAggregatorBase( dataRepository::DataContext const & ownerDataContext ): + m_ownerDataContext( ownerDataContext ) +{} + +template< typename Impl > +void +StatsAggregatorBase< Impl >::initStatisticsAggregation( dataRepository::Group & meshBodies, + SolverType & solver ) +{ + m_meshBodies = &meshBodies; + + solver.forDiscretizationOnMeshTargets( meshBodies, [&] ( string const & meshBodyName, + MeshLevel & mesh, + string_array const & regionNames ) + { + // getting the container of all requesters statistics groups (can be already initialized) + dataRepository::Group * meshStatsGroup = mesh.getGroupPointer( ViewKeys::statisticsString() ); + if( meshStatsGroup == nullptr ) + meshStatsGroup = &mesh.registerGroup( ViewKeys::statisticsString() ); + + // registering the container of instance statistics groups (must be unique for this instance) + string const & ownerName = getOwnerName(); + GEOS_ERROR_IF_NE_MSG( meshStatsGroup->hasGroup( ownerName ), false, + GEOS_FMT( "A statistics aggregator have already been requested for '{}'.", + ownerName ), + m_ownerDataContext ); + meshStatsGroup->registerGroup( ownerName ); + + // remembering the path of this discretization + MeshBody const & body = m_meshBodies->getGroup< MeshBody >( meshBodyName ); + DiscretizationGroupPath const path { + body.getIndexInParent(), + mesh.getIndexInParent(), + regionNames, + }; + m_discretizationsPaths.push_back( path ); + } ); +} + +template< typename Impl > +void +StatsAggregatorBase< Impl >::enableRegionStatisticsAggregation( RegionStatsRegisterFunc && registerStatsFunc ) +{ + if( m_meshBodies == nullptr ) + return; + + integer regionCount = 0; + integer subRegionCount = 0; + + for( auto const & path : m_discretizationsPaths ) + { + MeshLevel & mesh = getMeshLevel( path ); + ElementRegionManager & elemManager = mesh.getElemManager(); + dataRepository::Group & statisticsGroup = getInstanceStatisticsGroup( mesh ); + StatsGroupType & meshRegionsStats = registerStatsFunc( statisticsGroup, + ViewKeys::regionsStatisticsString() ); + + for( size_t i = 0; i < path.regionNames.size(); ++i ) + { + CellElementRegion & region = elemManager.getRegion< CellElementRegion >( path.regionNames[i] ); + StatsGroupType & regionStats = registerStatsFunc( meshRegionsStats, + region.getName() ); + + region.forElementSubRegions< CellElementSubRegion >( [&] ( CellElementSubRegion & subRegion ) + { + registerStatsFunc( regionStats, + subRegion.getName() ); + ++subRegionCount; + } ); + ++regionCount; + } + } + + GEOS_ERROR_IF( regionCount == 0 || subRegionCount == 0, + GEOS_FMT( "Missing region for computing statistics:\n- {} regions,\n- {} sub-regions.", + getOwnerName(), regionCount, subRegionCount ), + m_ownerDataContext ); + + m_regionStatsState.m_isEnabled = true; + m_regionStatsState.m_isDirty = true; +} + +template< typename Impl > +void +StatsAggregatorBase< Impl >::forRegionStatistics( RegionStatsFunc< MeshLevel > const & func ) const +{ + for( auto const & path : m_discretizationsPaths ) + { + MeshLevel & mesh = getMeshLevel( path ); + StatsGroupType & meshRegionsStats = getRegionsStatistics( mesh ); + + func( mesh, meshRegionsStats ); + } +} + +template< typename Impl > +void +StatsAggregatorBase< Impl >::forRegionStatistics( MeshLevel & mesh, + StatsGroupType & meshRegionsStatistics, + RegionStatsFunc< CellElementRegion > const & func ) const +{ + ElementRegionManager & elemManager = mesh.getElemManager(); + meshRegionsStatistics.template forSubGroups< StatsGroupType >( [&] ( StatsGroupType & regionStatistics ) + { + string_view targetName = regionStatistics.getTargetName(); + CellElementRegion & region = elemManager.getRegion< CellElementRegion >( string( targetName ) ); + + func( region, regionStatistics ); + } ); +} + +template< typename Impl > +void +StatsAggregatorBase< Impl >::forRegionStatistics( CellElementRegion & region, + StatsGroupType & regionStatistics, + RegionStatsFunc< CellElementSubRegion > const & func ) const +{ + regionStatistics.template forSubGroups< StatsGroupType >( [&] ( StatsGroupType & subRegionStatistics ) + { + string_view targetName = subRegionStatistics.getTargetName(); + CellElementSubRegion & subRegion = region.getSubRegion< CellElementSubRegion >( string( targetName ) ); + func( subRegion, subRegionStatistics ); + } ); +} + +template< typename Impl > +bool +StatsAggregatorBase< Impl >::isComputed( real64 const timeRequest, StatsGroupType const & stats ) +{ + real64 const timePrecisionScale = LvArray::math::max( LvArray::math::abs( timeRequest ), + LvArray::math::abs( stats.m_time ) ); + static constexpr real64 timeRelTol = 1.0e-12; + + return + !m_regionStatsState.m_isDirty && + LvArray::math::abs( timeRequest - stats.m_time ) < timeRelTol * timePrecisionScale; +} + +template< typename Impl > +void +StatsAggregatorBase< Impl >::setDirty() +{ + m_regionStatsState.m_isDirty = true; +} + +template< typename Impl > +bool +StatsAggregatorBase< Impl >::computeRegionsStatistics( real64 const timeRequest ) +{ + GEOS_MARK_FUNCTION; + + m_warnings.clear(); + + // computation of sub region stats + forRegionStatistics( [&, timeRequest] ( MeshLevel & mesh, StatsGroupType & meshRegionsStats ) + { + forRegionStatistics( mesh, + meshRegionsStats, + [&, timeRequest] ( CellElementRegion & region, StatsGroupType & regionStats ) + { + forRegionStatistics( region, + regionStats, + [&, timeRequest] ( CellElementSubRegion & subRegion, StatsGroupType & subRegionStats ) + { + initStats( subRegionStats, timeRequest ); + computeSubRegionRankStats( subRegion, subRegionStats ); + } ); + } ); + } ); + + // aggregation of computations from the sub regions + forRegionStatistics( [&, timeRequest] ( MeshLevel & mesh, StatsGroupType & meshRegionsStats ) + { + initStats( meshRegionsStats, timeRequest ); + + forRegionStatistics( mesh, + meshRegionsStats, + [&, timeRequest] ( CellElementRegion & region, StatsGroupType & regionStats ) + { + initStats( regionStats, timeRequest ); + + forRegionStatistics( region, + regionStats, + [&, timeRequest] ( CellElementSubRegion &, StatsGroupType & subRegionStats ) + { + aggregateStats( regionStats, subRegionStats ); + + mpiAggregateStats( subRegionStats ); + postAggregateStats( subRegionStats ); + } ); + + aggregateStats( meshRegionsStats, regionStats ); + + mpiAggregateStats( regionStats ); + postAggregateStats( regionStats ); + } ); + + mpiAggregateStats( meshRegionsStats ); + postAggregateStats( meshRegionsStats ); + } ); + + m_regionStatsState.m_isDirty = false; + + return true; +} + +template< typename Impl > +MeshLevel & +StatsAggregatorBase< Impl >::getMeshLevel( DiscretizationGroupPath const & path ) const +{ + MeshBody & body = m_meshBodies->getGroup< MeshBody >( path.meshBody ); + MeshLevel & mesh = body.getMeshLevel( path.meshLevel ); + return mesh; +} + +template< typename Impl > +dataRepository::Group & +StatsAggregatorBase< Impl >::getInstanceStatisticsGroup( MeshLevel & mesh ) const +{ + // considering everything is initialized, or else, crash gracefully + dataRepository::Group & meshStatsGroup = mesh.getGroup( ViewKeys::statisticsString() ); + dataRepository::Group & instanceStatisticsGroup = meshStatsGroup.getGroup( getOwnerName() ); + return instanceStatisticsGroup; +} + +template< typename Impl > +typename StatsAggregatorBase< Impl >::StatsGroupType & +StatsAggregatorBase< Impl >:: +getRegionsStatistics( MeshLevel & mesh ) const +{ + // considering everything is initialized, or else, crash gracefully + dataRepository::Group & instanceStatisticsGroup = getInstanceStatisticsGroup( mesh ); + return instanceStatisticsGroup.getGroup< StatsGroupType >( ViewKeys::regionsStatisticsString() ); +} + +template< typename Impl > +typename StatsAggregatorBase< Impl >::StatsGroupType & +StatsAggregatorBase< Impl >::getRegionStatistics( MeshLevel & mesh, + string_view regionName ) const +{ + StatsGroupType & meshRegionsStats = getRegionsStatistics( mesh ); + StatsGroupType * const stats = meshRegionsStats.template getGroupPointer< StatsGroupType >( string( regionName ) ); + GEOS_THROW_IF( stats == nullptr, + GEOS_FMT( "Region '{}' not found to get region statistics, is it a target of the reservoir solver?\n" + "Available target regions:\n- {}", + regionName, stringutilities::join( meshRegionsStats.getSubGroupsNames(), "\n- " ) ), + InputError, m_ownerDataContext ); + return *stats; +} + +} /* namespace geos */ diff --git a/src/coreComponents/physicsSolvers/fluidFlow/CMakeLists.txt b/src/coreComponents/physicsSolvers/fluidFlow/CMakeLists.txt index 6eb55e9775a..da6001bf435 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 @@ -34,7 +35,8 @@ set( fluidFlowSolvers_headers SourceFluxStatistics.hpp SinglePhaseBase.hpp SinglePhaseBaseFields.hpp - SinglePhaseStatistics.hpp + SinglePhaseStatisticsAggregator.hpp + SinglePhaseStatisticsTask.hpp SinglePhaseFVM.hpp SinglePhaseHybridFVM.hpp SinglePhaseProppantBase.hpp @@ -142,13 +144,15 @@ set( fluidFlowSolvers_headers set( fluidFlowSolvers_sources CompositionalMultiphaseBase.cpp CompositionalMultiphaseFVM.cpp - CompositionalMultiphaseStatistics.cpp + CompositionalMultiphaseStatisticsAggregator.cpp + CompositionalMultiphaseStatisticsTask.cpp CompositionalMultiphaseHybridFVM.cpp ImmiscibleMultiphaseFlow.cpp ReactiveCompositionalMultiphaseOBL.cpp FlowSolverBase.cpp SinglePhaseBase.cpp - SinglePhaseStatistics.cpp + SinglePhaseStatisticsAggregator.cpp + SinglePhaseStatisticsTask.cpp SinglePhaseFVM.cpp SinglePhaseHybridFVM.cpp SinglePhaseProppantBase.cpp diff --git a/src/coreComponents/physicsSolvers/fluidFlow/CompositionalMultiphaseStatistics.cpp b/src/coreComponents/physicsSolvers/fluidFlow/CompositionalMultiphaseStatistics.cpp deleted file mode 100644 index c83695e5a98..00000000000 --- a/src/coreComponents/physicsSolvers/fluidFlow/CompositionalMultiphaseStatistics.cpp +++ /dev/null @@ -1,555 +0,0 @@ -/* - * ------------------------------------------------------------------------------------------------------------ - * 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 "CompositionalMultiphaseStatistics.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 "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" - - -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 ) -{ - 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(); - - if( dynamicCast< CompositionalMultiphaseHybridFVM * >( m_solver ) && m_computeCFLNumbers != 0 ) - { - GEOS_THROW( "The option to compute CFL numbers is incompatible with CompositionalMultiphaseHybridFVM", - InputError, getDataContext() ); - } -} - -void CompositionalMultiphaseStatistics::registerDataOnMesh( 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; - } - - m_solver->forDiscretizationOnMeshTargets( meshBodies, [&] ( string const &, - MeshLevel & mesh, - string_array const & regionNames ) - { - ElementRegionManager & elemManager = mesh.getElemManager(); - - 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 ) - { - 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( GEOS_FMT( "Time [{}]", units::getSymbol( units::Unit::Time ))), - 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 ) )), - } ); - - 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(); - } - } - } - } ); - - // if we have to compute CFL numbers later, we need to register additional variables - if( m_computeCFLNumbers ) - { - m_solver->registerDataForCFL( meshBodies ); - } -} - -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 ) -{ - 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 ); - } - } ); - - if( m_computeCFLNumbers ) - { - // current time is time_n + dt - computeCFLNumbers( time_n + dt, dt, domain ); - } - - return false; -} - -void CompositionalMultiphaseStatistics::computeRegionStatistics( real64 const time, - MeshLevel & mesh, - string_array const & regionNames ) const -{ - GEOS_MARK_FUNCTION; - - integer const numPhases = m_solver->numFluidPhases(); - integer const numComps = m_solver->numFluidComponents(); - - // Step 1: initialize the average/min/max quantities - ElementRegionManager & elemManager = mesh.getElemManager(); - for( size_t i = 0; i < regionNames.size(); ++i ) - { - ElementRegionBase & region = elemManager.getRegion( regionNames[i] ); - RegionStatistics & stats = region.getReference< RegionStatistics >( viewKeyStruct::regionStatisticsString() ); - - stats.averagePressure = 0.0; - stats.maxPressure = 0.0; - stats.minPressure = LvArray::NumericLimits< real64 >::max; - - stats.maxDeltaPressure = -LvArray::NumericLimits< real64 >::max; - stats.minDeltaPressure = LvArray::NumericLimits< real64 >::max; - - stats.averageTemperature = 0.0; - stats.maxTemperature = 0.0; - stats.minTemperature = LvArray::NumericLimits< real64 >::max; - - stats.totalPoreVolume = 0.0; - stats.totalUncompactedPoreVolume = 0.0; - stats.phasePoreVolume.setValues< serialPolicy >( 0.0 ); - - stats.phaseMass.setValues< serialPolicy >( 0.0 ); - stats.trappedPhaseMass.setValues< serialPolicy >( 0.0 ); - stats.immobilePhaseMass.setValues< serialPolicy >( 0.0 ); - stats.componentMass.setValues< serialPolicy >( 0.0 ); - } - - // Step 2: increment the average/min/max quantities for all the subRegions - elemManager.forElementSubRegions( regionNames, [&]( localIndex const, - ElementSubRegionBase & subRegion ) - { - - 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 >(); - arrayView2d< real64 const, compflow::USD_PHASE > const phaseVolFrac = - subRegion.getField< flow::phaseVolumeFraction >(); - arrayView1d< real64 const > const deltaPres = subRegion.getField< flow::deltaPressure >(); - - Group const & constitutiveModels = subRegion.getGroup( ElementSubRegionBase::groupKeyStruct::constitutiveModelsString() ); - - string const & solidName = subRegion.getReference< string >( CompositionalMultiphaseBase::viewKeyStruct::solidNamesString() ); - CoupledSolidBase const & solid = constitutiveModels.getGroup< CoupledSolidBase >( solidName ); - arrayView1d< real64 const > const refPorosity = solid.getReferencePorosity(); - arrayView2d< real64 const > const porosity = solid.getPorosity(); - - string const & fluidName = subRegion.getReference< string >( CompositionalMultiphaseBase::viewKeyStruct::fluidNamesString() ); - MultiFluidBase const & fluid = constitutiveModels.getGroup< MultiFluidBase >( fluidName ); - arrayView3d< real64 const, multifluid::USD_PHASE > const phaseDensity = fluid.phaseDensity(); - arrayView4d< real64 const, multifluid::USD_PHASE_COMP > const phaseCompFraction = fluid.phaseCompFraction(); - - - //get min vol fraction for each phase to dispactche immobile/mobile mass - string const & relpermName = subRegion.getReference< string >( CompositionalMultiphaseBase::viewKeyStruct::relPermNamesString() ); - RelativePermeabilityBase const & relperm = constitutiveModels.getGroup< RelativePermeabilityBase >( relpermName ); - arrayView3d< real64 const, relperm::USD_RELPERM > const phaseTrappedVolFrac = relperm.phaseTrappedVolFraction(); - arrayView3d< real64 const, relperm::USD_RELPERM > const phaseRelperm = relperm.phaseRelPerm(); - - real64 subRegionAvgPresNumerator = 0.0; - real64 subRegionMinPres = 0.0; - real64 subRegionMaxPres = 0.0; - real64 subRegionMinDeltaPres = 0.0; - real64 subRegionMaxDeltaPres = 0.0; - real64 subRegionAvgTempNumerator = 0.0; - real64 subRegionMinTemp = 0.0; - real64 subRegionMaxTemp = 0.0; - real64 subRegionTotalUncompactedPoreVol = 0.0; - array1d< real64 > subRegionPhaseDynamicPoreVol( numPhases ); - array1d< real64 > subRegionPhaseMass( numPhases ); - array1d< real64 > subRegionTrappedPhaseMass( numPhases ); - array1d< real64 > subRegionImmobilePhaseMass( numPhases ); - array1d< real64 > subRegionRelpermPhaseMass( numPhases ); - array2d< real64 > subRegionComponentMass( numPhases, numComps ); - - isothermalCompositionalMultiphaseBaseKernels:: - StatisticsKernel:: - launch< parallelDevicePolicy<> >( subRegion.size(), - numComps, - numPhases, - m_relpermThreshold, - elemGhostRank, - volume, - pres, - deltaPres, - temp, - refPorosity, - porosity, - phaseDensity, - phaseCompFraction, - phaseVolFrac, - phaseTrappedVolFrac, - phaseRelperm, - subRegionMinPres, - subRegionAvgPresNumerator, - subRegionMaxPres, - subRegionMinDeltaPres, - subRegionMaxDeltaPres, - subRegionMinTemp, - subRegionAvgTempNumerator, - subRegionMaxTemp, - subRegionTotalUncompactedPoreVol, - subRegionPhaseDynamicPoreVol.toView(), - subRegionPhaseMass.toView(), - subRegionTrappedPhaseMass.toView(), - subRegionImmobilePhaseMass.toView(), - subRegionComponentMass.toView() ); - - ElementRegionBase & region = elemManager.getRegion( ElementRegionBase::getParentRegion( subRegion ).getName() ); - RegionStatistics & stats = region.getReference< RegionStatistics >( viewKeyStruct::regionStatisticsString() ); - - stats.averagePressure += subRegionAvgPresNumerator; - if( subRegionMinPres < stats.minPressure ) - { - stats.minPressure = subRegionMinPres; - } - if( subRegionMaxPres > stats.maxPressure ) - { - stats.maxPressure = subRegionMaxPres; - } - - if( subRegionMinDeltaPres < stats.minDeltaPressure ) - { - stats.minDeltaPressure = subRegionMinDeltaPres; - } - if( subRegionMaxDeltaPres > stats.maxDeltaPressure ) - { - stats.maxDeltaPressure = subRegionMaxDeltaPres; - } - - stats.averageTemperature += subRegionAvgTempNumerator; - if( subRegionMinTemp < stats.minTemperature ) - { - stats.minTemperature = subRegionMinTemp; - } - if( subRegionMaxTemp > stats.maxTemperature ) - { - stats.maxTemperature = subRegionMaxTemp; - } - - stats.totalUncompactedPoreVolume += subRegionTotalUncompactedPoreVol; - for( integer ip = 0; ip < numPhases; ++ip ) - { - stats.phasePoreVolume[ip] += subRegionPhaseDynamicPoreVol[ip]; - stats.phaseMass[ip] += subRegionPhaseMass[ip]; - stats.trappedPhaseMass[ip] += subRegionTrappedPhaseMass[ip]; - stats.immobilePhaseMass[ip] += subRegionImmobilePhaseMass[ip]; - - for( integer ic = 0; ic < numComps; ++ic ) - { - stats.componentMass[ip][ic] += subRegionComponentMass[ip][ic]; - } - } - - } ); - - // Step 3: synchronize the results over the MPI ranks - for( size_t i = 0; i < regionNames.size(); ++i ) - { - ElementRegionBase & region = elemManager.getRegion( regionNames[i] ); - RegionStatistics & stats = region.getReference< RegionStatistics >( viewKeyStruct::regionStatisticsString() ); - - stats.minPressure = MpiWrapper::min( stats.minPressure ); - stats.maxPressure = MpiWrapper::max( stats.maxPressure ); - stats.minDeltaPressure = MpiWrapper::min( stats.minDeltaPressure ); - stats.maxDeltaPressure = MpiWrapper::max( stats.maxDeltaPressure ); - stats.minTemperature = MpiWrapper::min( stats.minTemperature ); - stats.maxTemperature = MpiWrapper::max( stats.maxTemperature ); - stats.totalUncompactedPoreVolume = MpiWrapper::sum( stats.totalUncompactedPoreVolume ); - stats.totalPoreVolume = 0.0; - for( integer ip = 0; ip < numPhases; ++ip ) - { - stats.phasePoreVolume[ip] = MpiWrapper::sum( stats.phasePoreVolume[ip] ); - stats.phaseMass[ip] = MpiWrapper::sum( stats.phaseMass[ip] ); - stats.trappedPhaseMass[ip] = MpiWrapper::sum( stats.trappedPhaseMass[ip] ); - stats.immobilePhaseMass[ip] = MpiWrapper::sum( stats.immobilePhaseMass[ip] ); - stats.totalPoreVolume += stats.phasePoreVolume[ip]; - for( integer ic = 0; ic < numComps; ++ic ) - { - stats.componentMass[ip][ic] = MpiWrapper::sum( stats.componentMass[ip][ic] ); - } - } - stats.averagePressure = MpiWrapper::sum( stats.averagePressure ); - stats.averageTemperature = MpiWrapper::sum( stats.averageTemperature ); - if( stats.totalUncompactedPoreVolume > 0 ) - { - float invTotalUncompactedPoreVolume = 1.0 / stats.totalUncompactedPoreVolume; - stats.averagePressure *= invTotalUncompactedPoreVolume; - stats.averageTemperature *= invTotalUncompactedPoreVolume; - } - else - { - 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] ) ); - } - - - // 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]; - } - - string_view massUnit = units::getSymbol( m_solver->getMassUnit() ); - - stdVector< string > phaseCompName; - phaseCompName.reserve( numPhases*numComps ); - stdVector< string > massValues; - phaseCompName.reserve( numPhases*numComps ); - - ConstitutiveManager const & constitutiveManager = this->getGroupByPath< ConstitutiveManager >( "/Problem/domain/Constitutive" ); - MultiFluidBase const & fluid = constitutiveManager.getGroup< MultiFluidBase >( m_solver->referenceFluidModelName() ); - auto const phaseNames = fluid.phaseNames(); - auto const componentNames = fluid.componentNames(); - for( integer ip = 0; ip < numPhases; ++ip ) - { - for( integer ic = 0; ic < numComps; ++ic ) - { - std::stringstream ss; - ss << phaseNames[ip] << "/" << componentNames[ic]; - phaseCompName.push_back( ss.str() ); - massValues.push_back( GEOS_FMT( "{}", stats.componentMass[ip][ic] ) ); - } - } - - if( isLogLevelActive< logInfo::Statistics >( this->getLogLevel() ) && MpiWrapper::commRank() == 0 ) - { - TableData compPhaseStatsData; - compPhaseStatsData.addRow( "Pressure [Pa]", stats.minPressure, stats.averagePressure, stats.maxPressure ); - compPhaseStatsData.addRow( "Delta pressure [Pa]", stats.minDeltaPressure, "/", stats.maxDeltaPressure ); - compPhaseStatsData.addRow( "Temperature [K]", stats.minTemperature, stats.averageTemperature, stats.maxTemperature ); - compPhaseStatsData.addSeparator(); - - compPhaseStatsData.addRow( "Total dynamic pore volume [rm^3]", CellType::MergeNext, CellType::MergeNext, stats.totalPoreVolume ); - compPhaseStatsData.addSeparator(); - compPhaseStatsData.addRow( "Phase dynamic pore volume [rm^3]", - stringutilities::joinLambda( phaseNames, "\n", []( auto data ) { return data[0]; } ), - CellType::MergeNext, - stringutilities::joinLambda( stats.phasePoreVolume, "\n", []( auto data ) { return data[0]; } ) ); - compPhaseStatsData.addSeparator(); - - compPhaseStatsData.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]; } ) ); - compPhaseStatsData.addSeparator(); - - compPhaseStatsData.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]; } ) ); - compPhaseStatsData.addSeparator(); - compPhaseStatsData.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]; } ) ); - compPhaseStatsData.addSeparator(); - - compPhaseStatsData.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]; } ) ); - compPhaseStatsData.addSeparator(); - compPhaseStatsData.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]; } ) ); - compPhaseStatsData.addSeparator(); - - compPhaseStatsData.addRow( GEOS_FMT( "Component mass [{}]", massUnit ), - stringutilities::join( phaseCompName, '\n' ), - CellType::MergeNext, - stringutilities::join( massValues, '\n' ) ); - - string const title = GEOS_FMT( "{}, {} (time {} s):", getName(), regionNames[i], time ); - TableLayout const compPhaseStatsLayout( title, { "statistics", "min", "average", "max" } ); - TableTextFormatter tableFormatter( compPhaseStatsLayout ); - GEOS_LOG_RANK_0( tableFormatter.toString( compPhaseStatsData ) ); - } - - if( m_writeCSV > 0 && MpiWrapper::commRank() == 0 ) - { - TableData tableData; - tableData.addRow( time, stats.minPressure, stats.averagePressure, stats.maxPressure, stats.minDeltaPressure, stats.maxDeltaPressure, - stats.minTemperature, stats.averageTemperature, stats.maxTemperature, stats.totalPoreVolume, - stringutilities::joinLambda( stats.phasePoreVolume, "\n", []( auto data ) { return data[0]; } ), - stringutilities::joinLambda( stats.phaseMass, "\n", []( auto data ) { return data[0]; } ), - stringutilities::joinLambda( stats.trappedPhaseMass, "\n", []( auto value ) { return value[0]; } ), - stringutilities::joinLambda( nonTrappedPhaseMass, "\n", []( auto value ) { return value[0]; } ), - stringutilities::joinLambda( stats.immobilePhaseMass, "\n", []( auto value ) { return value[0]; } ), - stringutilities::joinLambda( mobilePhaseMass, "\n", []( auto value ) { return value[0]; } ), - stringutilities::join( massValues, '\n' ) ); - - std::ofstream outputFile( m_outputDir + "/" + regionNames[i] + ".csv", std::ios_base::app ); - TableCSVFormatter const csvOutput; - outputFile << csvOutput.dataToString( tableData ); - outputFile.close(); - } - } - -} - -void CompositionalMultiphaseStatistics::computeCFLNumbers( real64 const time, - real64 const dt, - DomainPartition & domain ) const -{ - GEOS_MARK_FUNCTION; - real64 maxPhaseCFL, maxCompCFL; - m_solver->computeCFLNumbers( domain, dt, maxPhaseCFL, maxCompCFL ); - - GEOS_LOG_LEVEL_RANK_0( logInfo::CFL, - GEOS_FMT( "{} (time {} s): Max phase CFL number: {}", getName(), time, maxPhaseCFL ) ); - GEOS_LOG_LEVEL_RANK_0( logInfo::CFL, - GEOS_FMT( "{} (time {} s): Max component CFL number: {}", getName(), time, maxCompCFL ) ); -} - - -REGISTER_CATALOG_ENTRY( TaskBase, - CompositionalMultiphaseStatistics, - string const &, dataRepository::Group * const ) - -} /* namespace geos */ diff --git a/src/coreComponents/physicsSolvers/fluidFlow/CompositionalMultiphaseStatisticsAggregator.cpp b/src/coreComponents/physicsSolvers/fluidFlow/CompositionalMultiphaseStatisticsAggregator.cpp new file mode 100644 index 00000000000..734b6389e59 --- /dev/null +++ b/src/coreComponents/physicsSolvers/fluidFlow/CompositionalMultiphaseStatisticsAggregator.cpp @@ -0,0 +1,383 @@ +/* + * ------------------------------------------------------------------------------------------------------------ + * 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 + * @details Region statistics data is stored as follow: + + * Problem : ProblemManager + * |-> domain : DomainPartition + * |-> MeshBodies : Group + * |-> cartesianMesh : MeshBody + * |-> meshLevels : Group + * |-> Level0 : MeshLevel + * | |-> nodeManager : NodeManager + * | | |-> sets : Group + * | | | * all : Wrapper< index array > + * | | | * xneg : Wrapper< index array > + * | | [...] (other element sets) + * | | + * | |-> ElementRegions : ElementRegionManager + * | | |-> Channel : CellElementRegion + * | | | |-> cb-0_0_0 : CellElementSubRegion + * | | | | | * pressure : Wrapper< real64 array > + * | | | | | * temperature : Wrapper< real64 array > + * | | | | [...] (other fields) + * | | | | + * | | | |-> cb-0_0_1 : CellElementSubRegion + * | | | | | * pressure : Wrapper< real64 array > + * | | | | | * temperature : Wrapper< real64 array > + * | | | | [...] (other fields) + * | | | | + * | | | [...] (other sub-regions) + * | | | + * | | |-> Barrier : CellElementRegion + * | | |-> cb-1_0_0 : CellElementSubRegion + * | | |-> cb-1_0_1 : CellElementSubRegion + * | | [...] (other sub-regions) + * | | + * | [...] (other element managers) + * ____ | | + * | | |-> statistics : Group (storage for all stats) + * | | |-> compFlowStats : Group (storage for this instance stats) + * | | | |-> cflStatistics : CFLStatistics + * | | | |-> regionsStatistics : RegionStatistics (aggregate) + * | | | |-> Channel : RegionStatistics (aggregate, mpi reduced) + * | | | | |-> cb-0_0_0 : RegionStatistics (compute read-back) + * stats | | | | |-> cb-0_0_1 : RegionStatistics (compute read-back) + * data -> | | | | [...] (other sub-regions stats) + * | | | | + * | | | |-> Barrier : RegionStatistics (aggregate, mpi reduced) + * | | | |-> cb-1_0_0 : RegionStatistics (compute read-back) + * | | | |-> cb-1_0_1 : RegionStatistics (compute read-back) + * | | | [...] (other sub-regions stats) + * | | | + * |___ | [...] (other stats storages) + * | + * [...] (other discretizations) + */ + +#include "CompositionalMultiphaseStatisticsAggregator.hpp" + +#include "physicsSolvers/StatisticsAggregatorBaseHelpers.hpp" +#include "mesh/DomainPartition.hpp" +#include "physicsSolvers/fluidFlow/CompositionalMultiphaseBase.hpp" +#include "physicsSolvers/fluidFlow/kernels/compositional/StatisticsKernel.hpp" +#include "constitutive/relativePermeability/RelativePermeabilityBase.hpp" + +namespace geos +{ + +using namespace constitutive; +using namespace dataRepository; + +namespace compositionalMultiphaseStatistics +{ + +RegionStatistics::RegionStatistics( string const & name, dataRepository::Group * const parent, + integer const numPhases, integer const numComponents ): + RegionStatisticsBase( name, parent ), + m_phaseDynamicPoreVolume( numPhases ), + m_phaseMass( numPhases ), + m_trappedPhaseMass( numPhases ), + m_nonTrappedPhaseMass( numPhases ), + m_immobilePhaseMass( numPhases ), + m_mobilePhaseMass( numPhases ), + m_componentMass( numPhases, numComponents ) +{ + // TODO : registerWrappers to store results in HDF5 (but need repairing of 1D HDF5 outputs) +} + +CFLStatistics::CFLStatistics( string const & name, dataRepository::Group * const parent ): + RegionStatisticsBase( name, parent ) +{ + // TODO : registerWrappers to store results in HDF5 (but need repairing of 1D HDF5 outputs) +} + +StatsAggregator::StatsAggregator( DataContext const & ownerDataContext ): + Base( ownerDataContext ), + m_params() +{} + +void StatsAggregator::initStatisticsAggregation( dataRepository::Group & meshBodies, + CompositionalMultiphaseBase & solver ) +{ + m_solver = &solver; + m_numPhases = solver.numFluidPhases(); + m_numComponents = solver.numFluidComponents(); + + Base::initStatisticsAggregation( meshBodies, solver ); +} + +void StatsAggregator::enableRegionStatisticsAggregation() +{ + auto const registerStats = [=] ( Group & parent, + string const & targetName ) -> RegionStatistics & + { + return parent.registerGroup( targetName, + std::make_unique< RegionStatistics >( targetName, &parent, + m_numPhases, + m_numComponents ) ); + }; + + Base::enableRegionStatisticsAggregation( registerStats ); +} + +void StatsAggregator::enableCFLStatistics() +{ + if( m_solver == nullptr ) + return; + + m_solver->registerDataForCFL( *m_meshBodies ); + for( auto const & path : m_discretizationsPaths ) + { + MeshLevel & mesh = getMeshLevel( path ); + Group & statisticsGroup = getInstanceStatisticsGroup( mesh ); + statisticsGroup.registerGroup< CFLStatistics >( ViewKeys::cflStatisticsString() ); + } + + m_cflStatsState.m_isEnabled = true; + m_cflStatsState.m_isDirty = true; +} + +void StatsAggregator::setDirty() +{ + Base::setDirty(); + m_cflStatsState.m_isDirty = true; +} + +bool StatsAggregator::computeCFLNumbers( real64 const time, + real64 const dt, + DomainPartition & domain ) +{ + GEOS_MARK_FUNCTION; + real64 maxPhaseCFL, maxCompCFL; + CFLStatistics * stats = getCFLStatistics( domain ); + + m_warnings.clear(); + + if( !m_cflStatsState.m_isEnabled ) + { + m_warnings.emplace_back( "CFL numbers computation is not enabled." ); + return false; + } + if( stats == nullptr ) + { + m_warnings.emplace_back( GEOS_FMT( "No statistics structure to compute CFL numbers for domain '{}'.", domain.getName() )); + return false; + } + + m_solver->computeCFLNumbers( domain, dt, maxPhaseCFL, maxCompCFL ); + + stats->m_time = time; + stats->m_maxPhaseCFL = maxPhaseCFL; + stats->m_maxCompCFL = maxCompCFL; + + m_cflStatsState.m_isDirty = false; + + return true; +} + +CFLStatistics * StatsAggregator::getCFLStatistics( DomainPartition & domain ) const +{ + return domain.getGroupPointer< CFLStatistics >( ViewKeys::cflStatisticsString() ); +} + +CFLStatistics & StatsAggregator::getCflStatistics( MeshLevel & mesh ) const +{ + // considering everything is initialized, or else, crash gracefully + Group & statisticsGroup = getInstanceStatisticsGroup( mesh ); + return statisticsGroup.getGroup< CFLStatistics >( ViewKeys::cflStatisticsString() ); +} + +void StatsAggregator::initStats( RegionStatistics & stats, real64 const time ) const +{ + stats.m_time = time; + + stats.m_averagePressure = 0.0; + stats.m_maxPressure = 0.0; + stats.m_minPressure = LvArray::NumericLimits< real64 >::max; + + stats.m_maxDeltaPressure = -LvArray::NumericLimits< real64 >::max; + stats.m_minDeltaPressure = LvArray::NumericLimits< real64 >::max; + + stats.m_averageTemperature = 0.0; + stats.m_maxTemperature = 0.0; + stats.m_minTemperature = LvArray::NumericLimits< real64 >::max; + + stats.m_totalPoreVolume = 0.0; + stats.m_totalUncompactedPoreVolume = 0.0; + stats.m_phaseDynamicPoreVolume.setValues< serialPolicy >( 0.0 ); + + stats.m_phaseMass.setValues< serialPolicy >( 0.0 ); + stats.m_trappedPhaseMass.setValues< serialPolicy >( 0.0 ); + stats.m_nonTrappedPhaseMass.setValues< serialPolicy >( 0.0 ); + stats.m_immobilePhaseMass.setValues< serialPolicy >( 0.0 ); + stats.m_mobilePhaseMass.setValues< serialPolicy >( 0.0 ); + stats.m_componentMass.setValues< serialPolicy >( 0.0 ); +} + +void StatsAggregator::computeSubRegionRankStats( CellElementSubRegion & subRegion, + RegionStatistics & subRegionStats ) const +{ + arrayView1d< integer const > const elemGhostRank = subRegion.ghostRank(); + arrayView1d< real64 const > const volume = subRegion.getElementVolume(); + 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< fields::flow::phaseVolumeFraction >(); + arrayView1d< real64 const > const deltaPres = subRegion.getField< fields::flow::deltaPressure >(); + + Group const & constitutiveModels = subRegion.getGroup( ElementSubRegionBase::groupKeyStruct::constitutiveModelsString() ); + + string const & solidName = subRegion.getReference< string >( CompositionalMultiphaseBase::viewKeyStruct::solidNamesString() ); + CoupledSolidBase const & solid = constitutiveModels.getGroup< CoupledSolidBase >( solidName ); + arrayView1d< real64 const > const refPorosity = solid.getReferencePorosity(); + arrayView2d< real64 const > const porosity = solid.getPorosity(); + + string const & fluidName = subRegion.getReference< string >( CompositionalMultiphaseBase::viewKeyStruct::fluidNamesString() ); + MultiFluidBase const & fluid = constitutiveModels.getGroup< MultiFluidBase >( fluidName ); + arrayView3d< real64 const, multifluid::USD_PHASE > const phaseDensity = fluid.phaseDensity(); + arrayView4d< real64 const, multifluid::USD_PHASE_COMP > const phaseCompFraction = fluid.phaseCompFraction(); + + //get min vol fraction for each phase to dispactche immobile/mobile mass + string const & relpermName = subRegion.getReference< string >( CompositionalMultiphaseBase::viewKeyStruct::relPermNamesString() ); + RelativePermeabilityBase const & relperm = constitutiveModels.getGroup< RelativePermeabilityBase >( relpermName ); + arrayView3d< real64 const, relperm::USD_RELPERM > const phaseTrappedVolFrac = relperm.phaseTrappedVolFraction(); + arrayView3d< real64 const, relperm::USD_RELPERM > const phaseRelperm = relperm.phaseRelPerm(); + + isothermalCompositionalMultiphaseBaseKernels:: + StatisticsKernel:: + launch< parallelDevicePolicy<> >( subRegion.size(), + m_numComponents, + m_numPhases, + m_params.m_relpermThreshold, + elemGhostRank, + volume, + pres, + deltaPres, + temp, + refPorosity, + porosity, + phaseDensity, + phaseCompFraction, + phaseVolFrac, + phaseTrappedVolFrac, + phaseRelperm, + subRegionStats.m_minPressure, + subRegionStats.m_averagePressure, + subRegionStats.m_maxPressure, + subRegionStats.m_minDeltaPressure, + subRegionStats.m_maxDeltaPressure, + subRegionStats.m_minTemperature, + subRegionStats.m_averageTemperature, + subRegionStats.m_maxTemperature, + subRegionStats.m_totalUncompactedPoreVolume, + subRegionStats.m_phaseDynamicPoreVolume.toView(), + subRegionStats.m_phaseMass.toView(), + subRegionStats.m_trappedPhaseMass.toView(), + subRegionStats.m_immobilePhaseMass.toView(), + subRegionStats.m_componentMass.toView() ); +} + + +void StatsAggregator::aggregateStats( RegionStatistics & stats, + RegionStatistics const & other ) const +{ + stats.m_averagePressure += other.m_averagePressure; + stats.m_minPressure = LvArray::math::min( stats.m_minPressure, other.m_minPressure ); + stats.m_maxPressure = LvArray::math::max( stats.m_maxPressure, other.m_maxPressure ); + + stats.m_minDeltaPressure = LvArray::math::min( stats.m_minDeltaPressure, other.m_minDeltaPressure ); + stats.m_maxDeltaPressure = LvArray::math::max( stats.m_maxDeltaPressure, other.m_maxDeltaPressure ); + + stats.m_averageTemperature += other.m_averageTemperature; + stats.m_minTemperature = LvArray::math::min( stats.m_minTemperature, other.m_minTemperature ); + stats.m_maxTemperature = LvArray::math::max( stats.m_maxTemperature, other.m_maxTemperature ); + + stats.m_totalUncompactedPoreVolume += other.m_totalUncompactedPoreVolume; + + for( integer ip = 0; ip < m_numPhases; ++ip ) + { + stats.m_phaseDynamicPoreVolume[ip] += other.m_phaseDynamicPoreVolume[ip]; + stats.m_phaseMass[ip] += other.m_phaseMass[ip]; + stats.m_trappedPhaseMass[ip] += other.m_trappedPhaseMass[ip]; + stats.m_immobilePhaseMass[ip] += other.m_immobilePhaseMass[ip]; + + for( integer ic = 0; ic < m_numComponents; ++ic ) + { + stats.m_componentMass[ip][ic] += other.m_componentMass[ip][ic]; + } + } +} + +void StatsAggregator::mpiAggregateStats( RegionStatistics & stats ) const +{ + stats.m_averagePressure = MpiWrapper::sum( stats.m_averagePressure ); + stats.m_minPressure = MpiWrapper::min( stats.m_minPressure ); + stats.m_maxPressure = MpiWrapper::max( stats.m_maxPressure ); + + stats.m_minDeltaPressure = MpiWrapper::min( stats.m_minDeltaPressure ); + stats.m_maxDeltaPressure = MpiWrapper::max( stats.m_maxDeltaPressure ); + + stats.m_averageTemperature = MpiWrapper::sum( stats.m_averageTemperature ); + stats.m_minTemperature = MpiWrapper::min( stats.m_minTemperature ); + stats.m_maxTemperature = MpiWrapper::max( stats.m_maxTemperature ); + + stats.m_totalUncompactedPoreVolume = MpiWrapper::sum( stats.m_totalUncompactedPoreVolume ); + + for( integer ip = 0; ip < m_numPhases; ++ip ) + { + stats.m_phaseDynamicPoreVolume[ip] = MpiWrapper::sum( stats.m_phaseDynamicPoreVolume[ip] ); + stats.m_phaseMass[ip] = MpiWrapper::sum( stats.m_phaseMass[ip] ); + stats.m_trappedPhaseMass[ip] = MpiWrapper::sum( stats.m_trappedPhaseMass[ip] ); + stats.m_immobilePhaseMass[ip] = MpiWrapper::sum( stats.m_immobilePhaseMass[ip] ); + + for( integer ic = 0; ic < m_numComponents; ++ic ) + { + stats.m_componentMass[ip][ic] = MpiWrapper::sum( stats.m_componentMass[ip][ic] ); + } + } +} + +void StatsAggregator::postAggregateStats( RegionStatistics & stats ) +{ + if( stats.m_totalUncompactedPoreVolume > 0 ) + { + float invTotalUncompactedPoreVolume = 1.0 / stats.m_totalUncompactedPoreVolume; + stats.m_averagePressure *= invTotalUncompactedPoreVolume; + stats.m_averageTemperature *= invTotalUncompactedPoreVolume; + } + else + { + stats.m_averagePressure = 0.0; + stats.m_averageTemperature = 0.0; + m_warnings.emplace_back( GEOS_FMT( "Cannot compute average pressure for '{}' because pore volume is zero in '{}'.", + getOwnerName(), stats.getTargetName() ) ); + } + + for( integer ip = 0; ip < m_numPhases; ++ip ) + { + stats.m_totalPoreVolume += stats.m_phaseDynamicPoreVolume[ip]; + stats.m_nonTrappedPhaseMass[ip] = stats.m_phaseMass[ip] - stats.m_trappedPhaseMass[ip]; + stats.m_mobilePhaseMass[ip] = stats.m_phaseMass[ip] - stats.m_immobilePhaseMass[ip]; + } +} + +} /* namespace compositionalMultiphaseStatistics */ + +template class StatsAggregatorBase< compositionalMultiphaseStatistics::StatsAggregator >; + +} /* namespace geos */ diff --git a/src/coreComponents/physicsSolvers/fluidFlow/CompositionalMultiphaseStatisticsAggregator.hpp b/src/coreComponents/physicsSolvers/fluidFlow/CompositionalMultiphaseStatisticsAggregator.hpp new file mode 100644 index 00000000000..d5857d028f7 --- /dev/null +++ b/src/coreComponents/physicsSolvers/fluidFlow/CompositionalMultiphaseStatisticsAggregator.hpp @@ -0,0 +1,305 @@ +/* + * ------------------------------------------------------------------------------------------------------------ + * 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 CompositionalMultiphaseStatisticsAggregator.hpp + * @details Region statistics data is stored as follow: + + * Problem : ProblemManager + * |-> domain : DomainPartition + * |-> MeshBodies : Group + * |-> cartesianMesh : MeshBody + * |-> meshLevels : Group + * |-> Level0 : MeshLevel + * | |-> nodeManager : NodeManager + * | | |-> sets : Group + * | | | * all : Wrapper< index array > + * | | | * xneg : Wrapper< index array > + * | | [...] (other element sets) + * | | + * | |-> ElementRegions : ElementRegionManager + * | | |-> Channel : CellElementRegion + * | | | |-> cb-0_0_0 : CellElementSubRegion + * | | | | | * pressure : Wrapper< real64 array > + * | | | | | * temperature : Wrapper< real64 array > + * | | | | [...] (other fields) + * | | | | + * | | | |-> cb-0_0_1 : CellElementSubRegion + * | | | | | * pressure : Wrapper< real64 array > + * | | | | | * temperature : Wrapper< real64 array > + * | | | | [...] (other fields) + * | | | | + * | | | [...] (other sub-regions) + * | | | + * | | |-> Barrier : CellElementRegion + * | | |-> cb-1_0_0 : CellElementSubRegion + * | | |-> cb-1_0_1 : CellElementSubRegion + * | | [...] (other sub-regions) + * | | + * | [...] (other element managers) + * ____ | | + * | | |-> statistics : Group (storage for all stats) + * | | |-> compFlowStats : Group (storage for this instance stats) + * | | | |-> cflStatistics : CFLStatistics + * | | | |-> regionsStatistics : RegionStatistics (aggregate) + * | | | |-> Channel : RegionStatistics (aggregate, mpi reduced) + * | | | | |-> cb-0_0_0 : RegionStatistics (compute read-back) + * stats | | | | |-> cb-0_0_1 : RegionStatistics (compute read-back) + * data -> | | | | [...] (other sub-regions stats) + * | | | | + * | | | |-> Barrier : RegionStatistics (aggregate, mpi reduced) + * | | | |-> cb-1_0_0 : RegionStatistics (compute read-back) + * | | | |-> cb-1_0_1 : RegionStatistics (compute read-back) + * | | | [...] (other sub-regions stats) + * | | | + * |___ | [...] (other stats storages) + * | + * [...] (other discretizations) + */ + +#ifndef SRC_CORECOMPONENTS_PHYSICSSOLVERS_FLUIDFLOW_COMPOSITIONALMULTIPHASESTATISTICSAGGREGATOR_HPP_ +#define SRC_CORECOMPONENTS_PHYSICSSOLVERS_FLUIDFLOW_COMPOSITIONALMULTIPHASESTATISTICSAGGREGATOR_HPP_ + +#include "common/DataTypes.hpp" +#include "physicsSolvers/StatisticsAggregatorBase.hpp" + +namespace geos +{ + +class CompositionalMultiphaseBase; + +namespace compositionalMultiphaseStatistics +{ +class StatsAggregator; +class RegionStatistics; +} + +template<> +struct StatsAggregatorTraits< compositionalMultiphaseStatistics::StatsAggregator > +{ + using SolverType = CompositionalMultiphaseBase; + using StatsGroupType = compositionalMultiphaseStatistics::RegionStatistics; +}; + +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 +}; + +/** + * @brief Output data group to contain the result of a given stat aggregator on the dataRepository. + * Attributes are public since the class is a POD. + * @todo repair 1D HDF5 outputs to enable stats HDF5 outputs + */ +class RegionStatistics : public RegionStatisticsBase +{ +public: + + /// average region pressure (numerator value before postAggregateCompute()) + real64 m_averagePressure; + /// minimum region pressure + real64 m_minPressure; + /// maximum region pressure + real64 m_maxPressure; + + /// minimum region delta pressure + real64 m_minDeltaPressure; + /// maximum region delta pressure + real64 m_maxDeltaPressure; + + /// average region temperature (numerator value before postAggregateCompute()) + real64 m_averageTemperature; + /// minimum region temperature + real64 m_minTemperature; + /// maximum region temperature + real64 m_maxTemperature; + + /// total region pore volume + real64 m_totalPoreVolume; + /// total region uncompacted pore volume (not necessarily output, useful for weighting cell pressure data) + real64 m_totalUncompactedPoreVolume; + /// phase region dynamic pore volume + array1d< real64 > const m_phaseDynamicPoreVolume; + + /// region phase mass (trapped and non-trapped, immobile and mobile) + array1d< real64 > const m_phaseMass; + /// trapped region phase mass + array1d< real64 > const m_trappedPhaseMass; + /// non-trapped region phase mass (available after postAggregateCompute()) + array1d< real64 > const m_nonTrappedPhaseMass; + /// immobile region phase mass + array1d< real64 > const m_immobilePhaseMass; + /// mobile region phase mass (available after postAggregateCompute()) + array1d< real64 > const m_mobilePhaseMass; + /// region component mass + array2d< real64 > const m_componentMass; + + // TODO? -> split to struct PressureStats...MassStats: + // - optional computation of each stats + // - VKS for struct name ("pressureStats"..."massStats") + // - current RegionStatistics struct bits + + /** + * @brief Construct a new Region Statistics object + * @param targetName name of the data-repository object that is targeted by the statistics + * (mesh level / region / sub-region). + * @param parent the instance parent in data-repository + * @param numPhases Fluid phase count + * @param numComponents Fluid component count + */ + RegionStatistics( string const & targetName, dataRepository::Group * const parent, + integer numPhases, integer numComponents ); + + RegionStatistics( RegionStatistics && ) = default; + +}; + +/** + * @brief Output data group to contain the result of a given stat aggregator on the dataRepository. + * Attributes are public since the class is a POD. Can it be replaced by a wrapped-struct? + */ +class CFLStatistics : public RegionStatisticsBase +{ +public: + + /// Maximum Courant Friedrichs Lewy number in the grid for each phase + real64 m_maxPhaseCFL; + + /// Maximum Courant-Friedrichs-Lewy number in the grid for each component + real64 m_maxCompCFL; + + /** + * @brief Construct a new CFLStatistics object + * @param name instance name in data-repository + * @param parent the instance parent in data-repository + */ + CFLStatistics( const string & name, dataRepository::Group * const parent ); +}; + +/** + * @brief Reponsible of computing physical statistics over the grid, registering the result in the + * data repository, but not storing / outputing it by itself. It does not have mutable state + * except the encountered issues. + */ +class StatsAggregator : public StatsAggregatorBase< StatsAggregator > +{ +public: + + using Base = StatsAggregatorBase< StatsAggregator >; + + /** + * @brief the associated view keys + */ + struct ViewKeys + { + /// String for the cfl statistics group + constexpr static char const * cflStatisticsString() { return "cflStatistics"; } + }; + + /** + * @brief Construct a new Stats Aggregator object + * @param ownerName the unique name of the entity requesting the statistics. + * An error is thrown if not unique in this context. + */ + StatsAggregator( dataRepository::DataContext const & ownerDataContext ); + + /** + * @brief Enable the computation of any statistics, initialize data structure to collect them. + * Register the resulting data wrappers so they will be targeted by TimeHistory output + * @param solver flow solver object to retrieve: + - the simulated regions, + - fields for statistics computation. + * @param meshBodies The Group containing the MeshBody objects + */ + void initStatisticsAggregation( dataRepository::Group & meshBodies, + CompositionalMultiphaseBase & solver ); + + /** + * @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 or after the "registerDataOnMesh" initialization phase + * @param meshBodies The Group containing the MeshBody objects + */ + void enableRegionStatisticsAggregation(); + + /** + * @brief Register the results structs & wrappers so they will be targeted by TimeHistory output + * @note Must be called in or after the "registerDataOnMesh" initialization phase + * @param meshBodies The Group containing the MeshBody objects + */ + void enableCFLStatistics(); + + /** + * @brief set the statistics as dirty, ensuring isComputed() will be false until the next computation. + */ + void setDirty(); + + /** + * @brief Compute CFL numbers + * @param[in] time current time + * @param[in] dt the time step size + * @param[in] domain the domain partition + * @return false if there was a problem that prevented the statistics to be computed correctly. + */ + bool computeCFLNumbers( real64 const time, + real64 const dt, + DomainPartition & domain ); + + integer getNumPhases() const + { return m_numPhases; } + + integer getNumComponents() const + { return m_numComponents; } + + CFLStatistics * getCFLStatistics( DomainPartition & domain ) const; + + CFLStatistics & getCflStatistics( MeshLevel & mesh ) const; + + // template implementations + /// @cond DO_NOT_DOCUMENT + + void initStats( RegionStatistics & stats, real64 time ) const; + void computeSubRegionRankStats( CellElementSubRegion & subRegion, RegionStatistics & subRegionStats ) const; + void aggregateStats( RegionStatistics & stats, RegionStatistics const & other ) const; + void mpiAggregateStats( RegionStatistics & stats ) const; + void postAggregateStats( RegionStatistics & stats ); + + /// @endcond + +private: + + SolverType * m_solver = nullptr; + + AggregatorParameters m_params; + + StatsState m_cflStatsState; + + integer m_numPhases; + + integer m_numComponents; + +}; + +} /* 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..e6409a74d8c --- /dev/null +++ b/src/coreComponents/physicsSolvers/fluidFlow/CompositionalMultiphaseStatisticsTask.cpp @@ -0,0 +1,418 @@ +/* + * ------------------------------------------------------------------------------------------------------------ + * 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 "constitutive/fluid/multifluid/MultiFluidBase.hpp" +#include "physicsSolvers/LogLevelsInfo.hpp" +#include "physicsSolvers/fluidFlow/LogLevelsInfo.hpp" +#include "physicsSolvers/fluidFlow/CompositionalMultiphaseHybridFVM.hpp" + + +namespace geos +{ + +using namespace dataRepository; + +namespace compositionalMultiphaseStatistics +{ + +StatsTask::StatsTask( string const & 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(); + + GEOS_THROW_IF_EQ_MSG( m_solver, nullptr, + "To identify simulated regions, a solver must be provided.", + InputError, getWrapperDataContext( getSolverWrapperKey() ) ); + + if( !dynamicCast< CompositionalMultiphaseBase * >( m_solver ) ) + { + GEOS_THROW( "Incompatible solver selected, a compositional multiphase solver is expected", + InputError, getDataContext() ); + } + else if( dynamicCast< CompositionalMultiphaseHybridFVM * >( m_solver ) && m_computeCFLNumbers != 0 ) + { + GEOS_THROW( "The option to compute CFL numbers is incompatible with CompositionalMultiphaseHybridFVM", + InputError, getDataContext() ); + } + + m_aggregator = std::make_unique< StatsAggregator >( getDataContext() ); +} + +void StatsTask::registerDataOnMesh( Group & meshBodies ) +{ + // for now, this guard is needed to avoid breaking the xml schema generation + if( m_solver == nullptr || m_aggregator == nullptr ) + return; + + prepareFluidMetaData(); + + if( m_computeRegionStatistics || m_computeCFLNumbers ) + { + // expected to work as check is done in postInputInitialization() + CompositionalMultiphaseBase * castedSolver = dynamicCast< CompositionalMultiphaseBase * >( m_solver ); + GEOS_ERROR_IF_EQ_MSG( castedSolver, nullptr, + GEOS_FMT( "{} {}: Unexpected error (solver pointer changed?)", catalogName(), getDataContext() ) ); + m_aggregator->initStatisticsAggregation( meshBodies, *castedSolver ); + } + else + { + GEOS_WARNING( GEOS_FMT( "{} {}: No computing option enabled, no output is scheduled.", + catalogName(), getDataContext() ) ); + } + + if( m_computeRegionStatistics ) + m_aggregator->enableRegionStatisticsAggregation(); + + // if we have to compute CFL numbers later, we need to register additional variables + if( m_computeCFLNumbers ) + m_aggregator->enableCFLStatistics(); + + m_aggregator->forRegionStatistics( [&] ( MeshLevel & mesh, RegionStatistics & ) + { + prepareLogTableLayouts( mesh.getName() ); + prepareCsvTableLayouts( mesh.getName() ); + } ); +} + +void StatsTask::prepareFluidMetaData() +{ + using namespace constitutive; + + 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, string() ) ); + + 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_logFormatters.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( GEOS_FMT( "Time [{}]", units::getSymbol( units::Unit::Time ))), + TableLayout::Column( "Region" ), // TODO : mention this change in PR description + TableLayout::Column( GEOS_FMT( "Min pressure [{}]", units::getSymbol( units::Unit::Pressure ))), + TableLayout::Column( GEOS_FMT( "Average pressure [{}]", units::getSymbol( units::Unit::Pressure )) ), + TableLayout::Column( GEOS_FMT( "Max pressure [{}]", units::getSymbol( units::Unit::Pressure ) ) ), + TableLayout::Column( GEOS_FMT( "Min delta pressure [{}]", units::getSymbol( units::Unit::Pressure ))), + TableLayout::Column( GEOS_FMT( "Max delta pressure [{}]", units::getSymbol( units::Unit::Pressure ))), + TableLayout::Column( GEOS_FMT( "Min temperature [{}]", units::getSymbol( units::Unit::Temperature ) )), + TableLayout::Column( GEOS_FMT( "Average temperature [{}]", units::getSymbol( units::Unit::Temperature ) )), + TableLayout::Column( GEOS_FMT( "Max temperature [{}]", units::getSymbol( units::Unit::Temperature ) )), + TableLayout::Column( 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 ); + 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 = m_csvFormatters.get_inserted( string( meshName ) ); + csvFormatter = std::move( std::make_unique< TableCSVFormatter >( tableLayout ) ); + + // output CSV header + std::ofstream outputFile( getCsvFileName( meshName ) ); + outputFile << csvFormatter->headerToString(); + GEOS_LOG( GEOS_FMT( "table {} : {}", meshName, csvFormatter->headerToString() ) ); // TODO : remove this log +} + +string StatsTask::getCsvFileName( string_view meshName ) const +{ return GEOS_FMT( "{}/{}.csv", m_outputDir, meshName ); } + +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; + + GEOS_ERROR_IF( !m_aggregator, + "No statistics aggregator initialized!", getDataContext() ); + + m_aggregator->computeRegionsStatistics( statsTime ); + + m_aggregator->forRegionStatistics( [&] ( MeshLevel & mesh, RegionStatistics & meshRegionsStatistics ) + { + if( m_computeRegionStatistics ) + { + outputLogStats( statsTime, mesh, meshRegionsStatistics ); + outputCsvStats( statsTime, mesh, meshRegionsStatistics ); + } + } ); + + if( m_computeCFLNumbers ) + m_aggregator->computeCFLNumbers( statsTime, dt, domain ); + + return false; +} + +void StatsTask::outputLogStats( real64 const statsTime, + MeshLevel & mesh, + RegionStatistics & meshRegionsStatistics ) +{ + if( MpiWrapper::commRank() > 0 || !isLogLevelActive< logInfo::Statistics >( this->getLogLevel() ) ) + return; + + auto const formatterIter = m_logFormatters.find( mesh.getName() ); + if( formatterIter==m_logFormatters.end()) + return; + + TableTextFormatter const & formatter = *formatterIter->second; + TableData tableData; + static constexpr auto merge = CellType::MergeNext; + + 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 = units::getSymbol( units::ReservoirVolume ); + + tableData.getErrorsList().appendErrors( m_aggregator->getWarnings() ); + + tableData.addRow( "Statistics time", merge, merge, statsTime ); + + // lamda to apply for each region statistics + auto const outputRegionStats = [&] ( string_view targetName, RegionStatistics & stats ) + { + tableData.addSeparator(); + tableData.addRow( merge, merge, merge, "" ); + + tableData.addSeparator(); + tableData.addRow( merge, merge, merge, GEOS_FMT( "Region '{}'", targetName ) ); + tableData.addSeparator(); + tableData.addRow( "statistics", "min", "average", "max" ); + tableData.addSeparator(); + + tableData.addRow( GEOS_FMT( "Pressure [{}]", pressureUnit ), + stats.m_minPressure, stats.m_averagePressure, stats.m_maxPressure ); + tableData.addRow( GEOS_FMT( "Delta pressure [{}]", pressureUnit ), + stats.m_minDeltaPressure, "/", stats.m_maxDeltaPressure ); + tableData.addRow( GEOS_FMT( "Temperature [{}]", tempUnit ), + stats.m_minTemperature, stats.m_averageTemperature, stats.m_maxTemperature ); + + tableData.addSeparator(); + + tableData.addRow( GEOS_FMT( "Total dynamic pore volume [{}]", resVolUnit ), + "all", + CellType::MergeNext, + stats.m_totalPoreVolume ); + tableData.addRow( GEOS_FMT( "Phase dynamic pore volume [{}]", resVolUnit ), + stringutilities::joinLambda( m_fluid.m_phaseNames, "\n", []( auto data ) { return data[0]; } ), + CellType::MergeNext, + stringutilities::joinLambda( stats.m_phaseDynamicPoreVolume, "\n", []( auto data ) { return data[0]; } ) ); + + tableData.addSeparator(); + + tableData.addRow( GEOS_FMT( "Phase mass [{}]", massUnit ), + stringutilities::joinLambda( m_fluid.m_phaseNames, "\n", []( auto data ) { return data[0]; } ), + CellType::MergeNext, + stringutilities::joinLambda( stats.m_phaseMass, "\n", []( auto data ) { return data[0]; } ) ); + + tableData.addSeparator(); + + tableData.addRow( GEOS_FMT( "Trapped phase mass (metric 1) [{}]", massUnit ), + stringutilities::joinLambda( m_fluid.m_phaseNames, "\n", []( auto value ) { return value[0]; } ), + CellType::MergeNext, + stringutilities::joinLambda( stats.m_trappedPhaseMass, "\n", []( auto value ) { return value[0]; } ) ); + tableData.addRow( GEOS_FMT( "Non-trapped phase mass (metric 1) [{}]", massUnit ), + stringutilities::joinLambda( m_fluid.m_phaseNames, "\n", []( auto value ) { return value[0]; } ), + CellType::MergeNext, + stringutilities::joinLambda( stats.m_nonTrappedPhaseMass, "\n", []( auto value ) { return value[0]; } ) ); + + tableData.addSeparator(); + + tableData.addRow( GEOS_FMT( "Immobile phase mass (metric 2) [{}]", massUnit ), + stringutilities::joinLambda( m_fluid.m_phaseNames, "\n", []( auto value ) { return value[0]; } ), + CellType::MergeNext, + stringutilities::joinLambda( stats.m_immobilePhaseMass, "\n", []( auto value ) { return value[0]; } ) ); + tableData.addRow( GEOS_FMT( "Mobile phase mass (metric 2) [{}]", massUnit ), + stringutilities::joinLambda( m_fluid.m_phaseNames, "\n", []( auto value ) { return value[0]; } ), + CellType::MergeNext, + stringutilities::joinLambda( stats.m_mobilePhaseMass, "\n", []( auto value ) { return value[0]; } ) ); + + tableData.addSeparator(); + + tableData.addRow( GEOS_FMT( "Component mass [{}]", massUnit ), + "TODO" /*stringutilities::join( m_fluid.m_phaseCompNames, '\n' )*/, + CellType::MergeNext, + stringutilities::join( stats.m_componentMass, '\n' ) ); + }; + + // apply the lambda for each region and, finally, the mesh summary + m_aggregator->forRegionStatistics( mesh, meshRegionsStatistics, + [&] ( CellElementRegion & region, RegionStatistics & stats ) + { + outputRegionStats( region.getName(), stats ); + } ); + outputRegionStats( mesh.getName(), meshRegionsStatistics ); + + // output to log + GEOS_LOG_RANK_0( formatter.toString( tableData ) ); +} + +void StatsTask::outputCsvStats( real64 statsTime, + MeshLevel & mesh, + RegionStatistics & meshRegionsStatistics ) +{ + if( MpiWrapper::commRank() > 0 || m_writeCSV == 0 ) + return; + + auto const formatterIter = m_csvFormatters.find( mesh.getName() ); + if( formatterIter==m_csvFormatters.end()) + return; + + TableCSVFormatter const & formatter = *formatterIter->second; + TableData tableData; + + stdVector< string > row; + row.reserve( formatter.getLayout().getTotalLowermostColumnCount() ); + + // lamda to apply for each region statistics + auto const outputRegionStats = [&] ( string_view targetName, RegionStatistics & stats ) + { + + auto addPhaseValues = []( auto & list, auto const & values ) + { + for( auto value : values ) + list.emplace_back( std::to_string( value ) ); + }; + + row.clear(); + row.insert( row.begin(), + { std::to_string( statsTime ), + string( targetName ), + std::to_string( stats.m_minPressure ), + std::to_string( stats.m_averagePressure ), + std::to_string( stats.m_maxPressure ), + std::to_string( stats.m_minDeltaPressure ), + std::to_string( stats.m_maxDeltaPressure ), + std::to_string( stats.m_minTemperature ), + std::to_string( stats.m_averageTemperature ), + std::to_string( stats.m_maxTemperature ), + std::to_string( stats.m_totalPoreVolume ), + } ); + addPhaseValues( row, stats.m_phaseDynamicPoreVolume ); + addPhaseValues( row, stats.m_phaseMass ); + addPhaseValues( row, stats.m_trappedPhaseMass ); + addPhaseValues( row, stats.m_nonTrappedPhaseMass ); + addPhaseValues( row, stats.m_immobilePhaseMass ); + addPhaseValues( row, stats.m_mobilePhaseMass ); + addPhaseValues( row, stats.m_componentMass ); // TODO verify phase / comp ordering + + tableData.addRow( row ); + }; + + // apply the lambda for each region and, finally, the mesh summary + m_aggregator->forRegionStatistics( mesh, meshRegionsStatistics, + [&] ( CellElementRegion & region, RegionStatistics & stats ) + { + outputRegionStats( region.getName(), stats ); + } ); + outputRegionStats( mesh.getName(), meshRegionsStatistics ); + + // append to csv file + std::ofstream outputFile( getCsvFileName( mesh.getName() ), std::ios_base::app ); + outputFile << formatter.dataToString( tableData ); + outputFile.close(); +} + +REGISTER_CATALOG_ENTRY( TaskBase, + StatsTask, + string const &, dataRepository::Group * const ) + +} /* namespace compositionalMultiphaseStatistics */ + +} /* namespace geos */ diff --git a/src/coreComponents/physicsSolvers/fluidFlow/CompositionalMultiphaseStatistics.hpp b/src/coreComponents/physicsSolvers/fluidFlow/CompositionalMultiphaseStatisticsTask.hpp similarity index 55% rename from src/coreComponents/physicsSolvers/fluidFlow/CompositionalMultiphaseStatistics.hpp rename to src/coreComponents/physicsSolvers/fluidFlow/CompositionalMultiphaseStatisticsTask.hpp index dbe275b7b6f..6c337d49273 100644 --- a/src/coreComponents/physicsSolvers/fluidFlow/CompositionalMultiphaseStatistics.hpp +++ b/src/coreComponents/physicsSolvers/fluidFlow/CompositionalMultiphaseStatisticsTask.hpp @@ -14,25 +14,32 @@ */ /** - * @file CompositionalMultiphaseStatistics.hpp + * @file CompositionalMultiphaseStatisticsTask.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 "physicsSolvers/FieldStatisticsBase.hpp" +#include "physicsSolvers/fluidFlow/CompositionalMultiphaseStatisticsAggregator.hpp" +#include namespace geos { class CompositionalMultiphaseBase; +namespace compositionalMultiphaseStatistics +{ + /** - * @class CompositionalMultiphaseStatistics + * @class compositionalMultiphaseStatistics::StatsStats * * 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,12 @@ 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; - }; + StatsAggregator & getStatisticsAggregator() + { return *m_aggregator; } + + StatsAggregator const & getStatisticsAggregator() const + { return *m_aggregator; } + private: using Base = FieldStatisticsBase< CompositionalMultiphaseBase >; @@ -117,37 +88,38 @@ 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; + void registerDataOnMesh( Group & meshBodies ) 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 prepareFluidMetaData(); - /** - * @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 prepareLogTableLayouts( string_view tableName ); - void postInputInitialization() override; + void prepareCsvTableLayouts( string_view tableName ); - void registerDataOnMesh( Group & meshBodies ) override; + string getCsvFileName( string_view meshName ) const; + + void outputLogStats( real64 statsTime, + MeshLevel & mesh, + RegionStatistics & meshRegionsStatistics ); + + void outputCsvStats( real64 statsTime, + MeshLevel & mesh, + RegionStatistics & meshRegionsStatistics ); + + /// For each discretization (MeshLevel name), table formatter for log output. + stdMap< string, std::unique_ptr< TableTextFormatter > > m_logFormatters; + + /// For each discretization (MeshLevel name), table formatter for csv output. + stdMap< string, std::unique_ptr< TableCSVFormatter > > m_csvFormatters; + + // mesh statistics aggregator + std::unique_ptr< StatsAggregator > m_aggregator; /// Flag to decide whether CFL numbers are computed or not integer m_computeCFLNumbers; @@ -158,9 +130,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_ */ diff --git a/src/coreComponents/physicsSolvers/fluidFlow/SinglePhaseStatistics.cpp b/src/coreComponents/physicsSolvers/fluidFlow/SinglePhaseStatistics.cpp deleted file mode 100644 index 31cfee04b70..00000000000 --- a/src/coreComponents/physicsSolvers/fluidFlow/SinglePhaseStatistics.cpp +++ /dev/null @@ -1,303 +0,0 @@ -/* - * ------------------------------------------------------------------------------------------------------------ - * 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 SinglePhaseStatistics.cpp - */ - -#include "SinglePhaseStatistics.hpp" - -#include "mesh/DomainPartition.hpp" -#include "physicsSolvers/LogLevelsInfo.hpp" -#include "physicsSolvers/fluidFlow/SinglePhaseBase.hpp" -#include "physicsSolvers/fluidFlow/FlowSolverBaseFields.hpp" -#include "physicsSolvers/fluidFlow/SinglePhaseBaseFields.hpp" -#include "physicsSolvers/fluidFlow/kernels/singlePhase/StatisticsKernel.hpp" -#include "common/format/table/TableData.hpp" -#include "common/format/table/TableFormatter.hpp" -#include "common/format/table/TableLayout.hpp" - -namespace geos -{ - -using namespace constitutive; -using namespace fields; -using namespace dataRepository; - -SinglePhaseStatistics::SinglePhaseStatistics( const string & name, - Group * const parent ): - Base( name, parent ) -{ - addLogLevel< logInfo::Statistics >(); -} - -void SinglePhaseStatistics::registerDataOnMesh( 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; - } - - m_solver->forDiscretizationOnMeshTargets( meshBodies, [&] ( string const &, - MeshLevel & mesh, - string_array const & regionNames ) - { - ElementRegionManager & elemManager = mesh.getElemManager(); - - 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() } ); - - if( m_writeCSV > 0 && MpiWrapper::commRank() == 0 ) - { - TableLayout tableLayout( { - TableLayout::Column().setName( GEOS_FMT( "Time [{}]", units::getSymbol( units::Unit::Time ))), - 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 ) )), - TableLayout::Column().setName( GEOS_FMT( "Total fluid mass [{}]", units::getSymbol( units::Unit::Mass ))) - } ); - - TableCSVFormatter csvFormatter( tableLayout ); - std::ofstream outputFile( m_outputDir + "/" + regionNames[i] + ".csv" ); - outputFile << csvFormatter.headerToString(); - outputFile.close(); - } - } - } ); -} - -bool SinglePhaseStatistics::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 ) -{ - m_solver->forDiscretizationOnMeshTargets( domain.getMeshBodies(), [&] ( string const &, - MeshLevel & mesh, - string_array const & regionNames ) - { - // current time is time_n + dt - computeRegionStatistics( time_n + dt, mesh, regionNames ); - } ); - return false; -} - -void SinglePhaseStatistics::computeRegionStatistics( real64 const time, - MeshLevel & mesh, - string_array const & regionNames ) const -{ - GEOS_MARK_FUNCTION; - // Step 1: initialize the average/min/max quantities - ElementRegionManager & elemManager = mesh.getElemManager(); - for( size_t i = 0; i < regionNames.size(); ++i ) - { - ElementRegionBase & region = elemManager.getRegion( regionNames[i] ); - RegionStatistics & stats = region.getReference< RegionStatistics >( viewKeyStruct::regionStatisticsString() ); - - stats.averagePressure = 0.0; - stats.maxPressure = -LvArray::NumericLimits< real64 >::max; - stats.minPressure = LvArray::NumericLimits< real64 >::max; - - stats.maxDeltaPressure = -LvArray::NumericLimits< real64 >::max; - stats.minDeltaPressure = LvArray::NumericLimits< real64 >::max; - - stats.averageTemperature = 0.0; - stats.maxTemperature = -LvArray::NumericLimits< real64 >::max; - stats.minTemperature = LvArray::NumericLimits< real64 >::max; - - stats.totalPoreVolume = 0.0; - stats.totalUncompactedPoreVolume = 0.0; - stats.totalMass = 0.0; - } - - // Step 2: increment the average/min/max quantities for all the subRegions - elemManager.forElementSubRegions( regionNames, [&]( localIndex const, - ElementSubRegionBase & subRegion ) - { - - 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 deltaPres = subRegion.getField< flow::deltaPressure >(); - arrayView1d< real64 const > const temp = subRegion.getField< flow::temperature >(); - - string const & solidName = subRegion.getReference< string >( SinglePhaseBase::viewKeyStruct::solidNamesString() ); - Group const & constitutiveModels = subRegion.getGroup( ElementSubRegionBase::groupKeyStruct::constitutiveModelsString() ); - CoupledSolidBase const & solid = constitutiveModels.getGroup< CoupledSolidBase >( solidName ); - arrayView1d< real64 const > const refPorosity = solid.getReferencePorosity(); - arrayView2d< real64 const > const porosity = solid.getPorosity(); - - string const & fluidName = subRegion.template getReference< string >( FlowSolverBase::viewKeyStruct::fluidNamesString() ); - SingleFluidBase const & fluid = constitutiveModels.getGroup< SingleFluidBase >( fluidName ); - arrayView2d< real64 const, constitutive::singlefluid::USD_FLUID > const densities = fluid.density(); - - real64 subRegionAvgPresNumerator = 0.0; - real64 subRegionMinPres = 0.0; - real64 subRegionMaxPres = 0.0; - real64 subRegionMinDeltaPres = 0.0; - real64 subRegionMaxDeltaPres = 0.0; - real64 subRegionAvgTempNumerator = 0.0; - real64 subRegionMinTemp = 0.0; - real64 subRegionMaxTemp = 0.0; - real64 subRegionTotalUncompactedPoreVol = 0.0; - real64 subRegionTotalPoreVol = 0.0; - real64 subRegionTotalMass = 0.0; - - singlePhaseBaseKernels::StatisticsKernel:: - launch( subRegion.size(), - elemGhostRank, - volume, - pres, - deltaPres, - temp, - refPorosity, - porosity, - densities, - subRegionMinPres, - subRegionAvgPresNumerator, - subRegionMaxPres, - subRegionMinDeltaPres, - subRegionMaxDeltaPres, - subRegionMinTemp, - subRegionAvgTempNumerator, - subRegionMaxTemp, - subRegionTotalUncompactedPoreVol, - subRegionTotalPoreVol, - subRegionTotalMass ); - - ElementRegionBase & region = elemManager.getRegion( ElementRegionBase::getParentRegion( subRegion ).getName() ); - RegionStatistics & stats = region.getReference< RegionStatistics >( viewKeyStruct::regionStatisticsString() ); - - stats.averagePressure += subRegionAvgPresNumerator; - if( subRegionMinPres < stats.minPressure ) - { - stats.minPressure = subRegionMinPres; - } - if( subRegionMaxPres > stats.maxPressure ) - { - stats.maxPressure = subRegionMaxPres; - } - - if( subRegionMinDeltaPres < stats.minDeltaPressure ) - { - stats.minDeltaPressure = subRegionMinDeltaPres; - } - if( subRegionMaxDeltaPres > stats.maxDeltaPressure ) - { - stats.maxDeltaPressure = subRegionMaxDeltaPres; - } - - stats.averageTemperature += subRegionAvgTempNumerator; - if( subRegionMinTemp < stats.minTemperature ) - { - stats.minTemperature = subRegionMinTemp; - } - if( subRegionMaxTemp > stats.maxTemperature ) - { - stats.maxTemperature = subRegionMaxTemp; - } - - stats.totalUncompactedPoreVolume += subRegionTotalUncompactedPoreVol; - stats.totalPoreVolume += subRegionTotalPoreVol; - stats.totalMass += subRegionTotalMass; - } ); - - // Step 3: synchronize the results over the MPI ranks - for( size_t i = 0; i < regionNames.size(); ++i ) - { - ElementRegionBase & region = elemManager.getRegion( regionNames[i] ); - RegionStatistics & stats = region.getReference< RegionStatistics >( viewKeyStruct::regionStatisticsString() ); - - stats.minPressure = MpiWrapper::min( stats.minPressure ); - stats.averagePressure = MpiWrapper::sum( stats.averagePressure ); - stats.maxPressure = MpiWrapper::max( stats.maxPressure ); - - stats.minDeltaPressure = MpiWrapper::min( stats.minDeltaPressure ); - stats.maxDeltaPressure = MpiWrapper::max( stats.maxDeltaPressure ); - - stats.minTemperature = MpiWrapper::min( stats.minTemperature ); - stats.averageTemperature = MpiWrapper::sum( stats.averageTemperature ); - stats.maxTemperature = MpiWrapper::max( stats.maxTemperature ); - - stats.totalUncompactedPoreVolume = MpiWrapper::sum( stats.totalUncompactedPoreVolume ); - stats.totalPoreVolume = MpiWrapper::sum( stats.totalPoreVolume ); - stats.totalMass = MpiWrapper::sum( stats.totalMass ); - - if( stats.totalUncompactedPoreVolume > 0 ) - { - float invTotalUncompactedPoreVolume = 1.0 / stats.totalUncompactedPoreVolume; - stats.averagePressure *= invTotalUncompactedPoreVolume; - stats.averageTemperature *= invTotalUncompactedPoreVolume; - } - else - { - stats.averagePressure = 0.0; - stats.averageTemperature = 0.0; - GEOS_WARNING( GEOS_FMT( "{}: Cannot compute average pressure & temperature because region pore volume is zero.", regionNames[i] ), - getDataContext() ); - } - - string_view massUnit = units::getSymbol( m_solver->getMassUnit() ); - - if( isLogLevelActive< logInfo::Statistics >( this->getLogLevel())&& MpiWrapper::commRank() == 0 ) - { - TableData singPhaseStatsData; - singPhaseStatsData.addRow( "Pressure[Pa]", stats.minPressure, stats.averagePressure, stats.maxPressure ); - singPhaseStatsData.addRow( "Delta pressure [Pa]", stats.minDeltaPressure, "/", stats.maxDeltaPressure ); - singPhaseStatsData.addRow( "Temperature [K]", stats.minTemperature, stats.averageTemperature, stats.maxTemperature ); - singPhaseStatsData.addSeparator(); - - singPhaseStatsData.addRow( "Total dynamic pore volume [rm^3]", CellType::MergeNext, CellType::MergeNext, stats.totalPoreVolume ); - singPhaseStatsData.addSeparator(); - singPhaseStatsData.addRow( GEOS_FMT( "Total fluid mass [{}]", massUnit ), CellType::MergeNext, CellType::MergeNext, stats.totalMass ); - - string const title = GEOS_FMT( "{}, {} (time {} s):", getName(), regionNames[i], time ); - TableLayout const singPhaseStatsLayout( title, { "statistics", "min", "average", "max" } ); - TableTextFormatter tableFormatter( singPhaseStatsLayout ); - GEOS_LOG_RANK_0( tableFormatter.toString( singPhaseStatsData ) ); - - if( m_writeCSV > 0 && MpiWrapper::commRank() == 0 ) - { - std::ofstream outputFile( m_outputDir + "/" + regionNames[i] + ".csv", std::ios_base::app ); - outputFile << time << "," << stats.minPressure << "," << stats.averagePressure << "," << stats.maxPressure << "," << - stats.minDeltaPressure << "," << stats.maxDeltaPressure << "," << - stats.minTemperature << "," << stats.averageTemperature << "," << stats.maxTemperature << "," << - stats.totalPoreVolume << "," << stats.totalMass << std::endl; - outputFile.close(); - } - } - } -} - -REGISTER_CATALOG_ENTRY( TaskBase, - SinglePhaseStatistics, - string const &, dataRepository::Group * const ) - -} /* namespace geos */ diff --git a/src/coreComponents/physicsSolvers/fluidFlow/SinglePhaseStatisticsAggregator.cpp b/src/coreComponents/physicsSolvers/fluidFlow/SinglePhaseStatisticsAggregator.cpp new file mode 100644 index 00000000000..1ece65e1213 --- /dev/null +++ b/src/coreComponents/physicsSolvers/fluidFlow/SinglePhaseStatisticsAggregator.cpp @@ -0,0 +1,189 @@ +/* + * ------------------------------------------------------------------------------------------------------------ + * 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 SinglePhaseStatistics.cpp + */ + +#include "SinglePhaseStatisticsAggregator.hpp" + +#include "mesh/DomainPartition.hpp" +#include "physicsSolvers/LogLevelsInfo.hpp" +#include "physicsSolvers/StatisticsAggregatorBaseHelpers.hpp" +#include "physicsSolvers/fluidFlow/SinglePhaseBase.hpp" +#include "physicsSolvers/fluidFlow/FlowSolverBaseFields.hpp" +#include "physicsSolvers/fluidFlow/SinglePhaseBaseFields.hpp" +#include "physicsSolvers/fluidFlow/kernels/singlePhase/StatisticsKernel.hpp" +#include "common/format/table/TableData.hpp" +#include "common/format/table/TableFormatter.hpp" +#include "common/format/table/TableLayout.hpp" + +namespace geos +{ + +namespace singlePhaseStatistics +{ + +using namespace constitutive; +using namespace dataRepository; + +RegionStatistics::RegionStatistics( string const & name, dataRepository::Group * const parent ): + RegionStatisticsBase( name, parent ) +{} + +StatsAggregator::StatsAggregator( DataContext const & ownerDataContext ): + Base( ownerDataContext ) +{} + +void StatsAggregator::enableRegionStatisticsAggregation() +{ + auto const registerStats = [=] ( Group & parent, + string const & targetName ) -> RegionStatistics & + { + return parent.registerGroup( targetName, + std::make_unique< RegionStatistics >( targetName, &parent ) ); + }; + + Base::enableRegionStatisticsAggregation( registerStats ); +} + +void StatsAggregator::initStats( RegionStatistics & stats, real64 const time ) const +{ + stats.m_time = time; + + stats.m_averagePressure = 0.0; + stats.m_maxPressure = 0.0; + stats.m_minPressure = LvArray::NumericLimits< real64 >::max; + + stats.m_maxDeltaPressure = -LvArray::NumericLimits< real64 >::max; + stats.m_minDeltaPressure = LvArray::NumericLimits< real64 >::max; + + stats.m_averageTemperature = 0.0; + stats.m_maxTemperature = 0.0; + stats.m_minTemperature = LvArray::NumericLimits< real64 >::max; + + stats.m_totalDynamicPoreVolume = 0.0; + stats.m_totalUncompactedPoreVolume = 0.0; + + stats.m_totalMass = 0.0; +} + +void StatsAggregator::computeSubRegionRankStats( CellElementSubRegion & subRegion, + RegionStatistics & subRegionStats ) const +{ + static constexpr string_view solidNamesVK = SinglePhaseBase::viewKeyStruct::solidNamesString(); + static constexpr string_view fluidNamesVK = FlowSolverBase::viewKeyStruct::fluidNamesString(); + static constexpr string_view modelsVK = ElementSubRegionBase::groupKeyStruct::constitutiveModelsString(); + + arrayView1d< integer const > const elemGhostRank = subRegion.ghostRank(); + arrayView1d< real64 const > const volume = subRegion.getElementVolume(); + arrayView1d< real64 const > const pres = subRegion.getField< fields::flow::pressure >(); + arrayView1d< real64 const > const deltaPres = subRegion.getField< fields::flow::deltaPressure >(); + arrayView1d< real64 const > const temp = subRegion.getField< fields::flow::temperature >(); + + string const & solidName = subRegion.getReference< string >( string( solidNamesVK ) ); + Group const & constitutiveModels = subRegion.getGroup( string( modelsVK ) ); + CoupledSolidBase const & solid = constitutiveModels.getGroup< CoupledSolidBase >( solidName ); + arrayView1d< real64 const > const refPorosity = solid.getReferencePorosity(); + arrayView2d< real64 const > const porosity = solid.getPorosity(); + + string const & fluidName = subRegion.template getReference< string >( string( fluidNamesVK ) ); + SingleFluidBase const & fluid = constitutiveModels.getGroup< SingleFluidBase >( fluidName ); + arrayView2d< real64 const, constitutive::singlefluid::USD_FLUID > const densities = fluid.density(); + + singlePhaseBaseKernels::StatisticsKernel::launch( subRegion.size(), + elemGhostRank, + volume, + pres, + deltaPres, + temp, + refPorosity, + porosity, + densities, + subRegionStats.m_minPressure, + subRegionStats.m_averagePressure, + subRegionStats.m_maxPressure, + subRegionStats.m_minDeltaPressure, + subRegionStats.m_maxDeltaPressure, + subRegionStats.m_minTemperature, + subRegionStats.m_averageTemperature, + subRegionStats.m_maxTemperature, + subRegionStats.m_totalUncompactedPoreVolume, + subRegionStats.m_totalDynamicPoreVolume, + subRegionStats.m_totalMass ); +} + + +void StatsAggregator::aggregateStats( RegionStatistics & stats, + RegionStatistics const & other ) const +{ + stats.m_averagePressure += other.m_averagePressure; + stats.m_minPressure = LvArray::math::min( stats.m_minPressure, other.m_minPressure ); + stats.m_maxPressure = LvArray::math::max( stats.m_maxPressure, other.m_maxPressure ); + + stats.m_minDeltaPressure = LvArray::math::min( stats.m_minDeltaPressure, other.m_minDeltaPressure ); + stats.m_maxDeltaPressure = LvArray::math::max( stats.m_maxDeltaPressure, other.m_maxDeltaPressure ); + + stats.m_averageTemperature += other.m_averageTemperature; + stats.m_minTemperature = LvArray::math::min( stats.m_minTemperature, other.m_minTemperature ); + stats.m_maxTemperature = LvArray::math::max( stats.m_maxTemperature, other.m_maxTemperature ); + + stats.m_totalUncompactedPoreVolume += other.m_totalUncompactedPoreVolume; + stats.m_totalDynamicPoreVolume += other.m_totalDynamicPoreVolume; + + stats.m_totalMass += other.m_totalMass; +} + +void StatsAggregator::mpiAggregateStats( RegionStatistics & stats ) const +{ + stats.m_averagePressure = MpiWrapper::sum( stats.m_averagePressure ); + stats.m_minPressure = MpiWrapper::min( stats.m_minPressure ); + stats.m_maxPressure = MpiWrapper::max( stats.m_maxPressure ); + + stats.m_minDeltaPressure = MpiWrapper::min( stats.m_minDeltaPressure ); + stats.m_maxDeltaPressure = MpiWrapper::max( stats.m_maxDeltaPressure ); + + stats.m_averageTemperature = MpiWrapper::sum( stats.m_averageTemperature ); + stats.m_minTemperature = MpiWrapper::min( stats.m_minTemperature ); + stats.m_maxTemperature = MpiWrapper::max( stats.m_maxTemperature ); + + stats.m_totalUncompactedPoreVolume = MpiWrapper::sum( stats.m_totalUncompactedPoreVolume ); + stats.m_totalDynamicPoreVolume = MpiWrapper::sum( stats.m_totalDynamicPoreVolume ); + + stats.m_totalMass = MpiWrapper::sum( stats.m_totalMass ); +} + +void StatsAggregator::postAggregateStats( RegionStatistics & stats ) +{ + if( stats.m_totalUncompactedPoreVolume > 0 ) + { + float invTotalUncompactedPoreVolume = 1.0 / stats.m_totalUncompactedPoreVolume; + stats.m_averagePressure *= invTotalUncompactedPoreVolume; + stats.m_averageTemperature *= invTotalUncompactedPoreVolume; + } + else + { + stats.m_averagePressure = 0.0; + stats.m_averageTemperature = 0.0; + m_warnings.emplace_back( GEOS_FMT( "Cannot compute average pressure for '{}' because pore volume is zero in '{}'.", + getOwnerName(), stats.getTargetName() ) ); + } +} + +} /* namespace singlePhaseStatistics */ + +template class StatsAggregatorBase< singlePhaseStatistics::StatsAggregator >; + +} /* namespace geos */ diff --git a/src/coreComponents/physicsSolvers/fluidFlow/SinglePhaseStatisticsAggregator.hpp b/src/coreComponents/physicsSolvers/fluidFlow/SinglePhaseStatisticsAggregator.hpp new file mode 100644 index 00000000000..28b55442c8b --- /dev/null +++ b/src/coreComponents/physicsSolvers/fluidFlow/SinglePhaseStatisticsAggregator.hpp @@ -0,0 +1,197 @@ +/* + * ------------------------------------------------------------------------------------------------------------ + * 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 SinglePhaseStatisticsAggregator.hpp + * @details Region statistics data is stored as follow: + + * Problem : ProblemManager + * |-> domain : DomainPartition + * |-> MeshBodies : Group + * |-> cartesianMesh : MeshBody + * |-> meshLevels : Group + * |-> Level0 : MeshLevel + * | |-> nodeManager : NodeManager + * | | |-> sets : Group + * | | | * all : Wrapper< index array > + * | | | * xneg : Wrapper< index array > + * | | [...] (other element sets) + * | | + * | |-> ElementRegions : ElementRegionManager + * | | |-> Channel : CellElementRegion + * | | | |-> cb-0_0_0 : CellElementSubRegion + * | | | | | * pressure : Wrapper< real64 array > + * | | | | | * temperature : Wrapper< real64 array > + * | | | | [...] (other fields) + * | | | | + * | | | |-> cb-0_0_1 : CellElementSubRegion + * | | | | | * pressure : Wrapper< real64 array > + * | | | | | * temperature : Wrapper< real64 array > + * | | | | [...] (other fields) + * | | | | + * | | | [...] (other sub-regions) + * | | | + * | | |-> Barrier : CellElementRegion + * | | |-> cb-1_0_0 : CellElementSubRegion + * | | |-> cb-1_0_1 : CellElementSubRegion + * | | [...] (other sub-regions) + * | | + * | [...] (other element managers) + * ____ | | + * | | |-> statistics : Group (storage for all stats) + * | | |-> flowStats : Group (storage for this instance stats) + * | | | |-> regionsStatistics : RegionStatistics (aggregate) + * | | | |-> Channel : RegionStatistics (aggregate, mpi reduced) + * | | | | |-> cb-0_0_0 : RegionStatistics (compute read-back) + * stats | | | | |-> cb-0_0_1 : RegionStatistics (compute read-back) + * data -> | | | | [...] (other sub-regions stats) + * | | | | + * | | | |-> Barrier : RegionStatistics (aggregate, mpi reduced) + * | | | |-> cb-1_0_0 : RegionStatistics (compute read-back) + * | | | |-> cb-1_0_1 : RegionStatistics (compute read-back) + * | | | [...] (other sub-regions stats) + * | | | + * |___ | [...] (other stats storages) + * | + * [...] (other discretizations) + */ + +#ifndef SRC_CORECOMPONENTS_PHYSICSSOLVERS_FLUIDFLOW_SINGLEPHASESTATISTICSAGGREGATOR_HPP_ +#define SRC_CORECOMPONENTS_PHYSICSSOLVERS_FLUIDFLOW_SINGLEPHASESTATISTICSAGGREGATOR_HPP_ + +#include "physicsSolvers/StatisticsAggregatorBase.hpp" + +namespace geos +{ + +class SinglePhaseBase; + +namespace singlePhaseStatistics +{ +class StatsAggregator; +class RegionStatistics; +} + +template<> +struct StatsAggregatorTraits< singlePhaseStatistics::StatsAggregator > +{ + using SolverType = SinglePhaseBase; + using StatsGroupType = singlePhaseStatistics::RegionStatistics; +}; + +namespace singlePhaseStatistics +{ + +/** + * @brief Output data group to contain the result of a given stat aggregator on the dataRepository. + * Attributes are public since the class is a POD. + * @todo repair 1D HDF5 outputs to enable stats HDF5 outputs + */ +class RegionStatistics : public RegionStatisticsBase +{ +public: + + /// Time of statistics computation + real64 m_time; + + /// average region pressure (numerator value before postAggregateCompute()) + real64 m_averagePressure; + /// minimum region pressure + real64 m_minPressure; + /// maximum region pressure + real64 m_maxPressure; + + /// minimum region delta pressure + real64 m_minDeltaPressure; + /// maximum region delta pressure + real64 m_maxDeltaPressure; + + /// average region temperature (numerator value before postAggregateCompute()) + real64 m_averageTemperature; + /// minimum region temperature + real64 m_minTemperature; + /// maximum region temperature + real64 m_maxTemperature; + + /// total region pore volume + real64 m_totalDynamicPoreVolume; + /// total region uncompacted pore volume (not necessarily output, useful for weighting cell pressure data) + real64 m_totalUncompactedPoreVolume; + + // fluid mass + real64 m_totalMass; + + // TODO? -> split to struct PressureStats...MassStats: + // - optional computation of each stats + // - VKS for struct name ("pressureStats"..."massStats") + // - current RegionStatistics struct bits + + /** + * @brief Construct a new Region Statistics object + * @param targetName name of the data-repository object that is targeted by the statistics + * (mesh level / region / sub-region). + * @param parent the instance parent in data-repository + */ + RegionStatistics( string const & targetName, dataRepository::Group * const parent ); + + RegionStatistics( RegionStatistics && ) = default; + +}; + +/** + * @brief Reponsible of computing physical statistics over the grid, registering the result in the + * data repository, but not storing / outputing it by itself. It does not have mutable state + * except the encountered issues. + * @todo repair 1D HDF5 outputs to enable stats HDF5 outputs + */ +class StatsAggregator : public StatsAggregatorBase< StatsAggregator > +{ +public: + + using Base = StatsAggregatorBase< StatsAggregator >; + + /** + * @brief Construct a new Stats Aggregator object + * @param ownerName the unique name of the entity requesting the statistics. + * An error is thrown if not unique in this context. + */ + StatsAggregator( dataRepository::DataContext const & ownerDataContext ); + + /** + * @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 or after the "registerDataOnMesh" initialization phase + * @param meshBodies The Group containing the MeshBody objects + */ + void enableRegionStatisticsAggregation(); + + // template implementations + /// @cond DO_NOT_DOCUMENT + + void initStats( RegionStatistics & stats, real64 time ) const; + void computeSubRegionRankStats( CellElementSubRegion & subRegion, RegionStatistics & subRegionStats ) const; + void aggregateStats( RegionStatistics & stats, RegionStatistics const & other ) const; + void mpiAggregateStats( RegionStatistics & stats ) const; + void postAggregateStats( RegionStatistics & stats ); + + /// @endcond + +}; + +} /* namespace singlePhaseStatistics */ + +} /* namespace geos */ + +#endif /* SRC_CORECOMPONENTS_PHYSICSSOLVERS_FLUIDFLOW_SINGLEPHASESTATISTICSAGGREGATOR_HPP_ */ diff --git a/src/coreComponents/physicsSolvers/fluidFlow/SinglePhaseStatisticsTask.cpp b/src/coreComponents/physicsSolvers/fluidFlow/SinglePhaseStatisticsTask.cpp new file mode 100644 index 00000000000..3e918e6777f --- /dev/null +++ b/src/coreComponents/physicsSolvers/fluidFlow/SinglePhaseStatisticsTask.cpp @@ -0,0 +1,281 @@ +/* + * ------------------------------------------------------------------------------------------------------------ + * 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 SinglePhaseStatisticsTask.cpp + */ + +#include "SinglePhaseStatisticsTask.hpp" + +#include "physicsSolvers/LogLevelsInfo.hpp" +#include "physicsSolvers/fluidFlow/SinglePhaseBase.hpp" + +namespace geos +{ + +using namespace dataRepository; + +namespace singlePhaseStatistics +{ + +StatsTask::StatsTask( const string & name, Group * const parent ): + Base( name, parent ) +{ + addLogLevel< logInfo::Statistics >(); +} + +void StatsTask::postInputInitialization() +{ + Base::postInputInitialization(); + + GEOS_THROW_IF_EQ_MSG( m_solver, nullptr, + "To identify simulated regions, a solver must be provided.", + InputError, getWrapperDataContext( getSolverWrapperKey() ) ); + + if( !dynamicCast< SinglePhaseBase * >( m_solver ) ) + { + GEOS_THROW( "Incompatible solver selected, a single-phase solver is expected", + InputError, getDataContext() ); + } + + m_aggregator = std::make_unique< StatsAggregator >( getDataContext() ); +} + +void StatsTask::registerDataOnMesh( Group & meshBodies ) +{ + // for now, this guard is needed to avoid breaking the xml schema generation + if( m_solver == nullptr || m_aggregator == nullptr ) + return; + + if( m_writeCSV || isLogLevelActive< logInfo::Statistics >( this->getLogLevel()) ) + { + // expected to work as check is done in postInputInitialization() + SinglePhaseBase * castedSolver = dynamicCast< SinglePhaseBase * >( m_solver ); + GEOS_ERROR_IF_EQ_MSG( castedSolver, nullptr, + GEOS_FMT( "{} {}: Unexpected error (solver pointer changed?)", catalogName(), getDataContext() ) ); + m_aggregator->initStatisticsAggregation( meshBodies, *castedSolver ); + } + else + { + GEOS_WARNING( GEOS_FMT( "{} {}: No computing option enabled, no output is scheduled.", + catalogName(), getDataContext() ) ); + } + + m_aggregator->enableRegionStatisticsAggregation(); + + m_aggregator->forRegionStatistics( [&] ( MeshLevel & mesh, RegionStatistics & ) + { + prepareLogTableLayouts( mesh.getName() ); + prepareCsvTableLayouts( mesh.getName() ); + } ); +} + +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_logFormatters.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; + + string_view massUnit = units::getSymbol( m_solver->getMassUnit() ); + + TableLayout tableLayout( { + TableLayout::Column( GEOS_FMT( "Time [{}]", units::getSymbol( units::Unit::Time ))), + TableLayout::Column( "Region" ), // TODO : mention this change in PR description + TableLayout::Column( GEOS_FMT( "Min pressure [{}]", units::getSymbol( units::Unit::Pressure ))), + TableLayout::Column( GEOS_FMT( "Average pressure [{}]", units::getSymbol( units::Unit::Pressure )) ), + TableLayout::Column( GEOS_FMT( "Max pressure [{}]", units::getSymbol( units::Unit::Pressure ) ) ), + TableLayout::Column( GEOS_FMT( "Min delta pressure [{}]", units::getSymbol( units::Unit::Pressure ))), + TableLayout::Column( GEOS_FMT( "Max delta pressure [{}]", units::getSymbol( units::Unit::Pressure ))), + TableLayout::Column( GEOS_FMT( "Min temperature [{}]", units::getSymbol( units::Unit::Temperature ) )), + TableLayout::Column( GEOS_FMT( "Average temperature [{}]", units::getSymbol( units::Unit::Temperature ) )), + TableLayout::Column( GEOS_FMT( "Max temperature [{}]", units::getSymbol( units::Unit::Temperature ) )), + TableLayout::Column( GEOS_FMT( "Total dynamic pore volume [{}]", units::getSymbol( units::Unit::ReservoirVolume ) )), + TableLayout::Column( GEOS_FMT( "Total fluid mass [{}]", massUnit )), + } ); + + auto & csvFormatter = m_csvFormatters.get_inserted( string( meshName ) ); + csvFormatter = std::move( std::make_unique< TableCSVFormatter >( tableLayout ) ); + + // output CSV header + std::ofstream outputFile( getCsvFileName( meshName ) ); + outputFile << csvFormatter->headerToString(); + GEOS_LOG( GEOS_FMT( "table {} : {}", meshName, csvFormatter->headerToString() ) ); // TODO : remove this log +} + +string StatsTask::getCsvFileName( string_view meshName ) const +{ return GEOS_FMT( "{}/{}.csv", m_outputDir, meshName ); } + +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 & ) +{ + // current time is time_n + dt. TODO: verify implication of events ordering in 'time_n+dt' validity + real64 statsTime = time_n + dt; + + GEOS_ERROR_IF( !m_aggregator, + "No statistics aggregator initialized!", getDataContext() ); + + m_aggregator->computeRegionsStatistics( statsTime ); + + m_aggregator->forRegionStatistics( [&] ( MeshLevel & mesh, RegionStatistics & meshRegionsStatistics ) + { + outputLogStats( statsTime, mesh, meshRegionsStatistics ); + outputCsvStats( statsTime, mesh, meshRegionsStatistics ); + } ); + + return false; +} + +void StatsTask::outputLogStats( real64 const statsTime, + MeshLevel & mesh, + RegionStatistics & meshRegionsStatistics ) +{ + if( MpiWrapper::commRank() > 0 || !isLogLevelActive< logInfo::Statistics >( this->getLogLevel() ) ) + return; + + auto const formatterIter = m_logFormatters.find( mesh.getName() ); + if( formatterIter==m_logFormatters.end()) + return; + + TableTextFormatter const & formatter = *formatterIter->second; + TableData tableData; + static constexpr auto merge = CellType::MergeNext; + + 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 = units::getSymbol( units::ReservoirVolume ); + + tableData.getErrorsList().appendErrors( m_aggregator->getWarnings() ); + + tableData.addRow( "Statistics time", merge, merge, statsTime ); + + // lamda to apply for each region statistics + auto const outputRegionStats = [&] ( string_view targetName, RegionStatistics & stats ) + { + tableData.addSeparator(); + tableData.addRow( merge, merge, merge, "" ); + + tableData.addSeparator(); + tableData.addRow( merge, merge, merge, targetName ); + tableData.addSeparator(); + tableData.addRow( "statistics", "min", "average", "max" ); + tableData.addSeparator(); + + tableData.addRow( GEOS_FMT( "Pressure [{}]", pressureUnit ), + stats.m_minPressure, stats.m_averagePressure, stats.m_maxPressure ); + tableData.addRow( GEOS_FMT( "Delta pressure [{}]", pressureUnit ), + stats.m_minDeltaPressure, "/", stats.m_maxDeltaPressure ); + tableData.addRow( GEOS_FMT( "Temperature [{}]", tempUnit ), + stats.m_minTemperature, stats.m_averageTemperature, stats.m_maxTemperature ); + + tableData.addSeparator(); + + tableData.addRow( GEOS_FMT( "Total dynamic pore volume [{}]", resVolUnit ), + "all", CellType::MergeNext, stats.m_totalDynamicPoreVolume ); + + tableData.addRow( GEOS_FMT( "Total fluid mass [{}]", massUnit ), + "all", CellType::MergeNext, stats.m_totalMass ); + }; + + // apply the output lambda for the mesh then each regions + outputRegionStats( GEOS_FMT( "Discretization '{}'", mesh.getName() ), meshRegionsStatistics ); + + m_aggregator->forRegionStatistics( mesh, meshRegionsStatistics, + [&] ( CellElementRegion & region, RegionStatistics & stats ) + { + outputRegionStats( GEOS_FMT( "Region '{}'", region.getName() ), stats ); + } ); + + // output to log + GEOS_LOG_RANK_0( formatter.toString( tableData ) ); +} + +void StatsTask::outputCsvStats( real64 statsTime, + MeshLevel & mesh, + RegionStatistics & meshRegionsStatistics ) +{ + if( MpiWrapper::commRank() > 0 || m_writeCSV == 0 ) + return; + + auto const formatterIter = m_csvFormatters.find( mesh.getName() ); + if( formatterIter==m_csvFormatters.end()) + return; + + TableCSVFormatter const & formatter = *formatterIter->second; + TableData tableData; + + stdVector< string > row; + row.reserve( formatter.getLayout().getTotalLowermostColumnCount() ); + + // lamda to apply for each region statistics + auto const outputRegionStats = [&] ( string_view targetName, RegionStatistics & stats ) + { + row.clear(); + row.insert( row.begin(), + { std::to_string( statsTime ), + string( targetName ), + std::to_string( stats.m_minPressure ), + std::to_string( stats.m_averagePressure ), + std::to_string( stats.m_maxPressure ), + std::to_string( stats.m_minDeltaPressure ), + std::to_string( stats.m_maxDeltaPressure ), + std::to_string( stats.m_minTemperature ), + std::to_string( stats.m_averageTemperature ), + std::to_string( stats.m_maxTemperature ), + std::to_string( stats.m_totalDynamicPoreVolume ), + std::to_string( stats.m_totalMass ), + } ); + + tableData.addRow( row ); + }; + + // apply the output lambda for the mesh then each regions + outputRegionStats( mesh.getName(), meshRegionsStatistics ); + + m_aggregator->forRegionStatistics( mesh, meshRegionsStatistics, + [&] ( CellElementRegion & region, RegionStatistics & stats ) + { + outputRegionStats( region.getName(), stats ); + } ); + + // append to csv file + std::ofstream outputFile( getCsvFileName( mesh.getName() ), std::ios_base::app ); + outputFile << formatter.dataToString( tableData ); + outputFile.close(); +} + +REGISTER_CATALOG_ENTRY( TaskBase, + StatsTask, + string const &, dataRepository::Group * const ) + +} /* namespace singlePhaseStatistics */ + +} /* namespace geos */ diff --git a/src/coreComponents/physicsSolvers/fluidFlow/SinglePhaseStatistics.hpp b/src/coreComponents/physicsSolvers/fluidFlow/SinglePhaseStatisticsTask.hpp similarity index 54% rename from src/coreComponents/physicsSolvers/fluidFlow/SinglePhaseStatistics.hpp rename to src/coreComponents/physicsSolvers/fluidFlow/SinglePhaseStatisticsTask.hpp index 4906790b026..2d83a6b5a3f 100644 --- a/src/coreComponents/physicsSolvers/fluidFlow/SinglePhaseStatistics.hpp +++ b/src/coreComponents/physicsSolvers/fluidFlow/SinglePhaseStatisticsTask.hpp @@ -21,18 +21,22 @@ #define SRC_CORECOMPONENTS_PHYSICSSOLVERS_FLUIDFLOW_SINGLEPHASESTATISTICS_HPP_ #include "physicsSolvers/FieldStatisticsBase.hpp" +#include "physicsSolvers/fluidFlow/SinglePhaseStatisticsAggregator.hpp" namespace geos { class SinglePhaseBase; +namespace singlePhaseStatistics +{ + /** - * @class SinglePhaseStatistics + * @class singlePhaseStatistics::Task * - * Task class allowing for the computation of aggregate statistics in single-phase simulations + * Task class allowing for the computation of aggregate statistics in single phase simulations */ -class SinglePhaseStatistics : public FieldStatisticsBase< SinglePhaseBase > +class StatsTask : public FieldStatisticsBase< SinglePhaseBase > { public: @@ -41,14 +45,11 @@ class SinglePhaseStatistics : public FieldStatisticsBase< SinglePhaseBase > * @param[in] name the name of the task coming from the xml * @param[in] parent the parent group of the task */ - SinglePhaseStatistics( const string & name, - Group * const parent ); + StatsTask( const string & name, dataRepository::Group * const parent ); /// Accessor for the catalog name static string catalogName() { return "SinglePhaseStatistics"; } - /// Accessor for the region statistics catalog name - static string regionStatisticsName() { return "regionStatistics"; } /** * @defgroup Tasks Interface Functions * @@ -65,64 +66,46 @@ class SinglePhaseStatistics : public FieldStatisticsBase< SinglePhaseBase > /**@}*/ + StatsAggregator & getStatisticsAggregator() + { return *m_aggregator; } - /** - * @struct viewKeyStruct holds char strings and viewKeys for fast lookup - */ - struct viewKeyStruct - { - /// String for the region statistics - constexpr static char const * regionStatisticsString() { return "regionStatistics"; } - }; - - 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; - - // fluid mass - real64 totalMass; - - /// 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; - }; + StatsAggregator const & getStatisticsAggregator() const + { return *m_aggregator; } private: using Base = FieldStatisticsBase< SinglePhaseBase >; - /** - * @brief Compute some statistics on the reservoir (average field pressure, etc) - * @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 postInputInitialization() override; void registerDataOnMesh( Group & meshBodies ) override; + void prepareLogTableLayouts( string_view tableName ); + + void prepareCsvTableLayouts( string_view tableName ); + + string getCsvFileName( string_view meshName ) const; + + void outputLogStats( real64 statsTime, + MeshLevel & mesh, + RegionStatistics & meshRegionsStatistics ); + + void outputCsvStats( real64 statsTime, + MeshLevel & mesh, + RegionStatistics & meshRegionsStatistics ); + + /// For each discretization (MeshLevel name), table formatter for log output. + stdMap< string, std::unique_ptr< TableTextFormatter > > m_logFormatters; + + /// For each discretization (MeshLevel name), table formatter for csv output. + stdMap< string, std::unique_ptr< TableCSVFormatter > > m_csvFormatters; + + // mesh statistics aggregator + std::unique_ptr< StatsAggregator > m_aggregator; + }; +} /* namespace singlePhaseStatistics */ } /* namespace geos */ diff --git a/src/coreComponents/physicsSolvers/fluidFlow/unitTests/testFlowStatistics.cpp b/src/coreComponents/physicsSolvers/fluidFlow/unitTests/testFlowStatistics.cpp index 356661e434e..389cb79770b 100644 --- a/src/coreComponents/physicsSolvers/fluidFlow/unitTests/testFlowStatistics.cpp +++ b/src/coreComponents/physicsSolvers/fluidFlow/unitTests/testFlowStatistics.cpp @@ -17,8 +17,13 @@ #include "integrationTests/testingUtilities/TestingTasks.hpp" #include "mainInterface/initialization.hpp" #include "mainInterface/GeosxState.hpp" +#include "physicsSolvers/fluidFlow/CompositionalMultiphaseStatisticsAggregator.hpp" +#include "physicsSolvers/fluidFlow/CompositionalMultiphaseStatisticsTask.hpp" +#include "physicsSolvers/fluidFlow/SinglePhaseBase.hpp" +#include "physicsSolvers/fluidFlow/CompositionalMultiphaseBase.hpp" +#include "physicsSolvers/fluidFlow/SinglePhaseStatisticsAggregator.hpp" #include "physicsSolvers/fluidFlow/SourceFluxStatistics.hpp" -#include "physicsSolvers/fluidFlow/SinglePhaseStatistics.hpp" +#include "physicsSolvers/fluidFlow/SinglePhaseStatisticsTask.hpp" #include @@ -46,6 +51,7 @@ struct TestInputs string timeStepFluxStatsPath; string wholeSimFluxStatsPath; string flowSolverPath; + string statsTaskPath; // rates for each timesteps, for each phases array2d< real64 > sourceRates; @@ -206,27 +212,42 @@ void setRateTable( array2d< real64 > & rateTable, std::initializer_list< std::in } } -real64 getTotalFluidMass( ProblemManager & problem, string_view flowSolverPath ) +template< typename SolverType > +real64 getTotalFluidMass( ProblemManager & problem, SolverType & solver, string_view statsTaskPath ); + +template<> +real64 getTotalFluidMass< SinglePhaseBase >( ProblemManager & problem, + SinglePhaseBase & solver, + string_view statsTaskPath ) { - real64 totalMass = 0.0; - PhysicsSolverBase const & solver = problem.getGroupByPath< PhysicsSolverBase >( string( flowSolverPath ) ); - solver.forDiscretizationOnMeshTargets( problem.getDomainPartition().getMeshBodies(), - [&] ( string const &, - MeshLevel & mesh, - string_array const & ) - { - mesh.getElemManager().forElementRegions( [&]( ElementRegionBase & region ) - { - SinglePhaseStatistics::RegionStatistics & regionStats = region.getReference< SinglePhaseStatistics::RegionStatistics >( - SinglePhaseStatistics::viewKeyStruct::regionStatisticsString() ); + using namespace singlePhaseStatistics; + MeshLevel & mesh = problem.getDomainPartition() + .getMeshBody( 0 ) + .getMeshLevel( solver.getDiscretizationName() ); + StatsTask const & statsTask = problem.getGroupByPath< StatsTask >( string( statsTaskPath ) ); + StatsAggregator const & statsAggregator = statsTask.getStatisticsAggregator(); + RegionStatistics const & stats = statsAggregator.getRegionsStatistics( mesh ); + return stats.m_totalMass; +} - totalMass += regionStats.totalMass; - } ); - } ); +template<> +real64 getTotalFluidMass< CompositionalMultiphaseBase >( ProblemManager & problem, + CompositionalMultiphaseBase & solver, + string_view statsTaskPath ) +{ + using namespace compositionalMultiphaseStatistics; + MeshLevel & mesh = problem.getDomainPartition() + .getMeshBody( 0 ) + .getMeshLevel( solver.getDiscretizationName() ); + StatsTask const & statsTask = problem.getGroupByPath< StatsTask >( string( statsTaskPath ) ); + StatsAggregator const & statsAggregator = statsTask.getStatisticsAggregator(); + RegionStatistics const & stats = statsAggregator.getRegionsStatistics( mesh ); + double totalMass = 0.0; + for( integer i = 0; i < statsAggregator.getNumPhases(); ++i ) + totalMass += stats.m_phaseMass[i]; return totalMass; } - /** * @brief Verification that the source flux statistics are correct for the current timestep * @param expectedMasses the expected mass values per phase @@ -529,6 +550,7 @@ TestSet getTestSet() testInputs.timeStepFluxStatsPath = "/Tasks/timeStepFluxStats"; testInputs.wholeSimFluxStatsPath = "/Tasks/wholeSimFluxStats"; testInputs.flowSolverPath = "/Solvers/testSolver"; + testInputs.statsTaskPath = "/Tasks/timeStepReservoirStats"; testInputs.dt = 500.0; testInputs.sourceElementsCount = 2; @@ -565,8 +587,10 @@ TEST_F( FlowStatisticsTest, checkSinglePhaseFluxStatistics ) setupProblemFromXML( problem, testSet.inputs.xmlInput.data() ); - real64 firstMass; + SinglePhaseBase & flowSolver = + problem.getGroupByPath< SinglePhaseBase >( testSet.inputs.flowSolverPath ); + real64 firstMass; TimeStepChecker & timeStepChecker = problem.getGroupByPath< TimeStepChecker >( testSet.inputs.timeStepCheckerPath ); timeStepChecker.setTimeStepCheckingFunction( [&]( real64 const time_n ) { @@ -578,7 +602,7 @@ TEST_F( FlowStatisticsTest, checkSinglePhaseFluxStatistics ) if( !passedFirstTimeStep ) { passedFirstTimeStep = true; - firstMass = getTotalFluidMass( problem, testSet.inputs.flowSolverPath ); + firstMass = getTotalFluidMass( problem, flowSolver, testSet.inputs.statsTaskPath ); } } ); @@ -589,12 +613,12 @@ TEST_F( FlowStatisticsTest, checkSinglePhaseFluxStatistics ) checkWholeSimTimeStepStats( problem, testSet, timeStepChecker ); // check singlephasestatistics results - real64 const lastMass = getTotalFluidMass( problem, testSet.inputs.flowSolverPath ); + real64 const lastMass = getTotalFluidMass( problem, flowSolver, testSet.inputs.statsTaskPath ); real64 const massDiffTol = 1e-7; EXPECT_NEAR( lastMass - firstMass, -testSet.totalMassProd[0], massDiffTol * std::abs( testSet.totalMassProd[0] ) ) << GEOS_FMT( "{} total mass difference from start to end is not consistent with fluxes production.", - SinglePhaseStatistics::catalogName() ); + singlePhaseStatistics::StatsTask::catalogName() ); } @@ -806,6 +830,7 @@ TestSet getTestSet() testInputs.timeStepFluxStatsPath = "/Tasks/timeStepFluxStats"; testInputs.wholeSimFluxStatsPath = "/Tasks/wholeSimFluxStats"; testInputs.flowSolverPath = "/Solvers/testSolver"; + testInputs.statsTaskPath = "/Tasks/timeStepReservoirStats"; testInputs.dt = 500.0; testInputs.sourceElementsCount = 1; @@ -854,12 +879,23 @@ TEST_F( FlowStatisticsTest, checkMultiPhaseFluxStatisticsMass ) setupProblemFromXML( problem, testSet.inputs.xmlInput.data() ); + CompositionalMultiphaseBase & flowSolver = + problem.getGroupByPath< CompositionalMultiphaseBase >( testSet.inputs.flowSolverPath ); + + real64 firstMass; TimeStepChecker & timeStepChecker = problem.getGroupByPath< TimeStepChecker >( testSet.inputs.timeStepCheckerPath ); timeStepChecker.setTimeStepCheckingFunction( [&]( real64 const time_n ) { integer const timestepId = timeStepChecker.getTestedTimeStepCount(); checkTimeStepStats( testSet, time_n, timestepId ); checkTimeStepFluxStats( problem, testSet, time_n, timestepId ); + + static bool passedFirstTimeStep = false; + if( !passedFirstTimeStep ) + { + passedFirstTimeStep = true; + firstMass = getTotalFluidMass( problem, flowSolver, testSet.inputs.statsTaskPath ); + } } ); // run simulation @@ -867,6 +903,14 @@ TEST_F( FlowStatisticsTest, checkMultiPhaseFluxStatisticsMass ) checkWholeSimFluxStats( problem, testSet ); checkWholeSimTimeStepStats( problem, testSet, timeStepChecker ); + + // check compositionalmultiphasestatistics results + real64 const lastMass = getTotalFluidMass( problem, flowSolver, testSet.inputs.statsTaskPath ); + real64 const massDiffTol = 1e-5; + EXPECT_NEAR( lastMass - firstMass, + -( testSet.totalMassProd[0] + testSet.totalMassProd[1] ), + massDiffTol * std::abs( testSet.totalMassProd[0] + testSet.totalMassProd[1] ) ) << GEOS_FMT( "{} total mass difference from start to end is not consistent with fluxes production.", + singlePhaseStatistics::StatsTask::catalogName() ); } @@ -1078,6 +1122,7 @@ TestSet getTestSet() testInputs.timeStepFluxStatsPath = "/Tasks/timeStepFluxStats"; testInputs.wholeSimFluxStatsPath = "/Tasks/wholeSimFluxStats"; testInputs.flowSolverPath = "/Solvers/testSolver"; + testInputs.statsTaskPath = "/Tasks/timeStepReservoirStats"; testInputs.dt = 500.0; testInputs.sourceElementsCount = 1; @@ -1130,12 +1175,23 @@ TEST_F( FlowStatisticsTest, checkMultiPhaseFluxStatisticsMol ) setupProblemFromXML( problem, testSet.inputs.xmlInput.data() ); + CompositionalMultiphaseBase & flowSolver = + problem.getGroupByPath< CompositionalMultiphaseBase >( testSet.inputs.flowSolverPath ); + + real64 firstMass; TimeStepChecker & timeStepChecker = problem.getGroupByPath< TimeStepChecker >( testSet.inputs.timeStepCheckerPath ); timeStepChecker.setTimeStepCheckingFunction( [&]( real64 const time_n ) { integer const timestepId = timeStepChecker.getTestedTimeStepCount(); checkTimeStepStats( testSet, time_n, timestepId ); checkTimeStepFluxStats( problem, testSet, time_n, timestepId ); + + static bool passedFirstTimeStep = false; + if( !passedFirstTimeStep ) + { + passedFirstTimeStep = true; + firstMass = getTotalFluidMass( problem, flowSolver, testSet.inputs.statsTaskPath ); + } } ); // run simulation @@ -1143,6 +1199,14 @@ TEST_F( FlowStatisticsTest, checkMultiPhaseFluxStatisticsMol ) checkWholeSimFluxStats( problem, testSet ); checkWholeSimTimeStepStats( problem, testSet, timeStepChecker ); + + // check compositionalmultiphasestatistics results + real64 const lastMass = getTotalFluidMass( problem, flowSolver, testSet.inputs.statsTaskPath ); + real64 const massDiffTol = 1e-5; + EXPECT_NEAR( lastMass - firstMass, + -( testSet.totalMassProd[0] + testSet.totalMassProd[1] ), + massDiffTol * std::abs( testSet.totalMassProd[0] + testSet.totalMassProd[1] ) ) << GEOS_FMT( "{} total mass difference from start to end is not consistent with fluxes production.", + singlePhaseStatistics::StatsTask::catalogName() ); } diff --git a/src/coreComponents/physicsSolvers/fluidFlow/wells/CompositionalMultiphaseWell.cpp b/src/coreComponents/physicsSolvers/fluidFlow/wells/CompositionalMultiphaseWell.cpp index 275b0451d82..e687527f808 100644 --- a/src/coreComponents/physicsSolvers/fluidFlow/wells/CompositionalMultiphaseWell.cpp +++ b/src/coreComponents/physicsSolvers/fluidFlow/wells/CompositionalMultiphaseWell.cpp @@ -19,19 +19,23 @@ #include "CompositionalMultiphaseWell.hpp" +#include "LvArray/src/system.hpp" #include "codingUtilities/Utilities.hpp" #include "common/DataTypes.hpp" #include "common/TimingMacros.hpp" +#include "common/logger/Logger.hpp" #include "constitutive/ConstitutiveManager.hpp" #include "constitutive/fluid/multifluid/MultiFluidBase.hpp" #include "constitutive/fluid/multifluid/MultiFluidFields.hpp" #include "constitutive/fluid/multifluid/MultiFluidSelector.hpp" #include "dataRepository/Group.hpp" #include "mesh/DomainPartition.hpp" +#include "mesh/MeshBody.hpp" #include "mesh/PerforationFields.hpp" #include "mesh/WellElementSubRegion.hpp" #include "mesh/mpiCommunications/CommunicationTools.hpp" #include "physicsSolvers/LogLevelsInfo.hpp" +#include "physicsSolvers/fluidFlow/CompositionalMultiphaseBase.hpp" #include "physicsSolvers/fluidFlow/wells/LogLevelsInfo.hpp" #include "physicsSolvers/fluidFlow/wells/WellSolverBaseFields.hpp" #include "physicsSolvers/fluidFlow/wells/WellFields.hpp" @@ -48,7 +52,7 @@ #include "physicsSolvers/fluidFlow/kernels/compositional/PhaseVolumeFractionKernel.hpp" #include "physicsSolvers/fluidFlow/kernels/compositional/ThermalPhaseVolumeFractionKernel.hpp" #include "physicsSolvers/fluidFlow/kernels/compositional/FluidUpdateKernel.hpp" -#include "physicsSolvers/fluidFlow/CompositionalMultiphaseStatistics.hpp" +#include "physicsSolvers/fluidFlow/CompositionalMultiphaseStatisticsAggregator.hpp" #if defined( __INTEL_COMPILER ) #pragma GCC optimize "O0" @@ -60,6 +64,19 @@ namespace geos using namespace dataRepository; using namespace constitutive; using namespace fields; +using namespace compositionalMultiphaseStatistics; + +CompositionalMultiphaseBase & getFlowSolver( CompositionalMultiphaseWell & wellSolver ) +{ + // TODO: change the way we access the flowSolver here + return wellSolver.getParent().getGroup< CompositionalMultiphaseBase >( wellSolver.getFlowSolverName() ); +} + +CompositionalMultiphaseBase const & getFlowSolver( CompositionalMultiphaseWell const & wellSolver ) +{ + // TODO: change the way we access the flowSolver here + return wellSolver.getParent().getGroup< CompositionalMultiphaseBase >( wellSolver.getFlowSolverName() ); +} CompositionalMultiphaseWell::CompositionalMultiphaseWell( const string & name, Group * const parent ) @@ -366,7 +383,7 @@ void CompositionalMultiphaseWell::validateConstitutiveModels( DomainPartition co GEOS_MARK_FUNCTION; ConstitutiveManager const & cm = domain.getConstitutiveManager(); - CompositionalMultiphaseBase const & flowSolver = getParent().getGroup< CompositionalMultiphaseBase >( getFlowSolverName() ); + CompositionalMultiphaseBase const & flowSolver = getFlowSolver( *this ); string const referenceFluidName = flowSolver.referenceFluidModelName(); MultiFluidBase const & referenceFluid = cm.getConstitutiveRelation< MultiFluidBase >( m_referenceFluidModelName ); @@ -429,33 +446,21 @@ void CompositionalMultiphaseWell::validateWellConstraints( real64 const & time_n WellElementSubRegion const & subRegion ) { WellControls & wellControls = getWellControls( subRegion ); + if( !wellControls.useSurfaceConditions() ) { - bool const useSeg =wellControls.referenceReservoirRegion().empty(); - GEOS_WARNING_IF( useSeg, - GEOS_FMT( "WellControls {} not set and well constraint fluid property calculations will use " - "top segement pressure and temp ", - WellControls::viewKeyStruct::referenceReservoirRegionString() ) ); - if( useSeg ) - { - wellControls.setRegionAveragePressure( -1 ); - wellControls.setRegionAverageTemperature( -1 ); - } - else - { - // Check if region name exists in list of Reservoir's target regions - string const regionName = wellControls.referenceReservoirRegion(); - CompositionalMultiphaseBase const & flowSolver = getParent().getGroup< CompositionalMultiphaseBase >( getFlowSolverName() ); - string_array const & targetRegionsNames = flowSolver.getTargetRegionNames(); - auto const pos = std::find( targetRegionsNames.begin(), targetRegionsNames.end(), regionName ); - GEOS_ERROR_IF( pos == targetRegionsNames.end(), - GEOS_FMT( "Region {} is not a target of the reservoir solver and cannot be used for referenceReservoirRegion in WellControl {}.", - regionName, wellControls.getName() ), - getDataContext() ); - + bool const useSegmentValues = wellControls.referenceReservoirRegion().empty(); + static bool firstNoRefRegionMsg = true; + if( useSegmentValues && firstNoRefRegionMsg ) + { + GEOS_WARNING( WellControls::viewKeyStruct::referenceReservoirRegionString() << + " not set: well constraint fluid property calculations will use top segement pressure and temp ", + wellControls.getDataContext() ); + firstNoRefRegionMsg = false; } } + string const & fluidName = subRegion.getReference< string >( viewKeyStruct::fluidNamesString()); MultiFluidBase const & fluid = subRegion.getConstitutiveModel< MultiFluidBase >( fluidName ); @@ -662,7 +667,7 @@ void CompositionalMultiphaseWell::updateBHPForConstraint( WellElementSubRegion & } -void CompositionalMultiphaseWell::updateVolRatesForConstraint( ElementRegionManager const & elemManager, WellElementSubRegion const & subRegion ) +void CompositionalMultiphaseWell::updateVolRatesForConstraint( WellElementSubRegion const & subRegion ) { GEOS_MARK_FUNCTION; @@ -680,9 +685,6 @@ void CompositionalMultiphaseWell::updateVolRatesForConstraint( ElementRegionMana WellControls & wellControls = getWellControls( subRegion ); // subRegion data - - arrayView1d< real64 const > const & pres = subRegion.getField< well::pressure >(); - arrayView1d< real64 const > const & temp = subRegion.getField< well::temperature >(); arrayView1d< real64 const > const & connRate = subRegion.getField< well::mixtureConnectionRate >(); arrayView2d< real64 const, compflow::USD_COMP > const & compFrac = subRegion.getField< well::globalCompFraction >(); @@ -707,47 +709,8 @@ void CompositionalMultiphaseWell::updateVolRatesForConstraint( ElementRegionMana string const massUnit = m_useMass ? "kg" : "mol"; integer const useSurfaceConditions = wellControls.useSurfaceConditions(); - real64 flashPressure; - real64 flashTemperature; - if( useSurfaceConditions ) - { - // use surface conditions - flashPressure = wellControls.getSurfacePressure(); - flashTemperature = wellControls.getSurfaceTemperature(); - } - else - { - if( !wellControls.referenceReservoirRegion().empty() ) - { - ElementRegionBase const & region = elemManager.getRegion( wellControls.referenceReservoirRegion()); - GEOS_ERROR_IF ( !region.hasWrapper( CompositionalMultiphaseStatistics::regionStatisticsName() ), - GEOS_FMT( "WellControl {} referenceReservoirRegion field requires CompositionalMultiphaseStatistics to be configured for region {} ", - wellControls.getName(), wellControls.referenceReservoirRegion() ), - getDataContext() ); - - 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 {} ", - wellControls.getName(), wellControls.referenceReservoirRegion() ), - getDataContext()); - } - // If flashPressure is not set by region the value is defaulted to -1 and indicates to use top segment conditions - flashPressure = wellControls.getRegionAveragePressure(); - if( flashPressure < 0.0 ) - { - // region name not set, use segment conditions - flashPressure = pres[iwelemRef]; - flashTemperature = temp[iwelemRef]; - } - else - { - // use reservoir region averages - flashTemperature = wellControls.getRegionAverageTemperature(); - } - } + ReferenceConditions const refConditions = getReferenceConditions( subRegion ); + arrayView1d< real64 > const & currentPhaseVolRate = wellControls.getReference< array1d< real64 > >( CompositionalMultiphaseWell::viewKeyStruct::currentPhaseVolRateString() ); arrayView2d< real64 > const & dCurrentPhaseVolRate = @@ -777,8 +740,6 @@ void CompositionalMultiphaseWell::updateVolRatesForConstraint( ElementRegionMana forAll< serialPolicy >( 1, [&numComp, &numPhase, fluidSeparatorWrapper, - pres, - temp, compFrac, dCompFrac_dCompDens, connRate, @@ -790,8 +751,7 @@ void CompositionalMultiphaseWell::updateVolRatesForConstraint( ElementRegionMana dPhaseFrac, logSurfaceCondition, &useSurfaceConditions, - &flashPressure, - &flashTemperature, + &refConditions, ¤tTotalVolRate, dCurrentTotalVolRate, currentPhaseVolRate, @@ -811,24 +771,18 @@ void CompositionalMultiphaseWell::updateVolRatesForConstraint( ElementRegionMana // - Surface conditions: using the surface pressure provided by the user // - Segment conditions: using the pressure in the top element // - Reservoir conditions: using the average region pressure - if( useSurfaceConditions ) + fluidSeparatorWrapper.update( iwelemRef, 0, + refConditions.pressure, refConditions.temperature, + compFrac[iwelemRef] ); + if( useSurfaceConditions && logSurfaceCondition ) { - // we need to compute the surface density - fluidSeparatorWrapper.update( iwelemRef, 0, flashPressure, flashTemperature, compFrac[iwelemRef] ); - if( logSurfaceCondition ) - { - GEOS_LOG_RANK( GEOS_FMT( "{}: surface density computed with P_surface = {} Pa and T_surface = {} K", - wellControlsName, flashPressure, flashTemperature ) ); - } + GEOS_LOG_RANK( GEOS_FMT( "{}: surface density computed with P_surface = {} Pa and T_surface = {} K", + wellControlsName, refConditions.pressure, refConditions.temperature ) ); + } #ifdef GEOS_USE_HIP - GEOS_UNUSED_VAR( wellControlsName ); + GEOS_UNUSED_VAR( wellControlsName ); #endif - } - else - { - fluidSeparatorWrapper.update( iwelemRef, 0, flashPressure, flashTemperature, compFrac[iwelemRef] ); - } // Step 2: update the total volume rate real64 const currentTotalRate = connRate[iwelemRef]; @@ -916,7 +870,90 @@ void CompositionalMultiphaseWell::updateVolRatesForConstraint( ElementRegionMana } ); } +void CompositionalMultiphaseWell::precomputeReferenceConditions( real64 const time_n, + Group & meshBodies, + MeshBody & meshBody, + WellElementSubRegion const & subRegion ) +{ + WellControls & wellControls = getWellControls( subRegion ); + integer const useSurfaceConditions = wellControls.useSurfaceConditions(); + if( !useSurfaceConditions ) + { + string_view refRegionName = wellControls.referenceReservoirRegion(); + bool const useSegmentValues = refRegionName.empty(); + if( useSegmentValues ) + { + wellControls.setRegionAveragePressure( -1 ); + wellControls.setRegionAverageTemperature( -1 ); + } + else + { + auto & flowSolver = getParent().getGroup< CompositionalMultiphaseBase >( getFlowSolverName() ); + MeshLevel & flowMeshLevel = meshBody.getMeshLevel( flowSolver.getDiscretizationName() ); + + if( !m_reservoirStatsAggregator ) + { // lazily initialize the region statistics aggregator + m_reservoirStatsAggregator = std::make_unique< StatsAggregator >( wellControls.getDataContext() ); + m_reservoirStatsAggregator->initStatisticsAggregation( meshBodies, flowSolver ); + m_reservoirStatsAggregator->enableRegionStatisticsAggregation(); + } + + RegionStatistics & stats = m_reservoirStatsAggregator->getRegionStatistics( flowMeshLevel, refRegionName ); + + // compute region stats only if needed (could have already been done for another subRegion) + if( !m_reservoirStatsAggregator->isComputed( time_n, stats ) ) + m_reservoirStatsAggregator->computeRegionsStatistics( time_n ); + + GEOS_WARNING_IF( stats.m_averagePressure <= 0.0, + GEOS_FMT( "No region average quantities computed in reference region '{}'.", + wellControls.referenceReservoirRegion() ), + wellControls.getWrapperDataContext( WellControls::viewKeyStruct::referenceReservoirRegionString() ), + getDataContext() ); + + wellControls.setRegionAveragePressure( stats.m_averagePressure ); + wellControls.setRegionAverageTemperature( stats.m_averageTemperature ); + } + } +} + +CompositionalMultiphaseWell::ReferenceConditions +CompositionalMultiphaseWell::getReferenceConditions( WellElementSubRegion const & subRegion ) +{ + WellControls & wellControls = getWellControls( subRegion ); + integer const useSurfaceConditions = wellControls.useSurfaceConditions(); + if( useSurfaceConditions ) + { + // use surface conditions + return { + /* .pressure = */ wellControls.getSurfacePressure(), + /* .temperature = */ wellControls.getSurfaceTemperature(), + }; + } + else + { + if( wellControls.getRegionAveragePressure() > 0.0 && wellControls.getRegionAverageTemperature() > 0.0 ) + { // reference region condition properly computed, we can return them + return { + /* .pressure = */ wellControls.getRegionAveragePressure(), + /* .temperature = */ wellControls.getRegionAverageTemperature(), + }; + } + else + { // region average stats not initialized or initialized, fallback to top segment values + GEOS_WARNING( "CompositionalMultiphaseWell: region average statsistics of reference region not initialized," + " fallback to top segment values.", + wellControls.getDataContext() ); + arrayView1d< real64 const > const & pres = subRegion.getField< well::pressure >(); + arrayView1d< real64 const > const & temp = subRegion.getField< well::temperature >(); + localIndex const iwelemRef = subRegion.getTopWellElementIndex(); + return { + /* .pressure = */ pres[iwelemRef], + /* .temperature = */ temp[iwelemRef], + }; + } + } +} void CompositionalMultiphaseWell::updateFluidModel( WellElementSubRegion & subRegion ) { @@ -1008,7 +1045,7 @@ void CompositionalMultiphaseWell::updateState( DomainPartition & domain ) WellControls & wellControls = getWellControls( subRegion ); if( wellControls.getWellStatus() == WellControls::Status::OPEN ) { - real64 const maxRegionPhaseVolFrac = updateSubRegionState( elemManager, subRegion ); + real64 const maxRegionPhaseVolFrac = updateSubRegionState( subRegion ); maxPhaseVolFrac = LvArray::math::max( maxRegionPhaseVolFrac, maxPhaseVolFrac ); } } ); @@ -1021,14 +1058,14 @@ void CompositionalMultiphaseWell::updateState( DomainPartition & domain ) } -real64 CompositionalMultiphaseWell::updateSubRegionState( ElementRegionManager const & elemManager, WellElementSubRegion & subRegion ) +real64 CompositionalMultiphaseWell::updateSubRegionState( WellElementSubRegion & subRegion ) { // update properties updateGlobalComponentFraction( subRegion ); // update volumetric rates for the well constraints // note: this must be called before updateFluidModel - updateVolRatesForConstraint( elemManager, subRegion ); + updateVolRatesForConstraint( subRegion ); // update densities, phase fractions, phase volume fractions @@ -1047,16 +1084,17 @@ void CompositionalMultiphaseWell::initializeWells( DomainPartition & domain, rea integer const numComp = m_numComponents; integer const numPhase = m_numPhases; - // TODO: change the way we access the flowSolver here - CompositionalMultiphaseBase const & flowSolver = getParent().getGroup< CompositionalMultiphaseBase >( getFlowSolverName() ); + Group & meshBodies = domain.getMeshBodies(); + CompositionalMultiphaseBase const & flowSolver = getFlowSolver( *this ); // loop over the wells - forDiscretizationOnMeshTargets( domain.getMeshBodies(), [&] ( string const &, - MeshLevel & mesh, - string_array const & regionNames ) + forDiscretizationOnMeshTargets( meshBodies, [&] ( string const & meshBodyName, + MeshLevel & mesh, + string_array const & regionNames ) { - + MeshBody & meshBody = domain.getMeshBody( meshBodyName ); ElementRegionManager & elemManager = mesh.getElemManager(); + compositionalMultiphaseWellKernels::PresTempCompFracInitializationKernel::CompFlowAccessors resCompFlowAccessors( mesh.getElemManager(), flowSolver.getName() ); compositionalMultiphaseWellKernels::PresTempCompFracInitializationKernel::MultiFluidAccessors @@ -1146,7 +1184,8 @@ void CompositionalMultiphaseWell::initializeWells( DomainPartition & domain, rea wellElemCompDens ); // 5) Recompute the pressure-dependent properties - updateSubRegionState( elemManager, subRegion ); + precomputeReferenceConditions( time_n, meshBodies, meshBody, subRegion ); + updateSubRegionState( subRegion ); // 6) Estimate the well rates // TODO: initialize rates using perforation rates @@ -1726,7 +1765,7 @@ void CompositionalMultiphaseWell::computePerforationRates( real64 const & time_n { // TODO: change the way we access the flowSolver here - CompositionalMultiphaseBase const & flowSolver = getParent().getGroup< CompositionalMultiphaseBase >( getFlowSolverName() ); + CompositionalMultiphaseBase const & flowSolver = getFlowSolver( *this ); ElementRegionManager & elemManager = mesh.getElemManager(); elemManager.forElementSubRegions< WellElementSubRegion >( regionNames, [&]( localIndex const, @@ -1952,7 +1991,7 @@ void CompositionalMultiphaseWell::resetStateToBeginningOfStep( DomainPartition & if( wellControls.isWellOpen( ) ) { - updateSubRegionState( elemManager, subRegion ); + updateSubRegionState( subRegion ); } } ); } ); @@ -2090,11 +2129,14 @@ void CompositionalMultiphaseWell::implicitStepSetup( real64 const & time_n, { WellSolverBase::implicitStepSetup( time_n, dt, domain ); - forDiscretizationOnMeshTargets ( domain.getMeshBodies(), [&] ( string const &, - MeshLevel & mesh, - string_array const & regionNames ) - { + if( m_reservoirStatsAggregator ) + m_reservoirStatsAggregator->setDirty(); // TODO: is it useful? are timestep cut properly managed for this call? + Group & meshBodies = domain.getMeshBodies(); + forDiscretizationOnMeshTargets ( meshBodies, [&] ( string const &, + MeshLevel & mesh, + string_array const & regionNames ) + { ElementRegionManager & elemManager = mesh.getElemManager(); elemManager.forElementSubRegions< WellElementSubRegion >( regionNames, @@ -2146,7 +2188,7 @@ void CompositionalMultiphaseWell::implicitStepSetup( real64 const & time_n, validateWellConstraints( time_n, dt, subRegion ); - updateSubRegionState( elemManager, subRegion ); + updateSubRegionState( subRegion ); } } ) ; diff --git a/src/coreComponents/physicsSolvers/fluidFlow/wells/CompositionalMultiphaseWell.hpp b/src/coreComponents/physicsSolvers/fluidFlow/wells/CompositionalMultiphaseWell.hpp index cb8ed9f1a6f..677ea19d64c 100644 --- a/src/coreComponents/physicsSolvers/fluidFlow/wells/CompositionalMultiphaseWell.hpp +++ b/src/coreComponents/physicsSolvers/fluidFlow/wells/CompositionalMultiphaseWell.hpp @@ -20,8 +20,10 @@ #ifndef GEOS_PHYSICSSOLVERS_FLUIDFLOW_WELLS_COMPOSITIONALMULTIPHASEWELL_HPP_ #define GEOS_PHYSICSSOLVERS_FLUIDFLOW_WELLS_COMPOSITIONALMULTIPHASEWELL_HPP_ +#include "common/DataTypes.hpp" #include "constitutive/fluid/multifluid/Layouts.hpp" #include "constitutive/relativePermeability/Layouts.hpp" +#include "mesh/MeshLevel.hpp" #include "physicsSolvers/fluidFlow/wells/WellSolverBase.hpp" #include "physicsSolvers/fluidFlow/CompositionalMultiphaseBase.hpp" @@ -34,6 +36,11 @@ class ConstitutiveManager; class MultiFluidBase; } +namespace compositionalMultiphaseStatistics +{ +class StatsAggregator; +} + /** * @class CompositionalMultiphaseWell * @@ -143,7 +150,7 @@ class CompositionalMultiphaseWell : public WellSolverBase * @param elemManager the well region manager containing the well * @param subRegion the well subregion containing all the primary and dependent fields */ - void updateVolRatesForConstraint( ElementRegionManager const & elemManager, WellElementSubRegion const & subRegion ); + void updateVolRatesForConstraint( WellElementSubRegion const & subRegion ); /** * @brief Recompute the current BHP pressure @@ -185,7 +192,7 @@ class CompositionalMultiphaseWell : public WellSolverBase */ virtual void updateState( DomainPartition & domain ) override; - virtual real64 updateSubRegionState( ElementRegionManager const & elemManager, WellElementSubRegion & subRegion ) override; + virtual real64 updateSubRegionState( WellElementSubRegion & subRegion ) override; virtual string wellElementDofName() const override { return viewKeyStruct::dofFieldString(); } @@ -358,6 +365,12 @@ class CompositionalMultiphaseWell : public WellSolverBase private: + struct ReferenceConditions + { + real64 pressure; + real64 temperature; + }; + /** * @brief Initialize all the primary and secondary variables in all the wells * @param domain the domain containing the well manager to access individual wells @@ -366,7 +379,12 @@ class CompositionalMultiphaseWell : public WellSolverBase virtual void setConstitutiveNames( ElementSubRegionBase & subRegion ) const override; + void precomputeReferenceConditions( real64 time_n, + Group & meshBodies, + MeshBody & meshBody, + WellElementSubRegion const & subRegion ); + ReferenceConditions getReferenceConditions( WellElementSubRegion const & subRegion ); /// flag indicating whether mass or molar formulation should be used integer m_useMass; @@ -398,6 +416,8 @@ class CompositionalMultiphaseWell : public WellSolverBase /// index of the target phase, used to impose the phase rate constraint localIndex m_targetPhaseIndex; + /// optional statistics aggregator to get the average pressure of simulated region + std::unique_ptr< compositionalMultiphaseStatistics::StatsAggregator > m_reservoirStatsAggregator; }; diff --git a/src/coreComponents/physicsSolvers/fluidFlow/wells/SinglePhaseWell.cpp b/src/coreComponents/physicsSolvers/fluidFlow/wells/SinglePhaseWell.cpp index 6414693bf22..e6cc5a33092 100644 --- a/src/coreComponents/physicsSolvers/fluidFlow/wells/SinglePhaseWell.cpp +++ b/src/coreComponents/physicsSolvers/fluidFlow/wells/SinglePhaseWell.cpp @@ -27,6 +27,7 @@ #include "constitutive/fluid/singlefluid/SingleFluidSelector.hpp" #include "dataRepository/Group.hpp" #include "mesh/DomainPartition.hpp" +#include "mesh/ElementRegionManager.hpp" #include "mesh/WellElementSubRegion.hpp" #include "mesh/PerforationFields.hpp" #include "mesh/mpiCommunications/CommunicationTools.hpp" @@ -43,7 +44,7 @@ #include "physicsSolvers/fluidFlow/wells/kernels/SinglePhasePerforationFluxKernels.hpp" #include "physicsSolvers/fluidFlow/kernels/singlePhase/FluidUpdateKernel.hpp" #include "physicsSolvers/fluidFlow/kernels/singlePhase/SolutionCheckKernel.hpp" -#include "physicsSolvers/fluidFlow/SinglePhaseStatistics.hpp" +#include "physicsSolvers/fluidFlow/SinglePhaseStatisticsAggregator.hpp" namespace geos { @@ -52,6 +53,19 @@ using namespace dataRepository; using namespace constitutive; using namespace fields; using namespace singlePhaseWellKernels; +using namespace singlePhaseStatistics; + +SinglePhaseBase & getFlowSolver( SinglePhaseWell & wellSolver ) +{ + // TODO: change the way we access the flowSolver here + return wellSolver.getParent().getGroup< SinglePhaseBase >( wellSolver.getFlowSolverName() ); +} + +SinglePhaseBase const & getFlowSolver( SinglePhaseWell const & wellSolver ) +{ + // TODO: change the way we access the flowSolver here + return wellSolver.getParent().getGroup< SinglePhaseBase >( wellSolver.getFlowSolverName() ); +} SinglePhaseWell::SinglePhaseWell( const string & name, Group * const parent ): @@ -143,29 +157,21 @@ void SinglePhaseWell::validateWellConstraints( real64 const & time_n, WellElementSubRegion const & subRegion ) { WellControls & wellControls = getWellControls( subRegion ); + if( !wellControls.useSurfaceConditions() ) { - bool useSeg = wellControls.referenceReservoirRegion().empty(); - GEOS_WARNING_IF( useSeg, - "WellControls referenceReservoirRegion not set and well constraint fluid property calculations will use top segement pressure and temp " ); - if( useSeg ) - { - wellControls.setRegionAveragePressure( -1 ); - } - else - { - // Check if region name exists in list of Reservoir's target regions - string const regionName = wellControls.referenceReservoirRegion(); - SinglePhaseBase const & flowSolver = getParent().getGroup< SinglePhaseBase >( getFlowSolverName() ); - string_array const & targetRegionsNames = flowSolver.getTargetRegionNames(); - auto const pos = std::find( targetRegionsNames.begin(), targetRegionsNames.end(), regionName ); - GEOS_ERROR_IF( pos == targetRegionsNames.end(), - GEOS_FMT( "Region {} is not a target of the reservoir solver and cannot be used for referenceReservoirRegion in WellControl {}.", - regionName, wellControls.getName() ), - getDataContext() ); + bool useSegmentValues = wellControls.referenceReservoirRegion().empty(); + static bool firstNoRefRegionMsg = true; + if( useSegmentValues && firstNoRefRegionMsg ) + { + GEOS_WARNING( WellControls::viewKeyStruct::referenceReservoirRegionString() << + " not set: well constraint fluid property calculations will use top segement pressure and temp ", + wellControls.getDataContext() ); + firstNoRefRegionMsg = false; } } + WellControls::Control currentControl = wellControls.getControl(); real64 const targetTotalRate = wellControls.getTargetTotalRate( time_n ); real64 const targetPhaseRate = wellControls.getTargetPhaseRate( time_n ); @@ -249,7 +255,7 @@ void SinglePhaseWell::updateBHPForConstraint( WellElementSubRegion & subRegion ) } -void SinglePhaseWell::updateVolRateForConstraint( ElementRegionManager const & elemManager, WellElementSubRegion & subRegion ) +void SinglePhaseWell::updateVolRateForConstraint( WellElementSubRegion & subRegion ) { GEOS_MARK_FUNCTION; @@ -263,9 +269,6 @@ void SinglePhaseWell::updateVolRateForConstraint( ElementRegionManager const & e // subRegion data - arrayView1d< real64 const > const pres = - subRegion.getField< well::pressure >(); - arrayView1d< real64 const > const & connRate = subRegion.getField< well::connectionRate >(); @@ -282,39 +285,8 @@ void SinglePhaseWell::updateVolRateForConstraint( ElementRegionManager const & e string const wellControlsName = wellControls.getName(); bool const logSurfaceCondition = isLogLevelActive< logInfo::WellControl >( wellControls.getLogLevel()); integer const useSurfaceConditions = wellControls.useSurfaceConditions(); - real64 flashPressure; - if( useSurfaceConditions ) - { - // use surface conditions - flashPressure = wellControls.getSurfacePressure(); - } - else - { - if( !wellControls.referenceReservoirRegion().empty() ) - { - ElementRegionBase const & region = elemManager.getRegion( wellControls.referenceReservoirRegion() ); - GEOS_ERROR_IF ( !region.hasWrapper( SinglePhaseStatistics::regionStatisticsName()), - GEOS_FMT( "WellControl {} referenceReservoirRegion field requires SinglePhaseStatistics to be configured for region {} ", - wellControls.getName(), wellControls.referenceReservoirRegion() ), - getDataContext() ); - - 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 {} ", - wellControls.getName(), wellControls.referenceReservoirRegion() ), - getDataContext()); - wellControls.setRegionAveragePressure( stats.averagePressure ); - wellControls.setRegionAverageTemperature( stats.averageTemperature ); - } - // use region conditions - flashPressure = wellControls.getRegionAveragePressure(); - if( flashPressure < 0.0 ) - { - // use segment conditions - flashPressure = pres[iwelemRef]; - } - } + ReferenceConditions const refConditions = getReferenceConditions( subRegion ); + real64 & currentVolRate = wellControls.getReference< real64 >( SinglePhaseWell::viewKeyStruct::currentVolRateString() ); @@ -331,13 +303,12 @@ void SinglePhaseWell::updateVolRateForConstraint( ElementRegionManager const & e using COFFSET_WJ = singlePhaseWellKernels::ColOffset_WellJac< IS_THERMAL >; // bring everything back to host, capture the scalars by reference forAll< serialPolicy >( 1, [fluidWrapper, - pres, connRate, dens, dDens, logSurfaceCondition, &useSurfaceConditions, - &flashPressure, + refConditions, ¤tVolRate, dCurrentVolRate, &iwelemRef, @@ -347,28 +318,17 @@ void SinglePhaseWell::updateVolRateForConstraint( ElementRegionManager const & e // - Surface conditions: using the surface pressure provided by the user // - Reservoir conditions: using the pressure in the top element - if( useSurfaceConditions ) + fluidWrapper.update( iwelemRef, 0, refConditions.pressure ); + if( useSurfaceConditions && logSurfaceCondition ) { - // we need to compute the surface density - fluidWrapper.update( iwelemRef, 0, flashPressure ); - if( logSurfaceCondition ) - { - - GEOS_LOG_RANK( GEOS_FMT( "{}: surface density computed with P_surface = {} Pa", - wellControlsName, flashPressure ) ); - } + GEOS_LOG_RANK( GEOS_FMT( "{}: surface density computed with P_surface = {} Pa", + wellControlsName, refConditions.pressure ) ); + } #ifdef GEOS_USE_HIP - GEOS_UNUSED_VAR( wellControlsName ); + GEOS_UNUSED_VAR( wellControlsName ); #endif - } - else - { - real64 const refPres = pres[iwelemRef]; - fluidWrapper.update( iwelemRef, 0, refPres ); - } - real64 const densInv = 1.0 / dens[iwelemRef][0]; currentVolRate = connRate[iwelemRef] * densInv; @@ -388,6 +348,90 @@ void SinglePhaseWell::updateVolRateForConstraint( ElementRegionManager const & e } ); } +void SinglePhaseWell::precomputeReferenceConditions( real64 const time_n, + Group & meshBodies, + MeshBody & meshBody, + WellElementSubRegion const & subRegion ) +{ + WellControls & wellControls = getWellControls( subRegion ); + integer const useSurfaceConditions = wellControls.useSurfaceConditions(); + if( !useSurfaceConditions ) + { + string_view refRegionName = wellControls.referenceReservoirRegion(); + bool const useSegmentValues = refRegionName.empty(); + if( useSegmentValues ) + { + wellControls.setRegionAveragePressure( -1 ); + wellControls.setRegionAverageTemperature( -1 ); + } + else + { + auto & flowSolver = getParent().getGroup< SinglePhaseBase >( getFlowSolverName() ); + MeshLevel & flowMeshLevel = meshBody.getMeshLevel( flowSolver.getDiscretizationName() ); + + if( !m_reservoirStatsAggregator ) + { // lazily initialize the region statistics aggregator + m_reservoirStatsAggregator = std::make_unique< StatsAggregator >( wellControls.getDataContext() ); + m_reservoirStatsAggregator->initStatisticsAggregation( meshBodies, flowSolver ); + m_reservoirStatsAggregator->enableRegionStatisticsAggregation(); + } + + RegionStatistics & stats = m_reservoirStatsAggregator->getRegionStatistics( flowMeshLevel, refRegionName ); + + // compute region stats only if needed (could have already been done for another subRegion) + if( !m_reservoirStatsAggregator->isComputed( time_n, stats ) ) + m_reservoirStatsAggregator->computeRegionsStatistics( time_n ); + + GEOS_WARNING_IF( stats.m_averagePressure <= 0.0, + GEOS_FMT( "No region average quantities computed in reference region '{}'.", + wellControls.referenceReservoirRegion() ), + wellControls.getWrapperDataContext( WellControls::viewKeyStruct::referenceReservoirRegionString() ), + getDataContext() ); + + wellControls.setRegionAveragePressure( stats.m_averagePressure ); + wellControls.setRegionAverageTemperature( stats.m_averageTemperature ); + } + } +} + +SinglePhaseWell::ReferenceConditions +SinglePhaseWell::getReferenceConditions( WellElementSubRegion const & subRegion ) +{ + WellControls & wellControls = getWellControls( subRegion ); + integer const useSurfaceConditions = wellControls.useSurfaceConditions(); + if( useSurfaceConditions ) + { + // use surface conditions + return { + /* .pressure = */ wellControls.getSurfacePressure(), + /* .temperature = */ wellControls.getSurfaceTemperature(), + }; + } + else + { + if( wellControls.getRegionAveragePressure() > 0.0 && wellControls.getRegionAverageTemperature() > 0.0 ) + { // reference region condition properly computed, we can return them + return { + /* .pressure = */ wellControls.getRegionAveragePressure(), + /* .temperature = */ wellControls.getRegionAverageTemperature(), + }; + } + else + { // region average stats not initialized or initialized, fallback to top segment values + GEOS_WARNING( "Region average statistics of reference region not initialized, fallback to top segment values.", + wellControls.getDataContext() ); + + arrayView1d< real64 const > const & pres = subRegion.getField< well::pressure >(); + arrayView1d< real64 const > const & temp = subRegion.getField< well::temperature >(); + localIndex const iwelemRef = subRegion.getTopWellElementIndex(); + return { + /* .pressure = */ pres[iwelemRef], + /* .temperature = */ temp[iwelemRef], + }; + } + } +} + void SinglePhaseWell::updateFluidModel( WellElementSubRegion & subRegion ) const { GEOS_MARK_FUNCTION; @@ -405,11 +449,11 @@ void SinglePhaseWell::updateFluidModel( WellElementSubRegion & subRegion ) const } ); } -real64 SinglePhaseWell::updateSubRegionState( ElementRegionManager const & elemManager, WellElementSubRegion & subRegion ) +real64 SinglePhaseWell::updateSubRegionState( WellElementSubRegion & subRegion ) { // update volumetric rates for the well constraints // Warning! This must be called before updating the fluid model - updateVolRateForConstraint( elemManager, subRegion ); + updateVolRateForConstraint( subRegion ); // update density in the well elements updateFluidModel( subRegion ); @@ -426,12 +470,15 @@ void SinglePhaseWell::initializeWells( DomainPartition & domain, real64 const & GEOS_MARK_FUNCTION; GEOS_UNUSED_VAR( time_n ); + Group & meshBodies = domain.getMeshBodies(); + // loop over the wells - forDiscretizationOnMeshTargets( domain.getMeshBodies(), [&] ( string const &, + forDiscretizationOnMeshTargets( domain.getMeshBodies(), [&] ( string const & meshBodyName, MeshLevel & meshLevel, string_array const & regionNames ) { ElementRegionManager & elemManager = meshLevel.getElemManager(); + MeshBody & meshBody = domain.getMeshBody( meshBodyName ); elemManager.forElementSubRegions< WellElementSubRegion >( regionNames, [&]( localIndex const, @@ -467,7 +514,7 @@ void SinglePhaseWell::initializeWells( DomainPartition & domain, real64 const & if( wellControls.isWellOpen() && !hasNonZeroRate ) { // TODO: change the way we access the flowSolver here - SinglePhaseBase const & flowSolver = getParent().getGroup< SinglePhaseBase >( getFlowSolverName() ); + SinglePhaseBase const & flowSolver = getFlowSolver( *this ); PresTempInitializationKernel::SinglePhaseFlowAccessors resSinglePhaseFlowAccessors( meshLevel.getElemManager(), flowSolver.getName() ); PresTempInitializationKernel::SingleFluidAccessors resSingleFluidAccessors( meshLevel.getElemManager(), flowSolver.getName() ); @@ -495,7 +542,8 @@ 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( elemManager, subRegion ); + precomputeReferenceConditions( time_n, meshBodies, meshBody, subRegion ); + updateSubRegionState( subRegion ); string const & fluidName = subRegion.getReference< string >( viewKeyStruct::fluidNamesString() ); SingleFluidBase & fluid = subRegion.getConstitutiveModel< SingleFluidBase >( fluidName ); @@ -839,7 +887,7 @@ void SinglePhaseWell::computePerforationRates( real64 const & time_n, { // TODO: change the way we access the flowSolver here - SinglePhaseBase const & flowSolver = getParent().getGroup< SinglePhaseBase >( getFlowSolverName() ); + SinglePhaseBase const & flowSolver = getFlowSolver( *this ); PerforationKernel::SinglePhaseFlowAccessors resSinglePhaseFlowAccessors( mesh.getElemManager(), flowSolver.getName() ); PerforationKernel::SingleFluidAccessors resSingleFluidAccessors( mesh.getElemManager(), flowSolver.getName() ); ElementRegionManager & elemManager = mesh.getElemManager(); @@ -1156,7 +1204,7 @@ void SinglePhaseWell::resetStateToBeginningOfStep( DomainPartition & domain ) subRegion.getField< well::connectionRate_n >(); connRate.setValues< parallelDevicePolicy<> >( connRate_n ); - updateSubRegionState( elemManager, subRegion ); + updateSubRegionState( subRegion ); } ); } ); } @@ -1168,11 +1216,14 @@ void SinglePhaseWell::implicitStepSetup( real64 const & time, { WellSolverBase::implicitStepSetup( time, dt, domain ); - forDiscretizationOnMeshTargets( domain.getMeshBodies(), [&] ( string const &, - MeshLevel & mesh, - string_array const & regionNames ) - { + if( m_reservoirStatsAggregator ) + m_reservoirStatsAggregator->setDirty(); // TODO: is it useful? are timestep cut properly managed for this call? + Group & meshBodies = domain.getMeshBodies(); + forDiscretizationOnMeshTargets( meshBodies, [&] ( string const &, + MeshLevel & mesh, + string_array const & regionNames ) + { ElementRegionManager & elemManager = mesh.getElemManager(); elemManager.forElementSubRegions< WellElementSubRegion >( regionNames, [&]( localIndex const, @@ -1198,7 +1249,7 @@ void SinglePhaseWell::implicitStepSetup( real64 const & time, validateWellConstraints( time, dt, subRegion ); - updateSubRegionState( elemManager, subRegion ); + updateSubRegionState( subRegion ); } ); } ); } diff --git a/src/coreComponents/physicsSolvers/fluidFlow/wells/SinglePhaseWell.hpp b/src/coreComponents/physicsSolvers/fluidFlow/wells/SinglePhaseWell.hpp index 8c07e223fb7..b82ebe96342 100644 --- a/src/coreComponents/physicsSolvers/fluidFlow/wells/SinglePhaseWell.hpp +++ b/src/coreComponents/physicsSolvers/fluidFlow/wells/SinglePhaseWell.hpp @@ -36,6 +36,12 @@ namespace constitutive { class SingleFluidBase; } + +namespace singlePhaseStatistics +{ +class StatsAggregator; +} + class WellElementSubRegion; /** @@ -143,7 +149,7 @@ class SinglePhaseWell : public WellSolverBase * @param elemManager the well region manager * @param subRegion the well subregion containing all the primary and dependent fields */ - virtual void updateVolRateForConstraint( ElementRegionManager const & elemManager, WellElementSubRegion & subRegion ); + virtual void updateVolRateForConstraint( WellElementSubRegion & subRegion ); /** * @brief Recompute the BHP pressure that is used in the well constraints @@ -169,7 +175,7 @@ class SinglePhaseWell : public WellSolverBase * @param elemManager the elemManager containing the well * @param subRegion the well subRegion containing the well elements and their associated fields */ - virtual real64 updateSubRegionState( ElementRegionManager const & elemManager, WellElementSubRegion & subRegion ) override; + virtual real64 updateSubRegionState( WellElementSubRegion & subRegion ) override; /** * @brief function to assemble the linear system matrix and rhs @@ -279,6 +285,15 @@ class SinglePhaseWell : public WellSolverBase private: + struct ReferenceConditions + { + real64 pressure; + real64 temperature; + }; + + /// optional statistics aggregator to get the average pressure of simulated region + std::unique_ptr< singlePhaseStatistics::StatsAggregator > m_reservoirStatsAggregator; + virtual void setConstitutiveNames( ElementSubRegionBase & subRegion ) const override; /** @@ -287,6 +302,13 @@ class SinglePhaseWell : public WellSolverBase */ void initializeWells( DomainPartition & domain, real64 const & time_n ) override; + void precomputeReferenceConditions( real64 time_n, + Group & meshBodies, + MeshBody & meshBody, + WellElementSubRegion const & subRegion ); + + ReferenceConditions getReferenceConditions( WellElementSubRegion const & subRegion ); + /** * @brief Make sure that the well constraints are compatible * @param time_n the time at the beginning of the time step diff --git a/src/coreComponents/physicsSolvers/fluidFlow/wells/WellControls.hpp b/src/coreComponents/physicsSolvers/fluidFlow/wells/WellControls.hpp index 301c3c42f7c..a85e2868941 100644 --- a/src/coreComponents/physicsSolvers/fluidFlow/wells/WellControls.hpp +++ b/src/coreComponents/physicsSolvers/fluidFlow/wells/WellControls.hpp @@ -252,7 +252,7 @@ class WellControls : public dataRepository::Group * @brief Getter for the reservoir region associated with reservoir volume constraint * @return name of reservoir region */ - string referenceReservoirRegion() const { return m_referenceReservoirRegion; } + string const & referenceReservoirRegion() const { return m_referenceReservoirRegion; } /** * @brief Getter for the surface pressure when m_useSurfaceConditions == 1 @@ -316,12 +316,14 @@ class WellControls : public dataRepository::Group /** * @brief Getter for the reservoir average pressure when m_useSurfaceConditions == 0 + * @note When not available, value is less or equal to 0.0. * @return the pressure */ real64 getRegionAveragePressure() const { return m_regionAveragePressure; } /** * @brief Set the reservoir average pressure when m_useSurfaceConditions == 0 + * @note When not available, value is less or equal to 0.0. * @param[in] regionAveragePressure value for pressure */ void setRegionAveragePressure( real64 regionAveragePressure ) { m_regionAveragePressure = regionAveragePressure; } diff --git a/src/coreComponents/physicsSolvers/fluidFlow/wells/WellSolverBase.cpp b/src/coreComponents/physicsSolvers/fluidFlow/wells/WellSolverBase.cpp index 060ff895e73..61b9398d0e6 100644 --- a/src/coreComponents/physicsSolvers/fluidFlow/wells/WellSolverBase.cpp +++ b/src/coreComponents/physicsSolvers/fluidFlow/wells/WellSolverBase.cpp @@ -19,7 +19,9 @@ #include "WellSolverBase.hpp" +#include "dataRepository/Group.hpp" #include "mesh/DomainPartition.hpp" +#include "mesh/MeshBody.hpp" #include "mesh/PerforationFields.hpp" #include "mesh/WellElementRegion.hpp" #include "mesh/WellElementSubRegion.hpp" @@ -151,9 +153,10 @@ void WellSolverBase::initializePostSubGroups() { DomainPartition & domain = this->getGroupByPath< DomainPartition >( "/Problem/domain" ); FunctionManager & functionManager = FunctionManager::getInstance(); - forDiscretizationOnMeshTargets( domain.getMeshBodies(), [&] ( string const &, - MeshLevel & mesh, - string_array const & regionNames ) + Group & meshBodies = domain.getMeshBodies(); + forDiscretizationOnMeshTargets( meshBodies, [&] ( string const &, + MeshLevel & mesh, + string_array const & regionNames ) { ElementRegionManager & elemManager = mesh.getElemManager(); elemManager.forElementSubRegions< WellElementSubRegion >( regionNames, @@ -310,7 +313,7 @@ void WellSolverBase::updateState( DomainPartition & domain ) ElementRegionManager & elemManager = mesh.getElemManager(); elemManager.forElementSubRegions< WellElementSubRegion >( regionNames, [&]( localIndex const, WellElementSubRegion & subRegion ) - { updateSubRegionState( elemManager, subRegion ); } ); + { updateSubRegionState( subRegion ); } ); } ); } diff --git a/src/coreComponents/physicsSolvers/fluidFlow/wells/WellSolverBase.hpp b/src/coreComponents/physicsSolvers/fluidFlow/wells/WellSolverBase.hpp index 04fe58112b4..0f9bd2615b6 100644 --- a/src/coreComponents/physicsSolvers/fluidFlow/wells/WellSolverBase.hpp +++ b/src/coreComponents/physicsSolvers/fluidFlow/wells/WellSolverBase.hpp @@ -253,7 +253,7 @@ class WellSolverBase : public PhysicsSolverBase * @param elemManager the elemManager containing the well * @param subRegion the well subRegion containing the well elements and their associated fields */ - virtual real64 updateSubRegionState( ElementRegionManager const & elemManager, WellElementSubRegion & subRegion ) = 0; + virtual real64 updateSubRegionState( WellElementSubRegion & subRegion ) = 0; /** * @brief Recompute the perforation rates for all the wells