diff --git a/inputFiles/compositionalMultiphaseWell/simpleCo2InjTutorial_base.xml b/inputFiles/compositionalMultiphaseWell/simpleCo2InjTutorial_base.xml index 705b4934fb9..bcf4ce1e61b 100644 --- a/inputFiles/compositionalMultiphaseWell/simpleCo2InjTutorial_base.xml +++ b/inputFiles/compositionalMultiphaseWell/simpleCo2InjTutorial_base.xml @@ -191,15 +191,38 @@ + + + + + + + + + + + numElementsPerSegment="2" + logLevel="1"> @@ -34,6 +35,29 @@ + + + + + + + + + + + + + + + + + + diff --git a/src/coreComponents/physicsSolvers/LogLevelsInfo.hpp b/src/coreComponents/physicsSolvers/LogLevelsInfo.hpp index 73037d8bb04..b65b3d87af8 100644 --- a/src/coreComponents/physicsSolvers/LogLevelsInfo.hpp +++ b/src/coreComponents/physicsSolvers/LogLevelsInfo.hpp @@ -126,6 +126,12 @@ struct Statistics static constexpr std::string_view getDescription() { return "Print statistics"; } }; +struct StatisticsPerRegion +{ + static constexpr int getMinLogLevel() { return 2; } + static constexpr std::string_view getDescription() { return "Print per region statistics"; } +}; + struct SurfaceGenerator { static constexpr int getMinLogLevel() { return 1; } diff --git a/src/coreComponents/physicsSolvers/StatisticsAggregatorBase.hpp b/src/coreComponents/physicsSolvers/StatisticsAggregatorBase.hpp index 64c52358bc7..b4ad6706ff4 100644 --- a/src/coreComponents/physicsSolvers/StatisticsAggregatorBase.hpp +++ b/src/coreComponents/physicsSolvers/StatisticsAggregatorBase.hpp @@ -47,12 +47,15 @@ class RegionStatisticsBase : public dataRepository::Group * @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 statsOutputEnabled If true, the stats are saved in the output HDF5 - * (through dataRepository::RestartFlags, but not functional for this output for now). + * @param setName Optional name of a set, restricting the statistics on given mesh element. + * If none is specified (empty), "all" elements set is targeted. + * @param dataOutputEnabled If true, the stats are saved in the output HDF5 (dataRepository::RestartFlags, not functional for this output + * for now). */ RegionStatisticsBase( string const & targetName, dataRepository::Group * const parent, - bool statsOutputEnabled ); + string_view setName, + bool dataOutputEnabled ); /** * @return the name of the data-repository object that is targeted by the statistics @@ -61,6 +64,16 @@ class RegionStatisticsBase : public dataRepository::Group string_view getTargetName() const { return getName(); } + /** + * @return Optional name of a set, restricting the statistics on given mesh element. + * If none is specified (empty), "all" elements set is targeted. + */ + string_view getSetName() const + { return m_setName; } + +private: + /// @see getSetName() + string const m_setName; }; template< typename T > @@ -82,6 +95,10 @@ class StatsAggregatorBase using StatsGroupType = typename StatsAggregatorTraits< Impl >::StatsGroupType; + using SetType = SortedArray< localIndex >; + + using SetViewType = SortedArrayView< localIndex const >; + /** * @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, @@ -90,26 +107,41 @@ class StatsAggregatorBase * @tparam OwnerType the concrete type of the OwnerType param */ template< typename OwnerType > - using RegionStatsFunc = std::function< void ( OwnerType &, + using RegionStatsFunc = std::function< void ( OwnerType, StatsGroupType & ) >; /** * @brief A functor that can be used to register a statistics Group instance. Parameters: - * - dataRepository::Group &, the parent Group, - * - string const &, name of the statistics target (i.e. region name). + * - dataRepository::Group & targetName: the parent Group, + * - string const & parent: name of the statistics target (i.e. region name). + * - string_view setName: name of the target element set. If none is specified (empty), "all" elements set is targeted. + * - bool dataOutputEnabled: enable Group output features */ using RegionStatsRegisterFunc = std::function< StatsGroupType & ( dataRepository::Group &, - string const & ) >; + string const &, + string_view, + bool ) >; /** * @brief the associated view keys */ struct ViewKeys { - /// String for the discretization statistics group - constexpr static char const * statisticsString() { return "statistics"; } + /// String for the mesh-level statistics group + constexpr static string_view statisticsString() { return "statistics"; } /// String for the region statistics group - constexpr static char const * regionsStatisticsString() { return "regionsStatistics"; } + constexpr static string_view setsStatisticsString() { return "setsStatistics"; } + /// on-purpose generated compound set name + constexpr static string_view compoundSetNameString() { return "__compound"; } + }; + + /** + * @brief Allow to reference a set in a given MeshLevel structure. + */ + struct MeshLevelSet + { + MeshLevel & mesh; + string_view setName; }; /** @@ -121,8 +153,7 @@ class StatsAggregatorBase bool statsOutputEnabled ); /** - * @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 + * @brief Enable the computation of any statistics, explores the data structure to collect them. * @param solver flow solver object to retrieve: - the simulated regions, - fields for statistics computation. @@ -131,15 +162,98 @@ class StatsAggregatorBase void initStatisticsAggregation( dataRepository::Group & meshBodies, SolverType & solver ); - void forRegionStatistics( RegionStatsFunc< MeshLevel > const & functor ) const; + /** + * @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 ); + + /** + * @name Data Structure Traversal Strategy + */ + ///@{ + + /** + * @brief Execute the given functor on each mesh-level (discretization) available for this instance. + * @param functor the function to execute. + */ + void forRegionStatistics( RegionStatsFunc< MeshLevel & > const & functor ) const; + /** + * @brief Execute the given functor on each set present in the given mesh-level (discretization). + * @param mesh the target mesh-level. + * @param meshSetsStats the instance mesh-level statistics data structure. + * @param enableCompoundSet if true, the functor will also be applied on the compound set (which + * is an optional internal computation data structure). + * @param functor the function to execute. + */ void forRegionStatistics( MeshLevel & mesh, - StatsGroupType & meshRegionsStatistics, - RegionStatsFunc< CellElementRegion > const & functor ) const; + StatsGroupType & meshSetsStats, + bool enableCompoundSet, + RegionStatsFunc< MeshLevelSet > const & functor ) const; + /** + * @brief Execute the given functor on each region available in the given mesh-level for the given set. + * @param meshSet the target mesh-level set. + * @param setStats the mesh-level statistics data structure for the given set. + * @param functor the function to execute. + */ + void forRegionStatistics( MeshLevelSet meshSet, + StatsGroupType & setStats, + RegionStatsFunc< CellElementRegion & > const & functor ) const; + + /** + * @brief Execute the given functor on each sub-region available in the given region for the given set. + * @param meshSet the target region set. + * @param setRegionStats the region statistics data structure for the given set. + * @param functor the function to execute. + */ void forRegionStatistics( CellElementRegion & region, - StatsGroupType & regionStatistics, - RegionStatsFunc< CellElementSubRegion > const & functor ) const; + StatsGroupType & setRegionStats, + RegionStatsFunc< CellElementSubRegion & > const & functor ) const; + + /** + * @brief Get the Instance Statistics Group for the given mesh-level. + * @param mesh the mesh-level (discretization). + * @return A Group instance hierarchically holding all statistics data structures. + */ + dataRepository::Group & getInstanceStatisticsGroup( MeshLevel & mesh ) const; + + /** + * @return A statistics data structure (see StatsGroupType), aggregated from all targeted regions + * and sets in the given mesh-level. + * @param mesh the mesh-level (discretization). + */ + StatsGroupType & getMeshLevelStatistics( MeshLevel & mesh ) const; + + /** + * @return A statistics data structure (see StatsGroupType), aggregated from all targeted regions + * for the given set in the given mesh-level. + * @param mesh the mesh-level (discretization). + * @param setName the name of element set. If none is specified (empty), "all" elements set is targeted. + * @throw InputError if "all" is targeted, but the instance does not target "all" element set. + */ + StatsGroupType & getSetStatistics( MeshLevel & mesh, string_view setName = "" ) const; + + /** + * @return A statistics data structure (see StatsGroupType), for the given set, in the given region, + * in the given mesh-level. + * @param mesh The mesh-level (discretization) + * @param setName the name of element set. If none is specified (empty), "all" elements set is targeted. + * @param regionName The name of the desired region + * @throw InputError - if no statistics data is found for the given region name. + * - if "all" is targeted, but the instance does not target "all" element set. + */ + StatsGroupType & getRegionStatistics( MeshLevel & mesh, string_view regionName, string_view setName = "" ) const; + + ///@} /** * @param[in] timeRequest The time for which we want to know if the statistics are computed. @@ -149,17 +263,18 @@ class StatsAggregatorBase bool isComputed( real64 const timeRequest, StatsGroupType const & stats ); /** - * @brief set the statistics as dirty, ensuring isComputed() will be false until the next computation. + * @return true if the region statistics target given sets. + * false if the statistics are targeted on only one mesh-level-wide set */ - void setDirty(); + bool isRestrictedToSets() + {return !( m_setNames.empty() || ( m_setNames.size() == 1 && m_setNames.count( "all" ) > 0 ) ); } /** - * @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. + * @return true if the region statistics target multiple sets. + * false if the statistics are targeted on only one set */ - bool computeRegionsStatistics( real64 const timeRequest ); + bool isTargetingMultipleSets() + { return m_setNames.size() > 1; } /** * @return the name of the entity that needs the statistics. @@ -173,18 +288,6 @@ class StatsAggregatorBase stdVector< string > const & getWarnings() const { return m_warnings; } - dataRepository::Group & getInstanceStatisticsGroup( MeshLevel & mesh ) const; - - StatsGroupType & getRegionsStatistics( MeshLevel & mesh ) const; - - /** - * @return a specific statistics Group instance. - * @param mesh The desired mesh-level - * @param regionName The name of the desired region - * @throw InputError if no statistics data is found for the given region name. - */ - StatsGroupType & getRegionStatistics( MeshLevel & mesh, string_view regionName ) const; - protected: struct StatsState @@ -193,21 +296,39 @@ class StatsAggregatorBase bool m_isDirty = false; }; - struct DiscretizationGroupPath + using SetsCompound = stdMap< string, stdMap< string, SetType > >; + + struct DiscretizationSetPath { + /// The id of the mesh body in the "meshBodies" Group localIndex m_meshBody; + /// The id of the mesh level in the mesh body. localIndex m_meshLevel; + /// The names of this mesh-level element sets. If none is specified (empty), "all" elements set is targeted. + string_array m_setNames; + /// the regions names in the mesh-level / mesh level string_array m_regionNames; + /// if at least one set intersects with another, a compound of (local element index) set is needed to compute the mesh-level statistics. + SetsCompound m_setsCompound; }; /// @see getOwnerName() dataRepository::DataContext const & m_ownerDataContext; - bool const m_statsOutputEnabled; + bool const m_dataOutputEnabled; dataRepository::Group * m_meshBodies = nullptr; - std::vector< DiscretizationGroupPath > m_discretizationsPaths; + /// Data repository paths of the element sets per mesh levels + /// TODO: can it be generalized with the "all" set + std::vector< DiscretizationSetPath > m_discretizationsPaths; + + /// The list of mesh element sets to restrict the region statistics to. + /// Cannot be empty: if the whole mesh-level is to process, "all" must be targeted. + std::set< string > m_setNames; + + /// if true, at least ont set intersects with another, threfore we must compute mesh-level statistics with a coupound. + bool m_isAnySetsIntersecting; /// @see getWarnings() stdVector< string > m_warnings; @@ -220,14 +341,43 @@ class StatsAggregatorBase * 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 + * @param setNames The list of mesh element sets to restrict the statistics to. + * If empty, the whole mesh-level is processed. */ - void enableRegionStatisticsAggregation( RegionStatsRegisterFunc && registerStatsFunc ); + void enableRegionStatisticsAggregation( RegionStatsRegisterFunc && registerStatsFunc, + string_array const & setNames ); /** - * @param path the path of the discretization group in the data-repository. + * @param path the path of the mesh-level group in the data-repository. * @return MeshLevel& the MeshLevel Group for the given discretisation */ - MeshLevel & getMeshLevel( DiscretizationGroupPath const & path ) const; + MeshLevel & getMeshLevel( DiscretizationSetPath const & path ) const; + + /** + * @return An optional statistics data structure (see StatsGroupType), aggregated from all sets + * over all targeted regions, for the given set in the given mesh-level. + * @param mesh the mesh-level (discretization). + * @throw InputError if no compound set exists, for instance if isSetsCompoundEnabled == false. + */ + StatsGroupType & getCompoundSetStatistics( MeshLevel & mesh ) const; + + /** + * @return true if the RegionStatistics is processed from a compound set of other instances. + * @param stats a statistics data structure + */ + bool isCompoundSetStatistics( StatsGroupType const & stats ) const; + + /** + * @brief TODO + */ + static std::set< string > getMeshLevelPartialSetNames( MeshLevel const & mesh, + string_array const & regionNames ); + + /** + * @brief Set the statstics target sets to restrict the statistics to. + * @param setNames the list of target sets. If "all" is included, or none is present, all the discretisation is targeted. + */ + void setTargetSets( string_array const & setNames ); /** * @brief Initialize all statistics values to aggregable default values, @@ -243,10 +393,13 @@ class StatsAggregatorBase * @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 + * @param targetSet * @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 ); } + void computeSubRegionRankStats( CellElementSubRegion & subRegion, + StatsGroupType & subRegionStats, + SetType const & targetSet ) const + { static_cast< Impl const * >(this)->computeSubRegionRankStats( subRegion, subRegionStats, targetSet ); } /** * @brief Aggregate all instance statistics with those of another instance on the current rank. diff --git a/src/coreComponents/physicsSolvers/StatisticsAggregatorBaseHelpers.hpp b/src/coreComponents/physicsSolvers/StatisticsAggregatorBaseHelpers.hpp index aaf2165111e..3ff6a1e7b95 100644 --- a/src/coreComponents/physicsSolvers/StatisticsAggregatorBaseHelpers.hpp +++ b/src/coreComponents/physicsSolvers/StatisticsAggregatorBaseHelpers.hpp @@ -33,11 +33,13 @@ namespace geos inline RegionStatisticsBase::RegionStatisticsBase( string const & targetName, dataRepository::Group * const parent, - bool const statsOutputEnabled ): + string_view setName, + bool const dataOutputEnabled ): dataRepository::Group( targetName, parent ), - m_time( std::numeric_limits< double >::lowest() ) + m_time( std::numeric_limits< double >::lowest() ), + m_setName( setName ) { - Group::setRestartFlags( statsOutputEnabled ? + Group::setRestartFlags( dataOutputEnabled ? dataRepository::RestartFlags::WRITE_AND_READ : dataRepository::RestartFlags::NO_WRITE ); @@ -46,11 +48,35 @@ inline RegionStatisticsBase::RegionStatisticsBase( string const & targetName, template< typename Impl > StatsAggregatorBase< Impl >::StatsAggregatorBase( dataRepository::DataContext const & ownerDataContext, - bool const statsOutputEnabled ): + bool const dataOutputEnabled ): m_ownerDataContext( ownerDataContext ), - m_statsOutputEnabled( statsOutputEnabled ) + m_dataOutputEnabled( dataOutputEnabled ), + m_isAnySetsIntersecting( false ) {} +template< typename Impl > +std::set< string > +StatsAggregatorBase< Impl >::getMeshLevelPartialSetNames( MeshLevel const & mesh, + string_array const & regionNames ) +{ + std::set< string > foundSetNames; + + mesh.getElemManager().forElementRegions< CellElementRegion >( regionNames, + [&] ( localIndex &, + CellElementRegion const & region ) + { + region.forElementSubRegions< CellElementSubRegion >( [&] ( CellElementSubRegion const & subRegion ) { + dataRepository::Group const & sets = subRegion.sets(); + sets.forWrappers< SetType >( [&] ( dataRepository::Wrapper< SetType > const & set ) + { + foundSetNames.emplace( set.getName() ); + } ); + } ); + } ); + + return foundSetNames; +} + template< typename Impl > void StatsAggregatorBase< Impl >::initStatisticsAggregation( dataRepository::Group & meshBodies, @@ -63,9 +89,9 @@ StatsAggregatorBase< Impl >::initStatisticsAggregation( dataRepository::Group & string_array const & regionNames ) { // getting the container of all requesters statistics groups (can be already initialized) - dataRepository::Group * meshStatsGroup = mesh.getGroupPointer( ViewKeys::statisticsString() ); + dataRepository::Group * meshStatsGroup = mesh.getGroupPointer( string( ViewKeys::statisticsString()) ); if( meshStatsGroup == nullptr ) - meshStatsGroup = &mesh.registerGroup( ViewKeys::statisticsString() ); + meshStatsGroup = &mesh.registerGroup( string( ViewKeys::statisticsString() ) ); // registering the container of instance statistics groups (must be unique for this instance) string const & ownerName = getOwnerName(); @@ -75,99 +101,291 @@ StatsAggregatorBase< Impl >::initStatisticsAggregation( dataRepository::Group & m_ownerDataContext ); meshStatsGroup->registerGroup( ownerName ); - // remembering the path of this discretization MeshBody const & body = m_meshBodies->getGroup< MeshBody >( meshBodyName ); - DiscretizationGroupPath const path { - /* .m_meshBody = */ body.getIndexInParent(), - /* .m_meshLevel = */ mesh.getIndexInParent(), - /* .m_regionNames = */ regionNames, - }; - m_discretizationsPaths.push_back( path ); + + /// finding all sets in this mesh level + std::set< string > foundSetNames = getMeshLevelPartialSetNames( mesh, regionNames ); + + { // remembering the path of this discretization + DiscretizationSetPath const path { + /* .m_meshBody = */ body.getIndexInParent(), + /* .m_meshLevel = */ mesh.getIndexInParent(), + /* .m_setNames = */ string_array( foundSetNames.begin(), foundSetNames.end()), + /* .m_regionNames = */ regionNames, + /* .m_setsCompound = */ {}, + }; + m_discretizationsPaths.emplace_back( path ); + } } ); } template< typename Impl > void -StatsAggregatorBase< Impl >::enableRegionStatisticsAggregation( RegionStatsRegisterFunc && registerStatsFunc ) +StatsAggregatorBase< Impl >::enableRegionStatisticsAggregation( RegionStatsRegisterFunc && registerStatsFunc, + string_array const & setNames ) { if( m_meshBodies == nullptr ) return; - integer regionCount = 0; - integer subRegionCount = 0; + m_setNames.insert( setNames.begin(), setNames.end() ); + if( m_setNames.empty()) + m_setNames.emplace( "all" ); - for( auto const & path : m_discretizationsPaths ) + bool regionFound = false; + bool subRegionFound = false; + std::set< string > confirmedSets; + + for( auto & path : m_discretizationsPaths ) { MeshLevel & mesh = getMeshLevel( path ); ElementRegionManager & elemManager = mesh.getElemManager(); dataRepository::Group & statisticsGroup = getInstanceStatisticsGroup( mesh ); - StatsGroupType & meshRegionsStats = registerStatsFunc( statisticsGroup, - ViewKeys::regionsStatisticsString() ); + StatsGroupType & meshStats = registerStatsFunc( statisticsGroup, string( ViewKeys::setsStatisticsString() ), + "", m_dataOutputEnabled ); + SetsCompound meshSetsCompounds; + bool isAnySetIntersects = false; - for( size_t i = 0; i < path.m_regionNames.size(); ++i ) + // registering all stats Group in the dataRepository for each set + for( string const & setName : m_setNames ) { - CellElementRegion & region = elemManager.getRegion< CellElementRegion >( path.m_regionNames[i] ); - StatsGroupType & regionStats = registerStatsFunc( meshRegionsStats, - region.getName() ); + StatsGroupType & meshSetStats = registerStatsFunc( meshStats, setName, + setName, m_dataOutputEnabled ); - region.forElementSubRegions< CellElementSubRegion >( [&] ( CellElementSubRegion & subRegion ) + for( size_t regionId = 0; regionId < path.m_regionNames.size(); ++regionId ) { - registerStatsFunc( regionStats, - subRegion.getName() ); - ++subRegionCount; - } ); - ++regionCount; + CellElementRegion & region = elemManager.getRegion< CellElementRegion >( path.m_regionNames[regionId] ); + StatsGroupType & setRegionStats = registerStatsFunc( meshSetStats, region.getName(), + setName, m_dataOutputEnabled ); + stdMap< string, SetType > & regionSetsCompounds = meshSetsCompounds.get_inserted( region.getName() ); + regionFound = true; + + region.forElementSubRegions< CellElementSubRegion >( [&] ( CellElementSubRegion & subRegion ) + { + registerStatsFunc( setRegionStats, subRegion.getName(), + setName, false ); + subRegionFound = true; + + auto const * const subRegionSetWrapper = subRegion.sets().getWrapperPointer< SetType >( setName ); + if( subRegionSetWrapper != nullptr ) + { + confirmedSets.emplace( setName ); + + // Insert the set elements ids in the mesh-level compound. If doubles are found, sets intersect. + SetViewType const & subRegionSet = subRegionSetWrapper->reference(); + arrayView1d< integer const > const elemGhostRank = subRegion.ghostRank(); + SetType & subRegionSetsCompound = regionSetsCompounds.get_inserted( subRegion.getName() ); + for( localIndex setIter = 0; setIter < subRegionSet.size(); ++setIter ) + { + localIndex elemId = subRegionSet[setIter]; + if( elemGhostRank[elemId] < 0 ) + { + if( !subRegionSetsCompound.insert( elemId ) ) + isAnySetIntersects = true; + } + } + } + } ); + } + } + + // if needed, registering stats Group in the dataRepository for the compound set + isAnySetIntersects = (bool)( MpiWrapper::max( (unsigned char)( isAnySetIntersects ) ) ); + if( isAnySetIntersects ) + { + string_view setName = ViewKeys::compoundSetNameString(); + StatsGroupType & meshSetStats = registerStatsFunc( meshStats, string( setName ), + setName, false ); + for( size_t regionId = 0; regionId < path.m_regionNames.size(); ++regionId ) + { + CellElementRegion & region = elemManager.getRegion< CellElementRegion >( path.m_regionNames[regionId] ); + StatsGroupType & setRegionStats = registerStatsFunc( meshSetStats, region.getName(), + setName, false ); + + region.forElementSubRegions< CellElementSubRegion >( [&] ( CellElementSubRegion & subRegion ) + { + registerStatsFunc( setRegionStats, subRegion.getName(), + setName, false ); + + path.m_setsCompound + .get_inserted( region.getName()) + .insert( meshSetsCompounds + .at( region.getName() ) + .extract( subRegion.getName() ) ); + } ); + } } } - GEOS_ERROR_IF( regionCount == 0 || subRegionCount == 0, + m_regionStatsState.m_isEnabled = true; + m_regionStatsState.m_isDirty = true; + + GEOS_ERROR_IF( regionFound == 0 || subRegionFound == 0, GEOS_FMT( "Missing region for computing statistics:\n- {} regions,\n- {} sub-regions.", - getOwnerName(), regionCount, subRegionCount ), + getOwnerName(), + regionFound ? "found" : "no", + subRegionFound ? "found" : "no" ), m_ownerDataContext ); - m_regionStatsState.m_isEnabled = true; - m_regionStatsState.m_isDirty = true; + // at this point, m_setNames content is element sets user *requests*: they are the same on every ranks, + // but we don't know yet if any of these set is missing or empty (=> confirmedSets). + std::set< string > notFoundSets; + for( string const & requestedSetName : m_setNames ) + { + if( MpiWrapper::max( confirmedSets.count( requestedSetName ) ) == 0 ) + notFoundSets.emplace( requestedSetName ); + } + if( MpiWrapper::commRank() == 0 ) + GEOS_ERROR_IF( !notFoundSets.empty(), + GEOS_FMT( "During retrieval of statistics element sets, the following set(s) could not be " + "found or does not contain any element:\n- {}", + stringutilities::join( notFoundSets, "\n- " ) ), + m_ownerDataContext ); +} + +template< typename Impl > +MeshLevel & +StatsAggregatorBase< Impl >::getMeshLevel( DiscretizationSetPath const & path ) const +{ + MeshBody & body = m_meshBodies->getGroup< MeshBody >( path.m_meshBody ); + MeshLevel & mesh = body.getMeshLevel( path.m_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( string( ViewKeys::statisticsString() ) ); + dataRepository::Group & instanceStatisticsGroup = meshStatsGroup.getGroup( getOwnerName() ); + return instanceStatisticsGroup; +} + +template< typename Impl > +typename StatsAggregatorBase< Impl >::StatsGroupType & +StatsAggregatorBase< Impl >:: +getMeshLevelStatistics( MeshLevel & mesh ) const +{ + // considering everything is initialized, or else, crash gracefully + dataRepository::Group & instanceStatisticsGroup = getInstanceStatisticsGroup( mesh ); + return instanceStatisticsGroup.getGroup< StatsGroupType >( string( ViewKeys::setsStatisticsString()) ); +} + +template< typename Impl > +typename StatsAggregatorBase< Impl >::StatsGroupType & +StatsAggregatorBase< Impl >::getSetStatistics( MeshLevel & mesh, + string_view requestedSetName ) const +{ + // considering mesh-level stats structure is initialized, or else, crash gracefully + StatsGroupType & meshStats = getMeshLevelStatistics( mesh ); + string const setName = string( requestedSetName.empty() ? "all" : requestedSetName ); + StatsGroupType * const stats = meshStats.template getGroupPointer< StatsGroupType >( setName ); + GEOS_THROW_IF( stats == nullptr, + GEOS_FMT( "Element set '{}' statistics data structure not found, has the set been requested?\nRequested element sets:\n- {}", + setName, stringutilities::join( m_setNames, "\n- " ) ), + InputError, m_ownerDataContext ); + return *stats; +} + +template< typename Impl > +typename StatsAggregatorBase< Impl >::StatsGroupType & +StatsAggregatorBase< Impl >::getCompoundSetStatistics( MeshLevel & mesh ) const +{ + // considering mesh-level stats structure is initialized, or else, crash gracefully + StatsGroupType & meshStats = getMeshLevelStatistics( mesh ); + string const setName = string( ViewKeys::compoundSetNameString() ); + StatsGroupType * const stats = meshStats.template getGroupPointer< StatsGroupType >( setName ); + GEOS_THROW_IF( stats == nullptr, + GEOS_FMT( "Compound element set data structure not found ({}), is the aggregator initialized? Was the compound needed?\nRequested element sets:\n- {}", + setName, stringutilities::join( m_setNames, "\n- " ) ), + InputError, m_ownerDataContext ); + return *stats; +} + +template< typename Impl > +bool +StatsAggregatorBase< Impl >::isCompoundSetStatistics( StatsGroupType const & setStats ) const +{ + return setStats.getSetName() == ViewKeys::compoundSetNameString(); +} + +template< typename Impl > +typename StatsAggregatorBase< Impl >::StatsGroupType & +StatsAggregatorBase< Impl >::getRegionStatistics( MeshLevel & mesh, + string_view setName, + string_view regionName ) const +{ + StatsGroupType & setStats = getSetStatistics( mesh, setName ); + StatsGroupType * const stats = setStats.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( setStats.getSubGroupsNames(), "\n- " ) ), + InputError, m_ownerDataContext ); + return *stats; } template< typename Impl > void -StatsAggregatorBase< Impl >::forRegionStatistics( RegionStatsFunc< MeshLevel > const & func ) const +StatsAggregatorBase< Impl >::forRegionStatistics( RegionStatsFunc< MeshLevel & > const & func ) const { for( auto const & path : m_discretizationsPaths ) { MeshLevel & mesh = getMeshLevel( path ); - StatsGroupType & meshRegionsStats = getRegionsStatistics( mesh ); + StatsGroupType & meshLevelStats = getMeshLevelStatistics( mesh ); - func( mesh, meshRegionsStats ); + func( mesh, meshLevelStats ); } } template< typename Impl > void StatsAggregatorBase< Impl >::forRegionStatistics( MeshLevel & mesh, - StatsGroupType & meshRegionsStatistics, - RegionStatsFunc< CellElementRegion > const & func ) const + StatsGroupType & meshLevelStats, + bool enableCompoundSet, + RegionStatsFunc< MeshLevelSet > const & func ) const +{ + meshLevelStats.template forSubGroups< StatsGroupType >( [&] ( StatsGroupType & setStats ) + { + if( enableCompoundSet || setStats.getSetName() != ViewKeys::compoundSetNameString() ) + { + MeshLevelSet meshSet { + /* .mesh = */ mesh, + /* .setName = */ setStats.getTargetName(), + }; + + func( meshSet, setStats ); + } + } ); +} + +template< typename Impl > +void +StatsAggregatorBase< Impl >::forRegionStatistics( MeshLevelSet const meshSet, + StatsGroupType & setStats, + RegionStatsFunc< CellElementRegion & > const & func ) const { - ElementRegionManager & elemManager = mesh.getElemManager(); - meshRegionsStatistics.template forSubGroups< StatsGroupType >( [&] ( StatsGroupType & regionStatistics ) + ElementRegionManager & elemManager = meshSet.mesh.getElemManager(); + setStats.template forSubGroups< StatsGroupType >( [&] ( StatsGroupType & setRegionStats ) { - string_view targetName = regionStatistics.getTargetName(); - CellElementRegion & region = elemManager.getRegion< CellElementRegion >( string( targetName ) ); + string_view regionName = setRegionStats.getTargetName(); + CellElementRegion & region = elemManager.getRegion< CellElementRegion >( string( regionName ) ); - func( region, regionStatistics ); + func( region, setRegionStats ); } ); } template< typename Impl > void -StatsAggregatorBase< Impl >::forRegionStatistics( CellElementRegion & region, - StatsGroupType & regionStatistics, - RegionStatsFunc< CellElementSubRegion > const & func ) const +StatsAggregatorBase< Impl >::forRegionStatistics( CellElementRegion & setRegion, + StatsGroupType & setRegionStats, + RegionStatsFunc< CellElementSubRegion & > const & func ) const { - regionStatistics.template forSubGroups< StatsGroupType >( [&] ( StatsGroupType & subRegionStatistics ) + setRegionStats.template forSubGroups< StatsGroupType >( [&] ( StatsGroupType & subRegionStatistics ) { - string_view targetName = subRegionStatistics.getTargetName(); - CellElementSubRegion & subRegion = region.getSubRegion< CellElementSubRegion >( string( targetName ) ); + string_view subRegionName = subRegionStatistics.getTargetName(); + CellElementSubRegion & subRegion = setRegion.getSubRegion< CellElementSubRegion >( string( subRegionName ) ); + func( subRegion, subRegionStatistics ); } ); } @@ -200,101 +418,75 @@ StatsAggregatorBase< Impl >::computeRegionsStatistics( real64 const timeRequest m_warnings.clear(); - // computation of sub region stats - forRegionStatistics( [&, timeRequest] ( MeshLevel & mesh, StatsGroupType & meshRegionsStats ) + for( auto const & path : m_discretizationsPaths ) { - forRegionStatistics( mesh, - meshRegionsStats, - [&, timeRequest] ( CellElementRegion & region, StatsGroupType & regionStats ) + MeshLevel & mesh = getMeshLevel( path ); + StatsGroupType & meshLevelStats = getMeshLevelStatistics( mesh ); + bool const isSetsCompoundEnabled = !path.m_setsCompound.empty(); + + // computation of sub region stats for each selected set + initStats( meshLevelStats, timeRequest ); + forRegionStatistics( mesh, meshLevelStats, true, + [&, timeRequest] ( MeshLevelSet meshSet, StatsGroupType & setStats ) { - forRegionStatistics( region, - regionStats, - [&, timeRequest] ( CellElementSubRegion & subRegion, StatsGroupType & subRegionStats ) + initStats( setStats, timeRequest ); + forRegionStatistics( meshSet, setStats, + [&, timeRequest] ( CellElementRegion & region, StatsGroupType & regionStats ) { - initStats( subRegionStats, timeRequest ); - computeSubRegionRankStats( subRegion, subRegionStats ); + initStats( regionStats, timeRequest ); + forRegionStatistics( region, regionStats, + [&, timeRequest] ( CellElementSubRegion & subRegion, StatsGroupType & subRegionStats ) + { + initStats( subRegionStats, timeRequest ); + + SetType const & targetSet = meshSet.setName == ViewKeys::compoundSetNameString() ? + path.m_setsCompound.at( region.getName() ).at( subRegion.getName() ) : + subRegion.sets().getReference< SetType >( string( meshSet.setName ) ); + computeSubRegionRankStats( subRegion, subRegionStats, targetSet ); + } ); } ); } ); - } ); - - // aggregation of computations from the sub regions - forRegionStatistics( [&, timeRequest] ( MeshLevel & mesh, StatsGroupType & meshRegionsStats ) - { - initStats( meshRegionsStats, timeRequest ); - forRegionStatistics( mesh, - meshRegionsStats, - [&, timeRequest] ( CellElementRegion & region, StatsGroupType & regionStats ) + // aggregation of computations from the sub regions + forRegionStatistics( mesh, meshLevelStats, true, + [&, timeRequest] ( MeshLevelSet meshSet, StatsGroupType & setStats ) { - initStats( regionStats, timeRequest ); - - forRegionStatistics( region, - regionStats, - [&] ( CellElementSubRegion &, StatsGroupType & subRegionStats ) + forRegionStatistics( meshSet, setStats, + [&, timeRequest] ( CellElementRegion & region, StatsGroupType & regionStats ) { - aggregateStats( regionStats, subRegionStats ); - - mpiAggregateStats( subRegionStats ); - postAggregateStats( subRegionStats ); + forRegionStatistics( region, regionStats, + [&] ( CellElementSubRegion &, StatsGroupType & subRegionStats ) + { + aggregateStats( regionStats, subRegionStats ); + // sub-region stats finalization is disabled, it does not seem useful + // mpiAggregateStats( subRegionStats ); + // postAggregateStats( subRegionStats ); + } ); + + aggregateStats( setStats, regionStats ); + + mpiAggregateStats( regionStats ); + postAggregateStats( regionStats ); } ); - aggregateStats( meshRegionsStats, regionStats ); + // if there is no compound set, we can simply aggregate each set statistics to the mesh-level statistics. + // if there is one, we will only take into account the statistics of the compound set. + if( !isSetsCompoundEnabled || isCompoundSetStatistics( setStats ) ) + { + aggregateStats( meshLevelStats, setStats ); + } - mpiAggregateStats( regionStats ); - postAggregateStats( regionStats ); + mpiAggregateStats( setStats ); + postAggregateStats( setStats ); } ); - mpiAggregateStats( meshRegionsStats ); - postAggregateStats( meshRegionsStats ); - } ); + mpiAggregateStats( meshLevelStats ); + postAggregateStats( meshLevelStats ); + } 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.m_meshBody ); - MeshLevel & mesh = body.getMeshLevel( path.m_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/CompositionalMultiphaseStatisticsAggregator.cpp b/src/coreComponents/physicsSolvers/fluidFlow/CompositionalMultiphaseStatisticsAggregator.cpp index 29e1e1c5fcd..17fc9191292 100644 --- a/src/coreComponents/physicsSolvers/fluidFlow/CompositionalMultiphaseStatisticsAggregator.cpp +++ b/src/coreComponents/physicsSolvers/fluidFlow/CompositionalMultiphaseStatisticsAggregator.cpp @@ -13,62 +13,6 @@ * ------------------------------------------------------------------------------------------------------------ */ -/** - * @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" @@ -88,10 +32,11 @@ namespace compositionalMultiphaseStatistics RegionStatistics::RegionStatistics( string const & name, dataRepository::Group * const parent, - bool statsOutputEnabled, + string_view setName, + bool dataOutputEnabled, integer const numPhases, integer const numComponents ): - RegionStatisticsBase( name, parent, statsOutputEnabled ), + RegionStatisticsBase( name, parent, setName, dataOutputEnabled ), m_phaseDynamicPoreVolume( numPhases ), m_phaseMass( numPhases ), m_trappedPhaseMass( numPhases ), @@ -105,15 +50,18 @@ RegionStatistics::RegionStatistics( string const & name, CFLStatistics::CFLStatistics( string const & name, dataRepository::Group * const parent, - bool const statsOutputEnabled ): - RegionStatisticsBase( name, parent, statsOutputEnabled ) + bool const dataOutputEnabled ): + RegionStatisticsBase( name, + parent, + "", // for now, CFL numbers are not for given sets + dataOutputEnabled ) { // TODO : registerWrappers to store results in HDF5 (but need repairing of 1D HDF5 outputs) } StatsAggregator::StatsAggregator( DataContext const & ownerDataContext, - bool const statsOutputEnabled ): - Base( ownerDataContext, statsOutputEnabled ), + bool const dataOutputEnabled ): + Base( ownerDataContext, dataOutputEnabled ), m_params() {} @@ -127,20 +75,23 @@ void StatsAggregator::initStatisticsAggregation( dataRepository::Group & meshBod Base::initStatisticsAggregation( meshBodies, solver ); } -void StatsAggregator::enableRegionStatisticsAggregation() +void StatsAggregator::enableRegionStatisticsAggregation( string_array const & setNames ) { auto const registerStats = [=] ( Group & parent, - string const & targetName ) -> RegionStatistics & + string const & targetName, + string_view setName, + bool const dataOutputEnabled ) -> RegionStatistics & { return parent.registerGroup( targetName, std::make_unique< RegionStatistics >( targetName, &parent, - m_statsOutputEnabled, + setName, + dataOutputEnabled, m_numPhases, m_numComponents ) ); }; - Base::enableRegionStatisticsAggregation( registerStats ); + Base::enableRegionStatisticsAggregation( registerStats, setNames ); } void StatsAggregator::enableCFLStatistics() @@ -157,7 +108,7 @@ void StatsAggregator::enableCFLStatistics() statisticsGroup.registerGroup< CFLStatistics >( cflStatsName, std::make_unique< CFLStatistics >( cflStatsName, &statisticsGroup, - m_statsOutputEnabled ) ); + m_dataOutputEnabled ) ); } m_cflStatsState.m_isEnabled = true; @@ -218,6 +169,8 @@ void StatsAggregator::initStats( RegionStatistics & stats, real64 const time ) c { stats.m_time = time; + stats.m_elemCount = 0; + stats.m_averagePressure = 0.0; stats.m_maxPressure = 0.0; stats.m_minPressure = LvArray::NumericLimits< real64 >::max; @@ -244,7 +197,8 @@ void StatsAggregator::initStats( RegionStatistics & stats, real64 const time ) c } void StatsAggregator::computeSubRegionRankStats( CellElementSubRegion & subRegion, - RegionStatistics & subRegionStats ) const + RegionStatistics & subRegionStats, + SetType const & targetSet ) const { arrayView1d< integer const > const elemGhostRank = subRegion.ghostRank(); arrayView1d< real64 const > const volume = subRegion.getElementVolume(); @@ -274,7 +228,7 @@ void StatsAggregator::computeSubRegionRankStats( CellElementSubRegion & subRegio isothermalCompositionalMultiphaseBaseKernels:: StatisticsKernel:: - launch< parallelDevicePolicy<> >( subRegion.size(), + launch< parallelDevicePolicy<> >( targetSet.toViewConst(), m_numComponents, m_numPhases, m_params.m_relpermThreshold, @@ -290,6 +244,7 @@ void StatsAggregator::computeSubRegionRankStats( CellElementSubRegion & subRegio phaseVolFrac, phaseTrappedVolFrac, phaseRelperm, + subRegionStats.m_elemCount, subRegionStats.m_minPressure, subRegionStats.m_averagePressure, subRegionStats.m_maxPressure, @@ -309,6 +264,8 @@ void StatsAggregator::computeSubRegionRankStats( CellElementSubRegion & subRegio void StatsAggregator::aggregateStats( RegionStatistics & stats, RegionStatistics const & other ) const { + stats.m_elemCount += other.m_elemCount; + 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 ); @@ -338,6 +295,8 @@ void StatsAggregator::aggregateStats( RegionStatistics & stats, void StatsAggregator::mpiAggregateStats( RegionStatistics & stats ) const { + stats.m_elemCount = MpiWrapper::sum( stats.m_elemCount ); + stats.m_averagePressure = MpiWrapper::sum( stats.m_averagePressure ); stats.m_minPressure = MpiWrapper::min( stats.m_minPressure ); stats.m_maxPressure = MpiWrapper::max( stats.m_maxPressure ); diff --git a/src/coreComponents/physicsSolvers/fluidFlow/CompositionalMultiphaseStatisticsAggregator.hpp b/src/coreComponents/physicsSolvers/fluidFlow/CompositionalMultiphaseStatisticsAggregator.hpp index 002fd40d53a..ce84e5db6d6 100644 --- a/src/coreComponents/physicsSolvers/fluidFlow/CompositionalMultiphaseStatisticsAggregator.hpp +++ b/src/coreComponents/physicsSolvers/fluidFlow/CompositionalMultiphaseStatisticsAggregator.hpp @@ -20,49 +20,56 @@ * Problem : ProblemManager * |-> domain : DomainPartition * |-> MeshBodies : Group - * |-> cartesianMesh : MeshBody + * |-> "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 + * | | |-> "Channel" : CellElementRegion + * | | | |-> "cb-0_0_0" : CellElementSubRegion * | | | | | * pressure : Wrapper< real64 array > * | | | | | * temperature : Wrapper< real64 array > * | | | | [...] (other fields) * | | | | - * | | | |-> cb-0_0_1 : CellElementSubRegion + * | | | |-> "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 + * | | |-> "Barrier" : CellElementRegion + * | | |-> "cb-1_0_0" : CellElementSubRegion + * | | |-> "cb-1_0_1" : CellElementSubRegion * | | [...] (other sub-regions) * | | + * | |-> nodeManager : NodeManager + * | | |-> sets : Group + * | | | * all : Wrapper< index array > + * | | | * "myBox" : Wrapper< index array > + * | | [...] (other element sets) + * | | * | [...] (other element managers) * ____ | | * | | |-> statistics : Group (storage for all stats) - * | | |-> compFlowStats : Group (storage for this instance 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) + * | | | |-> setsStatistics : RegionStatistics (selected sets aggregate, mpi reduced) + * | | | |-> "..." : RegionStatistics (set aggregate, mpi reduced) + * | | | | | - all requested sets + * | | | | | - "all" if no set restriction or if requested + * | | | | | - "__compound" if a compoud set is needed to compute mesh-level stats + * | | | | | + * | | | | |-> "Channel" : RegionStatistics (region aggregate, mpi reduced) + * | | | | | |-> "cb-0_0_0" : RegionStatistics (sub-region compute read-back) + * stats | | | | | |-> "cb-0_0_1" : RegionStatistics (sub-region compute read-back) + * data -> | | | | | [...] (other sub-regions stats) + * | | | | | + * | | | | |-> "Barrier" : RegionStatistics (region aggregate, mpi reduced) + * | | | | |-> cb-1_0_0 : RegionStatistics (sub-region compute read-back) + * | | | | |-> cb-1_0_1 : RegionStatistics (sub-region compute read-back) + * | | | | [...] (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 set aggregates) * | | | * |___ | [...] (other stats storages) * | @@ -113,6 +120,9 @@ class RegionStatistics : public RegionStatisticsBase { public: + /// Element Count + globalIndex m_elemCount; + /// average region pressure (numerator value before postAggregateCompute()) real64 m_averagePressure; /// minimum region pressure @@ -170,6 +180,7 @@ class RegionStatistics : public RegionStatisticsBase */ RegionStatistics( string const & targetName, dataRepository::Group * const parent, + string_view setName, bool statsOutputEnabled, integer numPhases, integer numComponents ); @@ -236,7 +247,6 @@ class StatsAggregator : public StatsAggregatorBase< StatsAggregator > * @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 ); @@ -245,8 +255,10 @@ class StatsAggregator : public StatsAggregatorBase< StatsAggregator > * @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 setNames The list of mesh element sets to restrict the statistics to. + * If empty, the whole discretization is processed. */ - void enableRegionStatisticsAggregation(); + void enableRegionStatisticsAggregation( string_array const & setNames = string_array() ); /** * @brief Register the results structs & wrappers so they will be targeted by TimeHistory output @@ -284,7 +296,9 @@ class StatsAggregator : public StatsAggregatorBase< StatsAggregator > /// @cond DO_NOT_DOCUMENT void initStats( RegionStatistics & stats, real64 time ) const; - void computeSubRegionRankStats( CellElementSubRegion & subRegion, RegionStatistics & subRegionStats ) const; + void computeSubRegionRankStats( CellElementSubRegion & subRegion, + RegionStatistics & subRegionStats, + SetType const & targetSet ) const; void aggregateStats( RegionStatistics & stats, RegionStatistics const & other ) const; void mpiAggregateStats( RegionStatistics & stats ) const; void postAggregateStats( RegionStatistics & stats ); diff --git a/src/coreComponents/physicsSolvers/fluidFlow/CompositionalMultiphaseStatisticsTask.cpp b/src/coreComponents/physicsSolvers/fluidFlow/CompositionalMultiphaseStatisticsTask.cpp index be83abde139..70b33bbf81e 100644 --- a/src/coreComponents/physicsSolvers/fluidFlow/CompositionalMultiphaseStatisticsTask.cpp +++ b/src/coreComponents/physicsSolvers/fluidFlow/CompositionalMultiphaseStatisticsTask.cpp @@ -39,6 +39,16 @@ StatsTask::StatsTask( string const & name, Group * const parent ): m_computeCFLNumbers( 0 ), m_computeRegionStatistics( 1 ) { + registerWrapper( viewKeyStruct::setNamesString(), &m_setNames ). + setRTTypeName( rtTypes::CustomTypes::groupNameRefArray ). + setInputFlag( InputFlags::OPTIONAL ). + setSizedFromParent( 0 ). + setDescription( "Optional targeted mesh element set(s) for which the statistics will be restricted. " + "A set can be be defined by a 'Geometry' component, or correspond to an imported index set " + "of mesh elements in case of an external mesh. " + "If empty, all mesh regions will be processed. " + "Be aware that only the regions that are processed by the solver will be taken into account." ); + registerWrapper( viewKeyStruct::computeCFLNumbersString(), &m_computeCFLNumbers ). setApplyDefaultValue( 0 ). setInputFlag( InputFlags::OPTIONAL ). @@ -103,7 +113,7 @@ void StatsTask::registerDataOnMesh( Group & meshBodies ) } if( m_computeRegionStatistics ) - m_aggregator->enableRegionStatisticsAggregation(); + m_aggregator->enableRegionStatisticsAggregation( m_setNames ); // if we have to compute CFL numbers later, we need to register additional variables if( m_computeCFLNumbers ) @@ -143,7 +153,8 @@ void StatsTask::prepareLogTableLayouts( string_view meshName ) return; TableLayout const tableLayout = TableLayout() - .setTitle( GEOS_FMT( "{}: mesh {}", getName(), meshName ) ); + .setTitle( GEOS_FMT( "Statistics: {} / Discretization: {}", + getName(), meshName ) ); m_logFormatters.emplace( meshName, std::make_unique< TableTextFormatter >( tableLayout ) ); } @@ -177,7 +188,8 @@ void StatsTask::prepareCsvTableLayouts( string_view meshName ) TableLayout tableLayout( { TableLayout::Column( GEOS_FMT( "Time [{}]", units::getSymbol( units::Unit::Time ))), - TableLayout::Column( "Region" ), // TODO : mention this change in PR description + TableLayout::Column( "Set" ), + TableLayout::Column( "Region" ), 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 ) ) ), @@ -222,12 +234,12 @@ bool StatsTask::execute( real64 const time_n, m_aggregator->computeRegionsStatistics( statsTime ); - m_aggregator->forRegionStatistics( [&] ( MeshLevel & mesh, RegionStatistics & meshRegionsStatistics ) + m_aggregator->forRegionStatistics( [&] ( MeshLevel & mesh, RegionStatistics & meshLevelStats ) { if( m_computeRegionStatistics ) { - outputLogStats( statsTime, mesh, meshRegionsStatistics ); - outputCsvStats( statsTime, mesh, meshRegionsStatistics ); + outputLogStats( statsTime, mesh, meshLevelStats ); + outputCsvStats( statsTime, mesh, meshLevelStats ); } } ); @@ -239,7 +251,7 @@ bool StatsTask::execute( real64 const time_n, void StatsTask::outputLogStats( real64 const statsTime, MeshLevel & mesh, - RegionStatistics & meshRegionsStatistics ) + RegionStatistics & meshLevelStats ) { if( MpiWrapper::commRank() > 0 || !isLogLevelActive< logInfo::Statistics >( this->getLogLevel() ) ) return; @@ -262,11 +274,12 @@ void StatsTask::outputLogStats( real64 const statsTime, tableData.addRow( "Statistics time", merge, merge, statsTime ); // lamda to apply for each region statistics - auto const outputRegionStats = [&] ( string_view targetName, RegionStatistics & stats ) + auto const outputRegionStats = [&] ( string_view title, + RegionStatistics & stats ) { tableData.addSeparator(); tableData.addRow( merge, merge, merge, "" ); - tableData.addRow( merge, merge, merge, targetName ); + tableData.addRow( merge, merge, merge, title ); tableData.addSeparator(); tableData.addRow( "statistics", "min", "average", "max" ); @@ -332,13 +345,41 @@ void StatsTask::outputLogStats( real64 const statsTime, stringutilities::join( stats.m_componentMass, '\n' ) ); }; - // apply the lambda for each region and, finally, the mesh summary - outputRegionStats( GEOS_FMT( "Discretization '{}'", mesh.getName() ), meshRegionsStatistics ); + // apply the output lambda for the mesh sets and regions + string const meshTitle = GEOS_FMT( "Discretization: '{}'", mesh.getName() ); + bool const logPerRegion = isLogLevelActive< logInfo::StatisticsPerRegion >( this->getLogLevel() ); + + if( m_aggregator->isTargetingMultipleSets() ) + { + string const allSetsTitle = m_aggregator->isRestrictedToSets() ? "Selected sets" : "All elements"; + outputRegionStats( GEOS_FMT( "{} / {} ({} elements)", + meshTitle, allSetsTitle, meshLevelStats.m_elemCount ), + meshLevelStats ); + } - m_aggregator->forRegionStatistics( mesh, meshRegionsStatistics, - [&] ( CellElementRegion & region, RegionStatistics & stats ) + m_aggregator->forRegionStatistics( mesh, meshLevelStats, false, + [&] ( StatsAggregator::MeshLevelSet meshSet, + RegionStatistics & meshSetStats ) { - outputRegionStats( GEOS_FMT( "Region '{}'", region.getName() ), stats ); + string const setTitle = m_aggregator->isRestrictedToSets() ? + GEOS_FMT( "Element set: '{}'", meshSet.setName ) : + "All elements"; + + outputRegionStats( GEOS_FMT( "{} / {} ({} elements)", + meshTitle, setTitle, meshSetStats.m_elemCount ), + meshSetStats ); + + if( logPerRegion ) + { + m_aggregator->forRegionStatistics( meshSet, meshSetStats, + [&] ( CellElementRegion & region, + RegionStatistics & setRegionStats ) + { + outputRegionStats( GEOS_FMT( "{} / {} / Region: '{}' ({} elements)", + meshTitle, setTitle, region.getName(), setRegionStats.m_elemCount ), + setRegionStats ); + } ); + } } ); // output to log @@ -347,7 +388,7 @@ void StatsTask::outputLogStats( real64 const statsTime, void StatsTask::outputCsvStats( real64 statsTime, MeshLevel & mesh, - RegionStatistics & meshRegionsStatistics ) + RegionStatistics & meshLevelStats ) { if( MpiWrapper::commRank() > 0 || m_writeCSV == 0 ) return; @@ -363,7 +404,9 @@ void StatsTask::outputCsvStats( real64 statsTime, row.reserve( formatter.getLayout().getTotalLowermostColumnCount() ); // lamda to apply for each region statistics - auto const outputRegionStats = [&] ( string_view targetName, RegionStatistics & stats ) + auto const outputRegionStats = [&] ( string_view targetName, + string_view setName, + RegionStatistics & stats ) { auto addPhaseValues = []( auto & list, auto const & values ) @@ -375,6 +418,7 @@ void StatsTask::outputCsvStats( real64 statsTime, row.clear(); row.insert( row.begin(), { std::to_string( statsTime ), + string( setName ), string( targetName ), std::to_string( stats.m_minPressure ), std::to_string( stats.m_averagePressure ), @@ -398,12 +442,23 @@ void StatsTask::outputCsvStats( real64 statsTime, }; // apply the lambda for each region and, finally, the mesh summary - m_aggregator->forRegionStatistics( mesh, meshRegionsStatistics, - [&] ( CellElementRegion & region, RegionStatistics & stats ) + string_view allSetsTitle = "Selected sets"; + + outputRegionStats( mesh.getName(), allSetsTitle, meshLevelStats ); + + m_aggregator->forRegionStatistics( mesh, meshLevelStats, false, + [&] ( StatsAggregator::MeshLevelSet meshSet, + RegionStatistics & meshSetStats ) { - outputRegionStats( region.getName(), stats ); + outputRegionStats( meshSet.mesh.getName(), meshSet.setName, meshSetStats ); + + m_aggregator->forRegionStatistics( meshSet, meshSetStats, + [&] ( CellElementRegion & region, + RegionStatistics & setRegionStats ) + { + outputRegionStats( region.getName(), setRegionStats.getSetName(), setRegionStats ); + } ); } ); - outputRegionStats( mesh.getName(), meshRegionsStatistics ); // append to csv file std::ofstream outputFile( getCsvFileName( mesh.getName() ), std::ios_base::app ); diff --git a/src/coreComponents/physicsSolvers/fluidFlow/CompositionalMultiphaseStatisticsTask.hpp b/src/coreComponents/physicsSolvers/fluidFlow/CompositionalMultiphaseStatisticsTask.hpp index 6f29bbe844f..87a4cf32161 100644 --- a/src/coreComponents/physicsSolvers/fluidFlow/CompositionalMultiphaseStatisticsTask.hpp +++ b/src/coreComponents/physicsSolvers/fluidFlow/CompositionalMultiphaseStatisticsTask.hpp @@ -84,6 +84,8 @@ class StatsTask : public FieldStatisticsBase< CompositionalMultiphaseBase > */ struct viewKeyStruct { + /// String for optionnal targeted element set(s) + constexpr static char const * setNamesString() { return "setNames"; } /// String for the flag deciding the computation of the CFL numbers constexpr static char const * computeCFLNumbersString() { return "computeCFLNumbers"; } /// String for the flag deciding the computation of the region statistics @@ -106,11 +108,11 @@ class StatsTask : public FieldStatisticsBase< CompositionalMultiphaseBase > void outputLogStats( real64 statsTime, MeshLevel & mesh, - RegionStatistics & meshRegionsStatistics ); + RegionStatistics & meshSetsStatistics ); void outputCsvStats( real64 statsTime, MeshLevel & mesh, - RegionStatistics & meshRegionsStatistics ); + RegionStatistics & meshSetsStatistics ); /// For each discretization (MeshLevel name), table formatter for log output. stdMap< string, std::unique_ptr< TableTextFormatter > > m_logFormatters; @@ -121,6 +123,9 @@ class StatsTask : public FieldStatisticsBase< CompositionalMultiphaseBase > // mesh statistics aggregator std::unique_ptr< StatsAggregator > m_aggregator; + /// Optionnal targeted element set(s) + string_array m_setNames; + /// Flag to decide whether CFL numbers are computed or not integer m_computeCFLNumbers; diff --git a/src/coreComponents/physicsSolvers/fluidFlow/SinglePhaseStatisticsAggregator.cpp b/src/coreComponents/physicsSolvers/fluidFlow/SinglePhaseStatisticsAggregator.cpp index 3797dfbfd5a..d75c477be8b 100644 --- a/src/coreComponents/physicsSolvers/fluidFlow/SinglePhaseStatisticsAggregator.cpp +++ b/src/coreComponents/physicsSolvers/fluidFlow/SinglePhaseStatisticsAggregator.cpp @@ -41,33 +41,39 @@ using namespace dataRepository; RegionStatistics::RegionStatistics( string const & name, dataRepository::Group * const parent, - bool const statsOutputEnabled ): - RegionStatisticsBase( name, parent, statsOutputEnabled ) + string_view setName, + bool const dataOutputEnabled ): + RegionStatisticsBase( name, parent, setName, dataOutputEnabled ) {} StatsAggregator::StatsAggregator( DataContext const & ownerDataContext, - bool const statsOutputEnabled ): - Base( ownerDataContext, statsOutputEnabled ) + bool const dataOutputEnabled ): + Base( ownerDataContext, dataOutputEnabled ) {} -void StatsAggregator::enableRegionStatisticsAggregation() +void StatsAggregator::enableRegionStatisticsAggregation( string_array const & setNames ) { auto const registerStats = [=] ( Group & parent, - string const & targetName ) -> RegionStatistics & + string const & targetName, + string_view setName, + bool const dataOutputEnabled ) -> RegionStatistics & { return parent.registerGroup( targetName, std::make_unique< RegionStatistics >( targetName, &parent, - m_statsOutputEnabled ) ); + setName, + dataOutputEnabled ) ); }; - Base::enableRegionStatisticsAggregation( registerStats ); + Base::enableRegionStatisticsAggregation( registerStats, setNames ); } void StatsAggregator::initStats( RegionStatistics & stats, real64 const time ) const { stats.m_time = time; + stats.m_elemCount = 0; + stats.m_averagePressure = 0.0; stats.m_maxPressure = 0.0; stats.m_minPressure = LvArray::NumericLimits< real64 >::max; @@ -86,7 +92,8 @@ void StatsAggregator::initStats( RegionStatistics & stats, real64 const time ) c } void StatsAggregator::computeSubRegionRankStats( CellElementSubRegion & subRegion, - RegionStatistics & subRegionStats ) const + RegionStatistics & subRegionStats, + SetType const & targetSet ) const { static constexpr string_view solidNamesVK = SinglePhaseBase::viewKeyStruct::solidNamesString(); static constexpr string_view fluidNamesVK = FlowSolverBase::viewKeyStruct::fluidNamesString(); @@ -108,7 +115,7 @@ void StatsAggregator::computeSubRegionRankStats( CellElementSubRegion & subRegio SingleFluidBase const & fluid = constitutiveModels.getGroup< SingleFluidBase >( fluidName ); arrayView2d< real64 const, constitutive::singlefluid::USD_FLUID > const densities = fluid.density(); - singlePhaseBaseKernels::StatisticsKernel::launch( subRegion.size(), + singlePhaseBaseKernels::StatisticsKernel::launch( targetSet.toViewConst(), elemGhostRank, volume, pres, @@ -117,6 +124,7 @@ void StatsAggregator::computeSubRegionRankStats( CellElementSubRegion & subRegio refPorosity, porosity, densities, + subRegionStats.m_elemCount, subRegionStats.m_minPressure, subRegionStats.m_averagePressure, subRegionStats.m_maxPressure, @@ -133,6 +141,8 @@ void StatsAggregator::computeSubRegionRankStats( CellElementSubRegion & subRegio void StatsAggregator::aggregateStats( RegionStatistics & stats, RegionStatistics const & other ) const { + stats.m_elemCount += other.m_elemCount; + 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 ); @@ -152,6 +162,8 @@ void StatsAggregator::aggregateStats( RegionStatistics & stats, void StatsAggregator::mpiAggregateStats( RegionStatistics & stats ) const { + stats.m_elemCount = MpiWrapper::sum( stats.m_elemCount ); + stats.m_averagePressure = MpiWrapper::sum( stats.m_averagePressure ); stats.m_minPressure = MpiWrapper::min( stats.m_minPressure ); stats.m_maxPressure = MpiWrapper::max( stats.m_maxPressure ); diff --git a/src/coreComponents/physicsSolvers/fluidFlow/SinglePhaseStatisticsAggregator.hpp b/src/coreComponents/physicsSolvers/fluidFlow/SinglePhaseStatisticsAggregator.hpp index ed56cc705c3..7c8e1a120f6 100644 --- a/src/coreComponents/physicsSolvers/fluidFlow/SinglePhaseStatisticsAggregator.hpp +++ b/src/coreComponents/physicsSolvers/fluidFlow/SinglePhaseStatisticsAggregator.hpp @@ -20,48 +20,67 @@ * Problem : ProblemManager * |-> domain : DomainPartition * |-> MeshBodies : Group - * |-> cartesianMesh : MeshBody + * |-> "cartesianMesh" : MeshBody * |-> meshLevels : Group * |-> Level0 : MeshLevel + * | |-> ElementRegions : ElementRegionManager + * | | |-> elementRegionsGroup : Group + * | | |-> "Channel" : CellElementRegion + * | | | |-> elementSubRegions : Group + * | | | |-> "cb-0_0_0" : CellElementSubRegion + * | | | | | * pressure : Wrapper< real64 array > + * | | | | | * temperature : Wrapper< real64 array > + * | | | | |-> sets : Group + * | | | | | | * all : Wrapper< index array > + * | | | | | | * myBox : Wrapper< index array > + * | | | | | [...] (other element sets) + * | | | | | + * | | | | [...] (other fields) + * | | | | + * | | | |-> "cb-0_0_1" : CellElementSubRegion + * | | | | | * pressure : Wrapper< real64 array > + * | | | | | * temperature : Wrapper< real64 array > + * | | | | |-> sets : Group + * | | | | | | * all : Wrapper< index array > + * | | | | | | * myBox : Wrapper< index array > + * | | | | | [...] (other element sets) + * | | | | | + * | | | | [...] (other fields) + * | | | | + * | | | [...] (other sub-regions) + * | | | + * | | |-> "Barrier" : CellElementRegion + * | | |-> "cb-1_0_0" : CellElementSubRegion + * | | |-> "cb-1_0_1" : CellElementSubRegion + * | | [...] (other sub-regions) + * | | * | |-> nodeManager : NodeManager * | | |-> sets : Group * | | | * all : Wrapper< index array > - * | | | * xneg : Wrapper< index array > + * | | | * "myBox" : 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) + * | | |-> "flowStats" : Group (storage for this instance stats) + * | | | |-> setsStatistics : RegionStatistics (selected sets aggregate, mpi reduced) + * | | | |-> "..." : RegionStatistics (set aggregate, mpi reduced) + * | | | | | - all requested sets + * | | | | | - "all" if no set restriction or if requested + * | | | | | - "__compound" if a compoud set is needed to compute mesh-level stats + * | | | | | + * | | | | |-> "Channel" : RegionStatistics (region aggregate, mpi reduced) + * | | | | | |-> "cb-0_0_0" : RegionStatistics (sub-region compute read-back) + * stats | | | | | |-> "cb-0_0_1" : RegionStatistics (sub-region compute read-back) + * data -> | | | | | [...] (other sub-regions stats) + * | | | | | + * | | | | |-> "Barrier" : RegionStatistics (region aggregate, mpi reduced) + * | | | | |-> cb-1_0_0 : RegionStatistics (sub-region compute read-back) + * | | | | |-> cb-1_0_1 : RegionStatistics (sub-region compute read-back) + * | | | | [...] (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 set aggregates) * | | | * |___ | [...] (other stats storages) * | @@ -103,8 +122,8 @@ class RegionStatistics : public RegionStatisticsBase { public: - /// Time of statistics computation - real64 m_time; + /// Element Count + globalIndex m_elemCount; /// average region pressure (numerator value before postAggregateCompute()) real64 m_averagePressure; @@ -146,6 +165,7 @@ class RegionStatistics : public RegionStatisticsBase */ RegionStatistics( string const & targetName, dataRepository::Group * const parent, + string_view setName, bool statsOutputEnabled ); RegionStatistics( RegionStatistics && ) = default; @@ -176,14 +196,18 @@ class StatsAggregator : public StatsAggregatorBase< StatsAggregator > * @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 setNames The list of mesh element sets to restrict the statistics to. + * If empty, the whole discretization is processed. */ - void enableRegionStatisticsAggregation(); + void enableRegionStatisticsAggregation( string_array const & setNames = string_array() ); // template implementations /// @cond DO_NOT_DOCUMENT void initStats( RegionStatistics & stats, real64 time ) const; - void computeSubRegionRankStats( CellElementSubRegion & subRegion, RegionStatistics & subRegionStats ) const; + void computeSubRegionRankStats( CellElementSubRegion & subRegion, + RegionStatistics & subRegionStats, + SetType const & targetSet ) const; void aggregateStats( RegionStatistics & stats, RegionStatistics const & other ) const; void mpiAggregateStats( RegionStatistics & stats ) const; void postAggregateStats( RegionStatistics & stats ); diff --git a/src/coreComponents/physicsSolvers/fluidFlow/SinglePhaseStatisticsTask.cpp b/src/coreComponents/physicsSolvers/fluidFlow/SinglePhaseStatisticsTask.cpp index 8166e5c10e5..e3b70d8350b 100644 --- a/src/coreComponents/physicsSolvers/fluidFlow/SinglePhaseStatisticsTask.cpp +++ b/src/coreComponents/physicsSolvers/fluidFlow/SinglePhaseStatisticsTask.cpp @@ -33,6 +33,16 @@ namespace singlePhaseStatistics StatsTask::StatsTask( const string & name, Group * const parent ): Base( name, parent ) { + registerWrapper( viewKeyStruct::setNamesString(), &m_setNames ). + setRTTypeName( rtTypes::CustomTypes::groupNameRefArray ). + setInputFlag( InputFlags::OPTIONAL ). + setSizedFromParent( 0 ). + setDescription( "Optional targeted mesh element set(s) for which the statistics will be restricted. " + "A set can be be defined by a 'Geometry' component, or correspond to an imported index set " + "of mesh elements in case of an external mesh. " + "If empty, all mesh regions will be processed. " + "Be aware that only the regions that are processed by the solver will be taken into account." ); + addLogLevel< logInfo::Statistics >(); } @@ -73,7 +83,7 @@ void StatsTask::registerDataOnMesh( Group & meshBodies ) catalogName(), getDataContext() ) ); } - m_aggregator->enableRegionStatisticsAggregation(); + m_aggregator->enableRegionStatisticsAggregation( m_setNames ); m_aggregator->forRegionStatistics( [&] ( MeshLevel & mesh, RegionStatistics & ) { @@ -89,7 +99,8 @@ void StatsTask::prepareLogTableLayouts( string_view meshName ) return; TableLayout const tableLayout = TableLayout() - .setTitle( GEOS_FMT( "{}: mesh {}", getName(), meshName ) ); + .setTitle( GEOS_FMT( "Statistics: {} / Discretization: {}", + getName(), meshName ) ); m_logFormatters.emplace( meshName, std::make_unique< TableTextFormatter >( tableLayout ) ); } @@ -104,7 +115,8 @@ void StatsTask::prepareCsvTableLayouts( string_view meshName ) TableLayout tableLayout( { TableLayout::Column( GEOS_FMT( "Time [{}]", units::getSymbol( units::Unit::Time ))), - TableLayout::Column( "Region" ), // TODO : mention this change in PR description + TableLayout::Column( "Set" ), + TableLayout::Column( "Region" ), 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 ) ) ), @@ -143,10 +155,10 @@ bool StatsTask::execute( real64 const time_n, m_aggregator->computeRegionsStatistics( statsTime ); - m_aggregator->forRegionStatistics( [&] ( MeshLevel & mesh, RegionStatistics & meshRegionsStatistics ) + m_aggregator->forRegionStatistics( [&] ( MeshLevel & mesh, RegionStatistics & meshLevelStats ) { - outputLogStats( statsTime, mesh, meshRegionsStatistics ); - outputCsvStats( statsTime, mesh, meshRegionsStatistics ); + outputLogStats( statsTime, mesh, meshLevelStats ); + outputCsvStats( statsTime, mesh, meshLevelStats ); } ); return false; @@ -154,7 +166,7 @@ bool StatsTask::execute( real64 const time_n, void StatsTask::outputLogStats( real64 const statsTime, MeshLevel & mesh, - RegionStatistics & meshRegionsStatistics ) + RegionStatistics & meshLevelStats ) { if( MpiWrapper::commRank() > 0 || !isLogLevelActive< logInfo::Statistics >( this->getLogLevel() ) ) return; @@ -177,11 +189,12 @@ void StatsTask::outputLogStats( real64 const statsTime, tableData.addRow( "Statistics time", merge, merge, statsTime ); // lamda to apply for each region statistics - auto const outputRegionStats = [&] ( string_view targetName, RegionStatistics & stats ) + auto const outputRegionStats = [&] ( string_view title, + RegionStatistics & stats ) { tableData.addSeparator(); tableData.addRow( merge, merge, merge, "" ); - tableData.addRow( merge, merge, merge, targetName ); + tableData.addRow( merge, merge, merge, title ); tableData.addSeparator(); tableData.addRow( "statistics", "min", "average", "max" ); @@ -203,13 +216,41 @@ void StatsTask::outputLogStats( real64 const statsTime, "all", CellType::MergeNext, stats.m_totalMass ); }; - // apply the output lambda for the mesh then each regions - outputRegionStats( GEOS_FMT( "Discretization '{}'", mesh.getName() ), meshRegionsStatistics ); + // apply the output lambda for the mesh sets and regions + string const meshTitle = GEOS_FMT( "Discretization: '{}'", mesh.getName() ); + bool const logPerRegion = isLogLevelActive< logInfo::StatisticsPerRegion >( this->getLogLevel() ); + + if( m_aggregator->isTargetingMultipleSets() ) + { + string const allSetsTitle = m_aggregator->isRestrictedToSets() ? "Selected sets" : "All elements"; + outputRegionStats( GEOS_FMT( "{} / {} ({} elements)", + meshTitle, allSetsTitle, meshLevelStats.m_elemCount ), + meshLevelStats ); + } - m_aggregator->forRegionStatistics( mesh, meshRegionsStatistics, - [&] ( CellElementRegion & region, RegionStatistics & stats ) + m_aggregator->forRegionStatistics( mesh, meshLevelStats, false, + [&] ( StatsAggregator::MeshLevelSet meshSet, + RegionStatistics & meshSetStats ) { - outputRegionStats( GEOS_FMT( "Region '{}'", region.getName() ), stats ); + string const setTitle = m_aggregator->isRestrictedToSets() ? + GEOS_FMT( "Element set: '{}'", meshSet.setName ) : + "All elements"; + + outputRegionStats( GEOS_FMT( "{} / {} ({} elements)", + meshTitle, setTitle, meshSetStats.m_elemCount ), + meshSetStats ); + + if( logPerRegion ) + { + m_aggregator->forRegionStatistics( meshSet, meshSetStats, + [&] ( CellElementRegion & region, + RegionStatistics & setRegionStats ) + { + outputRegionStats( GEOS_FMT( "{} / {} / Region: '{}' ({} elements)", + meshTitle, setTitle, region.getName(), setRegionStats.m_elemCount ), + setRegionStats ); + } ); + } } ); // output to log @@ -218,7 +259,7 @@ void StatsTask::outputLogStats( real64 const statsTime, void StatsTask::outputCsvStats( real64 statsTime, MeshLevel & mesh, - RegionStatistics & meshRegionsStatistics ) + RegionStatistics & meshLevelStats ) { if( MpiWrapper::commRank() > 0 || m_writeCSV == 0 ) return; @@ -234,11 +275,14 @@ void StatsTask::outputCsvStats( real64 statsTime, row.reserve( formatter.getLayout().getTotalLowermostColumnCount() ); // lamda to apply for each region statistics - auto const outputRegionStats = [&] ( string_view targetName, RegionStatistics & stats ) + auto const outputRegionStats = [&] ( string_view targetName, + string_view setName, + RegionStatistics & stats ) { row.clear(); row.insert( row.begin(), { std::to_string( statsTime ), + string( setName ), string( targetName ), std::to_string( stats.m_minPressure ), std::to_string( stats.m_averagePressure ), @@ -255,13 +299,23 @@ void StatsTask::outputCsvStats( real64 statsTime, tableData.addRow( row ); }; - // apply the output lambda for the mesh then each regions - outputRegionStats( mesh.getName(), meshRegionsStatistics ); + // apply the output lambda for the mesh sets and regions + string_view allSetsTitle = "Selected sets"; + + outputRegionStats( mesh.getName(), allSetsTitle, meshLevelStats ); - m_aggregator->forRegionStatistics( mesh, meshRegionsStatistics, - [&] ( CellElementRegion & region, RegionStatistics & stats ) + m_aggregator->forRegionStatistics( mesh, meshLevelStats, false, + [&] ( StatsAggregator::MeshLevelSet meshSet, + RegionStatistics & meshSetStats ) { - outputRegionStats( region.getName(), stats ); + outputRegionStats( meshSet.mesh.getName(), meshSet.setName, meshSetStats ); + + m_aggregator->forRegionStatistics( meshSet, meshSetStats, + [&] ( CellElementRegion & region, + RegionStatistics & setRegionStats ) + { + outputRegionStats( region.getName(), setRegionStats.getSetName(), setRegionStats ); + } ); } ); // append to csv file diff --git a/src/coreComponents/physicsSolvers/fluidFlow/SinglePhaseStatisticsTask.hpp b/src/coreComponents/physicsSolvers/fluidFlow/SinglePhaseStatisticsTask.hpp index 2d83a6b5a3f..c42e45133c4 100644 --- a/src/coreComponents/physicsSolvers/fluidFlow/SinglePhaseStatisticsTask.hpp +++ b/src/coreComponents/physicsSolvers/fluidFlow/SinglePhaseStatisticsTask.hpp @@ -76,6 +76,15 @@ class StatsTask : public FieldStatisticsBase< SinglePhaseBase > using Base = FieldStatisticsBase< SinglePhaseBase >; + /** + * @struct viewKeyStruct holds char strings and viewKeys for fast lookup + */ + struct viewKeyStruct + { + /// String for optionnal targeted element set(s) + constexpr static char const * setNamesString() { return "setNames"; } + }; + void postInputInitialization() override; void registerDataOnMesh( Group & meshBodies ) override; @@ -88,11 +97,11 @@ class StatsTask : public FieldStatisticsBase< SinglePhaseBase > void outputLogStats( real64 statsTime, MeshLevel & mesh, - RegionStatistics & meshRegionsStatistics ); + RegionStatistics & meshSetsStatistics ); void outputCsvStats( real64 statsTime, MeshLevel & mesh, - RegionStatistics & meshRegionsStatistics ); + RegionStatistics & meshSetsStatistics ); /// For each discretization (MeshLevel name), table formatter for log output. stdMap< string, std::unique_ptr< TableTextFormatter > > m_logFormatters; @@ -103,6 +112,9 @@ class StatsTask : public FieldStatisticsBase< SinglePhaseBase > // mesh statistics aggregator std::unique_ptr< StatsAggregator > m_aggregator; + /// Optionnal targeted element set(s) + string_array m_setNames; + }; } /* namespace singlePhaseStatistics */ diff --git a/src/coreComponents/physicsSolvers/fluidFlow/kernels/compositional/StatisticsKernel.hpp b/src/coreComponents/physicsSolvers/fluidFlow/kernels/compositional/StatisticsKernel.hpp index 3e98c375613..65f53880255 100644 --- a/src/coreComponents/physicsSolvers/fluidFlow/kernels/compositional/StatisticsKernel.hpp +++ b/src/coreComponents/physicsSolvers/fluidFlow/kernels/compositional/StatisticsKernel.hpp @@ -51,7 +51,7 @@ struct StatisticsKernel template< typename POLICY > static void - launch( localIndex const size, + launch( SortedArrayView< localIndex const > const & targetSet, integer const numComps, integer const numPhases, real64 const relpermThreshold, @@ -67,6 +67,7 @@ struct StatisticsKernel arrayView2d< real64 const, compflow::USD_PHASE > const & phaseVolFrac, arrayView3d< real64 const, constitutive::relperm::USD_RELPERM > const & phaseTrappedVolFrac, arrayView3d< real64 const, constitutive::relperm::USD_RELPERM > const & phaseRelperm, + globalIndex & elemCount, real64 & minPres, real64 & avgPresNumerator, real64 & maxPres, @@ -82,6 +83,7 @@ struct StatisticsKernel arrayView1d< real64 > const & immobilePhaseMass, arrayView2d< real64 > const & dissolvedComponentMass ) { + RAJA::ReduceSum< parallelDeviceReduce, localIndex > subRegionElemCount( 0 ); RAJA::ReduceMin< parallelDeviceReduce, real64 > subRegionMinPres( LvArray::NumericLimits< real64 >::max ); RAJA::ReduceSum< parallelDeviceReduce, real64 > subRegionAvgPresNumerator( 0.0 ); RAJA::ReduceMax< parallelDeviceReduce, real64 > subRegionMaxPres( -LvArray::NumericLimits< real64 >::max ); @@ -96,45 +98,50 @@ struct StatisticsKernel // using an array of ReduceSum leads to a formal parameter overflow in CUDA. // As a workaround, we use a slice with RAJA::atomicAdd instead - forAll< parallelDevicePolicy<> >( size, [numComps, - numPhases, - relpermThreshold, - elemGhostRank, - volume, - refPorosity, - porosity, - pres, - deltaPres, - temp, - phaseDensity, - phaseVolFrac, - phaseTrappedVolFrac, - phaseRelperm, - phaseCompFraction, - subRegionMinPres, - subRegionAvgPresNumerator, - subRegionMaxPres, - subRegionMinDeltaPres, - subRegionMaxDeltaPres, - subRegionMinTemp, - subRegionAvgTempNumerator, - subRegionMaxTemp, - subRegionTotalUncompactedPoreVol, - phaseDynamicPoreVol, - phaseMass, - trappedPhaseMass, - immobilePhaseMass, - dissolvedComponentMass] GEOS_HOST_DEVICE ( localIndex const ei ) + forAll< parallelDevicePolicy<> >( targetSet.size(), + [targetSet, + numComps, + numPhases, + relpermThreshold, + elemGhostRank, + volume, + refPorosity, + porosity, + pres, + deltaPres, + temp, + phaseDensity, + phaseVolFrac, + phaseTrappedVolFrac, + phaseRelperm, + phaseCompFraction, + subRegionElemCount, + subRegionMinPres, + subRegionAvgPresNumerator, + subRegionMaxPres, + subRegionMinDeltaPres, + subRegionMaxDeltaPres, + subRegionMinTemp, + subRegionAvgTempNumerator, + subRegionMaxTemp, + subRegionTotalUncompactedPoreVol, + phaseDynamicPoreVol, + phaseMass, + trappedPhaseMass, + immobilePhaseMass, + dissolvedComponentMass] GEOS_HOST_DEVICE ( localIndex const setElemId ) { + localIndex ei = targetSet[setElemId]; + if( elemGhostRank[ei] >= 0 ) - { return; - } // To match our "reference", we have to use reference porosity here, not the actual porosity when we compute averages real64 const uncompactedPoreVol = volume[ei] * refPorosity[ei]; real64 const dynamicPoreVol = volume[ei] * porosity[ei][0]; - + + subRegionElemCount += 1; + subRegionMinPres.min( pres[ei] ); subRegionAvgPresNumerator += uncompactedPoreVol * pres[ei]; subRegionMaxPres.max( pres[ei] ); @@ -168,14 +175,18 @@ struct StatisticsKernel } ); + elemCount = globalIndex( subRegionElemCount.get() ); + minPres = subRegionMinPres.get(); avgPresNumerator = subRegionAvgPresNumerator.get(); maxPres = subRegionMaxPres.get(); minDeltaPres = subRegionMinDeltaPres.get(); maxDeltaPres = subRegionMaxDeltaPres.get(); + minTemp = subRegionMinTemp.get(); avgTempNumerator = subRegionAvgTempNumerator.get(); maxTemp = subRegionMaxTemp.get(); + totalUncompactedPoreVol = subRegionTotalUncompactedPoreVol.get(); // dummy loop to bring data back to the CPU diff --git a/src/coreComponents/physicsSolvers/fluidFlow/kernels/singlePhase/StatisticsKernel.hpp b/src/coreComponents/physicsSolvers/fluidFlow/kernels/singlePhase/StatisticsKernel.hpp index 96f8934a2ec..3cdca4b6e43 100644 --- a/src/coreComponents/physicsSolvers/fluidFlow/kernels/singlePhase/StatisticsKernel.hpp +++ b/src/coreComponents/physicsSolvers/fluidFlow/kernels/singlePhase/StatisticsKernel.hpp @@ -46,7 +46,7 @@ struct StatisticsKernel } static void - launch( localIndex const size, + launch( SortedArrayView< localIndex const > const & targetSet, arrayView1d< integer const > const & elemGhostRank, arrayView1d< real64 const > const & volume, arrayView1d< real64 const > const & pres, @@ -55,6 +55,7 @@ struct StatisticsKernel arrayView1d< real64 const > const & refPorosity, arrayView2d< real64 const > const & porosity, arrayView2d< real64 const, constitutive::singlefluid::USD_FLUID > const & density, + globalIndex & elemCount, real64 & minPres, real64 & avgPresNumerator, real64 & maxPres, @@ -67,6 +68,8 @@ struct StatisticsKernel real64 & totalPoreVol, real64 & totalMass ) { + RAJA::ReduceSum< parallelDeviceReduce, localIndex > subRegionElemCount( 0 ); + RAJA::ReduceMin< parallelDeviceReduce, real64 > subRegionMinPres( LvArray::NumericLimits< real64 >::max ); RAJA::ReduceSum< parallelDeviceReduce, real64 > subRegionAvgPresNumerator( 0.0 ); RAJA::ReduceMax< parallelDeviceReduce, real64 > subRegionMaxPres( -LvArray::NumericLimits< real64 >::max ); @@ -82,17 +85,20 @@ struct StatisticsKernel RAJA::ReduceSum< parallelDeviceReduce, real64 > subRegionTotalPoreVol( 0.0 ); RAJA::ReduceSum< parallelDeviceReduce, real64 > subRegionTotalMass( 0.0 ); - forAll< parallelDevicePolicy<> >( size, [=] GEOS_HOST_DEVICE ( localIndex const ei ) + forAll< parallelDevicePolicy<> >( targetSet.size(), + [=] GEOS_HOST_DEVICE ( localIndex const setElemId ) { + localIndex ei = targetSet[setElemId]; + if( elemGhostRank[ei] >= 0 ) - { return; - } // To match our "reference", we have to use reference porosity here, not the actual porosity when we compute averages real64 const uncompactedPoreVol = volume[ei] * refPorosity[ei]; real64 const dynamicPoreVol = volume[ei] * porosity[ei][0]; + subRegionElemCount += 1; + subRegionMinPres.min( pres[ei] ); subRegionAvgPresNumerator += uncompactedPoreVol * pres[ei]; subRegionMaxPres.max( pres[ei] ); @@ -109,6 +115,8 @@ struct StatisticsKernel subRegionTotalMass += dynamicPoreVol * density[ei][0]; } ); + elemCount = globalIndex( subRegionElemCount.get() ); + minPres = subRegionMinPres.get(); avgPresNumerator = subRegionAvgPresNumerator.get(); maxPres = subRegionMaxPres.get(); diff --git a/src/coreComponents/physicsSolvers/fluidFlow/unitTests/testFlowStatistics.cpp b/src/coreComponents/physicsSolvers/fluidFlow/unitTests/testFlowStatistics.cpp index e61551c70c3..ed50856e940 100644 --- a/src/coreComponents/physicsSolvers/fluidFlow/unitTests/testFlowStatistics.cpp +++ b/src/coreComponents/physicsSolvers/fluidFlow/unitTests/testFlowStatistics.cpp @@ -221,7 +221,7 @@ real64 getTotalFluidMass( ProblemManager & problem, .getMeshBody( 0 ) .getMeshLevel( solver.getDiscretizationName() ); auto const & statsAggregator = statsTask.getStatisticsAggregator(); - auto const & stats = statsAggregator.getRegionsStatistics( mesh ); + auto const & stats = statsAggregator.getMeshLevelStatistics( mesh ); return stats.m_totalMass; } diff --git a/src/coreComponents/schema/schema.xsd b/src/coreComponents/schema/schema.xsd index 728d7ec8eff..45ffe698e1c 100644 --- a/src/coreComponents/schema/schema.xsd +++ b/src/coreComponents/schema/schema.xsd @@ -6324,6 +6324,8 @@ Information output from lower logLevels is added with the desired log level + + @@ -6542,6 +6544,8 @@ Information output from lower logLevels is added with the desired log level 1 - Print statistics--> + +