From 4d7e205366597f1c3257e1711f372e1c7a41c4c8 Mon Sep 17 00:00:00 2001 From: dkachuma Date: Mon, 11 Mar 2024 12:21:03 -0500 Subject: [PATCH 1/7] First pass --- .../physicsSolvers/CMakeLists.txt | 2 + ...rmalCompositionalMultiphaseBaseKernels.hpp | 34 ++++ .../fluidFlow/RegionMultiphaseStatistics.cpp | 136 +++++++++++++++ .../fluidFlow/RegionMultiphaseStatistics.hpp | 157 ++++++++++++++++++ .../docs/RegionMultiphaseStatistics.rst | 15 ++ .../docs/RegionMultiphaseStatistics_other.rst | 9 + src/coreComponents/schema/docs/Tasks.rst | 1 + .../schema/docs/Tasks_other.rst | 1 + src/coreComponents/schema/schema.xsd | 21 +++ src/coreComponents/schema/schema.xsd.other | 2 + src/docs/sphinx/CompleteXMLSchema.rst | 14 ++ 11 files changed, 392 insertions(+) create mode 100644 src/coreComponents/physicsSolvers/fluidFlow/RegionMultiphaseStatistics.cpp create mode 100644 src/coreComponents/physicsSolvers/fluidFlow/RegionMultiphaseStatistics.hpp create mode 100644 src/coreComponents/schema/docs/RegionMultiphaseStatistics.rst create mode 100644 src/coreComponents/schema/docs/RegionMultiphaseStatistics_other.rst diff --git a/src/coreComponents/physicsSolvers/CMakeLists.txt b/src/coreComponents/physicsSolvers/CMakeLists.txt index c7d26a6bfd3..22c3c312b58 100644 --- a/src/coreComponents/physicsSolvers/CMakeLists.txt +++ b/src/coreComponents/physicsSolvers/CMakeLists.txt @@ -18,6 +18,7 @@ set( physicsSolvers_headers fluidFlow/CompositionalMultiphaseBase.hpp fluidFlow/CompositionalMultiphaseBaseFields.hpp fluidFlow/CompositionalMultiphaseStatistics.hpp + fluidFlow/RegionMultiphaseStatistics.hpp fluidFlow/IsothermalCompositionalMultiphaseBaseKernels.hpp fluidFlow/ThermalCompositionalMultiphaseBaseKernels.hpp fluidFlow/CompositionalMultiphaseFVM.hpp @@ -160,6 +161,7 @@ set( physicsSolvers_sources fluidFlow/CompositionalMultiphaseBase.cpp fluidFlow/CompositionalMultiphaseFVM.cpp fluidFlow/CompositionalMultiphaseStatistics.cpp + fluidFlow/RegionMultiphaseStatistics.cpp fluidFlow/IsothermalCompositionalMultiphaseFVMKernels.cpp fluidFlow/CompositionalMultiphaseHybridFVM.cpp fluidFlow/CompositionalMultiphaseHybridFVMKernels.cpp diff --git a/src/coreComponents/physicsSolvers/fluidFlow/IsothermalCompositionalMultiphaseBaseKernels.hpp b/src/coreComponents/physicsSolvers/fluidFlow/IsothermalCompositionalMultiphaseBaseKernels.hpp index 8adafb2945a..87f6a48c75d 100644 --- a/src/coreComponents/physicsSolvers/fluidFlow/IsothermalCompositionalMultiphaseBaseKernels.hpp +++ b/src/coreComponents/physicsSolvers/fluidFlow/IsothermalCompositionalMultiphaseBaseKernels.hpp @@ -2085,6 +2085,40 @@ struct StatisticsKernel } }; +/******************************** RegionStatisticsKernel ********************************/ + +struct RegionStatisticsKernel +{ + template< typename POLICY > + static void + launch( localIndex const size, + arrayView1d< integer const > const & elementGhostRank, + arrayView1d< real64 const > const & elementVolume, + arrayView1d< real64 const > const & pressure, + arrayView1d< real64 const > const & temperature, + arrayView1d< real64 const > const & referencePorosity, + arrayView2d< real64 const > const & porosity, + arrayView3d< real64 const, multifluid::USD_PHASE > const & phaseDensity, + arrayView4d< real64 const, multifluid::USD_PHASE_COMP > const & phaseComponentFraction, + arrayView2d< real64 const, compflow::USD_PHASE > const & phaseVolumeFraction, + arrayView3d< real64 const, relperm::USD_RELPERM > const & phaseTrappedVolumeFraction, + arrayView3d< real64 const, relperm::USD_RELPERM > const & phaseRelativePermeability ) + { +GEOS_UNUSED_VAR(elementGhostRank); +GEOS_UNUSED_VAR(elementVolume); +GEOS_UNUSED_VAR(pressure); +GEOS_UNUSED_VAR(temperature); +GEOS_UNUSED_VAR(referencePorosity); +GEOS_UNUSED_VAR(porosity); +GEOS_UNUSED_VAR(phaseDensity); +GEOS_UNUSED_VAR(phaseComponentFraction); +GEOS_UNUSED_VAR(phaseVolumeFraction); +GEOS_UNUSED_VAR(phaseTrappedVolumeFraction); +GEOS_UNUSED_VAR(phaseRelativePermeability); + } +}; + + /******************************** HydrostaticPressureKernel ********************************/ struct HydrostaticPressureKernel diff --git a/src/coreComponents/physicsSolvers/fluidFlow/RegionMultiphaseStatistics.cpp b/src/coreComponents/physicsSolvers/fluidFlow/RegionMultiphaseStatistics.cpp new file mode 100644 index 00000000000..6774f32beaa --- /dev/null +++ b/src/coreComponents/physicsSolvers/fluidFlow/RegionMultiphaseStatistics.cpp @@ -0,0 +1,136 @@ +/* + * ------------------------------------------------------------------------------------------------------------ + * SPDX-License-Identifier: LGPL-2.1-only + * + * Copyright (c) 2018-2020 Lawrence Livermore National Security LLC + * Copyright (c) 2018-2020 The Board of Trustees of the Leland Stanford Junior University + * Copyright (c) 2018-2020 TotalEnergies + * Copyright (c) 2019- GEOSX Contributors + * All rights reserved + * + * See top level LICENSE, COPYRIGHT, CONTRIBUTORS, NOTICE, and ACKNOWLEDGEMENTS files for details. + * ------------------------------------------------------------------------------------------------------------ + */ + +/** + * @file RegionMultiphaseStatistics.cpp + */ + +#include "RegionMultiphaseStatistics.hpp" +//#include "constitutive/fluid/multifluid/MultiFluidBase.hpp" +//#include "constitutive/relativePermeability/RelativePermeabilityBase.hpp" +//#include "constitutive/solid/CoupledSolidBase.hpp" +//#include "finiteVolume/FiniteVolumeManager.hpp" +//#include "finiteVolume/FluxApproximationBase.hpp" +//#include "mainInterface/ProblemManager.hpp" +//#include "physicsSolvers/PhysicsSolverManager.hpp" +#include "physicsSolvers/fluidFlow/CompositionalMultiphaseBase.hpp" +#include "mesh/MeshFields.hpp" +//#include "physicsSolvers/fluidFlow/CompositionalMultiphaseBaseFields.hpp" +//#include "physicsSolvers/fluidFlow/CompositionalMultiphaseHybridFVM.hpp" +//#include "physicsSolvers/fluidFlow/FlowSolverBaseFields.hpp" +//#include "physicsSolvers/fluidFlow/IsothermalCompositionalMultiphaseBaseKernels.hpp" +//#include "physicsSolvers/fluidFlow/IsothermalCompositionalMultiphaseFVMKernels.hpp" + +namespace geos +{ + +using namespace constitutive; +using namespace dataRepository; + +RegionMultiphaseStatistics::RegionMultiphaseStatistics( const string & name, + Group * const parent ): + Base( name, parent ) +{ + registerWrapper( viewKeyStruct::regionNamesString(), &m_regionNames ). + setInputFlag( InputFlags::REQUIRED ). + setDescription( "The list of names of the regions" ); + + registerWrapper( viewKeyStruct::regionIdentifiersString(), &m_regionIdentifiers ). + setInputFlag( InputFlags::REQUIRED ). + setDescription( "The list of integer identifiers of the regions" ); + + registerWrapper( viewKeyStruct::propertyNamesString(), &m_propertyNames ). + setInputFlag( InputFlags::OPTIONAL ). + setDescription( "The list of names of properties to report" ); +} + +void RegionMultiphaseStatistics::postProcessInput() +{ + Base::postProcessInput(); + + // Make sure we have the same number of region names and region indices + GEOS_ERROR_IF_NE(m_regionNames.size(), m_regionIdentifiers.size(), + GEOS_FMT( "The number of region names {} is different from the number of region identifiers {}.", + m_regionNames.size(), m_regionIdentifiers.size() ) ); + + // Convert the property names to integers + m_propertyNameTypes.clear(); + localIndex const propSize = m_propertyNames.size(); + if( propSize == 0 ) + { +m_propertyNameTypes.emplace_back(PropertyNameType::Pressure); +m_propertyNameTypes.emplace_back(PropertyNameType::PoreVolume); +m_propertyNameTypes.emplace_back(PropertyNameType::Saturation); + } + else + { + m_propertyNameTypes.resize( propSize ); + for( integer i = 0; i < propSize; ++i ) + { + m_propertyNameTypes[i] = EnumStrings< PropertyNameType >::fromString( m_propertyNames[i] ); + } + } +} +void RegionMultiphaseStatistics::registerDataOnMesh( Group & meshBodies ) +{ + Base::registerDataOnMesh( 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; + } + + string const fieldName = GEOS_FMT( "{}{}", getName(), viewKeyStruct::fieldNameString() ); + + m_solver->forDiscretizationOnMeshTargets( meshBodies, [&] ( string const &, + MeshLevel & mesh, + arrayView1d< string const > const & regionNames ) + { + ElementRegionManager & elemManager = mesh.getElemManager(); +elemManager.forElementSubRegions( regionNames, [&]( localIndex const, + ElementSubRegionBase & subRegion ) + { + subRegion.registerWrapper>( fieldName ). + setPlotLevel( PlotLevel::NOPLOT ). + setRestartFlags( RestartFlags::NO_WRITE ). + setDescription( fieldName ). + setRegisteringObjects( getName() ); + }); + }); +} + +bool RegionMultiphaseStatistics::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 ) +{ + GEOS_UNUSED_VAR( time_n ); + GEOS_UNUSED_VAR( dt ); + GEOS_UNUSED_VAR( domain ); + std::cout << "*** Executing " << getName() << "\n"; + return false; +} + +REGISTER_CATALOG_ENTRY( TaskBase, + RegionMultiphaseStatistics, + string const &, dataRepository::Group * const ) + +} /* namespace geos */ + diff --git a/src/coreComponents/physicsSolvers/fluidFlow/RegionMultiphaseStatistics.hpp b/src/coreComponents/physicsSolvers/fluidFlow/RegionMultiphaseStatistics.hpp new file mode 100644 index 00000000000..49e1bf72571 --- /dev/null +++ b/src/coreComponents/physicsSolvers/fluidFlow/RegionMultiphaseStatistics.hpp @@ -0,0 +1,157 @@ +/* + * ------------------------------------------------------------------------------------------------------------ + * SPDX-License-Identifier: LGPL-2.1-only + * + * Copyright (c) 2018-2020 Lawrence Livermore National Security LLC + * Copyright (c) 2018-2020 The Board of Trustees of the Leland Stanford Junior University + * Copyright (c) 2018-2020 TotalEnergies + * Copyright (c) 2019- GEOSX Contributors + * All rights reserved + * + * See top level LICENSE, COPYRIGHT, CONTRIBUTORS, NOTICE, and ACKNOWLEDGEMENTS files for details. + * ------------------------------------------------------------------------------------------------------------ + */ + +/** + * @file RegionMultiphaseStatistics.hpp + */ + +#ifndef GEOS_CORECOMPONENTS_PHYSICSSOLVERS_FLUIDFLOW_REGIONMULTIPHASESTATISTICS_HPP_ +#define GEOS_CORECOMPONENTS_PHYSICSSOLVERS_FLUIDFLOW_REGIONMULTIPHASESTATISTICS_HPP_ + +#include "physicsSolvers/FieldStatisticsBase.hpp" + +#include "codingUtilities/EnumStrings.hpp" + +namespace geos +{ + +class CompositionalMultiphaseBase; + +/** + * @class RegionMultiphaseStatistics + * + * Task class allowing for the computation of aggregate statistics in compositional multiphase simulations + */ +class RegionMultiphaseStatistics : public FieldStatisticsBase< CompositionalMultiphaseBase > +{ +public: + enum PropertyNameType : integer + { + Pressure, + PoreVolume, + Saturation + }; + +public: + + /** + * @brief Constructor for the statistics class + * @param[in] name the name of the task coming from the xml + * @param[in] parent the parent group of the task + */ + RegionMultiphaseStatistics( const string & name, + Group * const parent ); + + /// Accessor for the catalog name + static string catalogName() { return "RegionMultiphaseStatistics"; } + + /** + * @defgroup Tasks Interface Functions + * + * This function implements the interface defined by the abstract TaskBase class + */ + /**@{*/ + + virtual bool execute( real64 const time_n, + real64 const dt, + integer const cycleNumber, + integer const eventCounter, + real64 const eventProgress, + DomainPartition & domain ) override; + + /**@}*/ + +private: + + using Base = FieldStatisticsBase< CompositionalMultiphaseBase >; + + /** + * @struct viewKeyStruct holds char strings and viewKeys for fast lookup + */ + struct viewKeyStruct + { + /// String for the list of region names + constexpr static char const * regionNamesString() { return "regionNames"; } + /// String for the list of region indices + constexpr static char const * regionIdentifiersString() { return "regionIdentifiers"; } + /// String for the list of properties to report + constexpr static char const * propertyNamesString() { return "propertyNames"; } + /// String for the suffix of the field name + constexpr static char const * fieldNameString() { return "_region"; } + }; + + struct RegionStatistics + { + explicit RegionStatistics(integer regionCount); + /// 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; + }; + + void postProcessInput() override; + + void registerDataOnMesh( Group & meshBodies ) override; + + // The names of regions + string_array m_regionNames; + + // The ids of regions + array1d m_regionIdentifiers; + + // The names of properties to output + string_array m_propertyNames; + + // The names converted into integers + array1d< PropertyNameType > m_propertyNameTypes; +}; + +ENUM_STRINGS( RegionMultiphaseStatistics::PropertyNameType, + "pressure", + "poreVolume", + "saturation" ); + +} /* namespace geos */ + +#endif /* GEOS_CORECOMPONENTS_PHYSICSSOLVERS_FLUIDFLOW_REGIONMULTIPHASESTATISTICS_HPP_ */ diff --git a/src/coreComponents/schema/docs/RegionMultiphaseStatistics.rst b/src/coreComponents/schema/docs/RegionMultiphaseStatistics.rst new file mode 100644 index 00000000000..245e94393dc --- /dev/null +++ b/src/coreComponents/schema/docs/RegionMultiphaseStatistics.rst @@ -0,0 +1,15 @@ + + +================= ============= ======== ============================================== +Name Type Default Description +================= ============= ======== ============================================== +flowSolverName groupNameRef required Name of the flow solver +logLevel integer 0 Log level +name groupName required A name is required for any non-unique nodes +propertyNames string_array {} The list of names of properties to report +regionIdentifiers integer_array required The list of integer identifiers of the regions +regionNames string_array required The list of names of the regions +writeCSV integer 0 Write statistics into a CSV file +================= ============= ======== ============================================== + + diff --git a/src/coreComponents/schema/docs/RegionMultiphaseStatistics_other.rst b/src/coreComponents/schema/docs/RegionMultiphaseStatistics_other.rst new file mode 100644 index 00000000000..adf1c1b8aec --- /dev/null +++ b/src/coreComponents/schema/docs/RegionMultiphaseStatistics_other.rst @@ -0,0 +1,9 @@ + + +==== ==== ============================ +Name Type Description +==== ==== ============================ + (no documentation available) +==== ==== ============================ + + diff --git a/src/coreComponents/schema/docs/Tasks.rst b/src/coreComponents/schema/docs/Tasks.rst index 95301497c5c..2cbee4ce5c3 100644 --- a/src/coreComponents/schema/docs/Tasks.rst +++ b/src/coreComponents/schema/docs/Tasks.rst @@ -9,6 +9,7 @@ MultiphasePoromechanicsInitialization node :ref:`X PVTDriver node :ref:`XML_PVTDriver` PackCollection node :ref:`XML_PackCollection` ReactiveFluidDriver node :ref:`XML_ReactiveFluidDriver` +RegionMultiphaseStatistics node :ref:`XML_RegionMultiphaseStatistics` RelpermDriver node :ref:`XML_RelpermDriver` SinglePhasePoromechanicsInitialization node :ref:`XML_SinglePhasePoromechanicsInitialization` SinglePhaseReservoirPoromechanicsInitialization node :ref:`XML_SinglePhaseReservoirPoromechanicsInitialization` diff --git a/src/coreComponents/schema/docs/Tasks_other.rst b/src/coreComponents/schema/docs/Tasks_other.rst index dca37a26a99..742de7a8c31 100644 --- a/src/coreComponents/schema/docs/Tasks_other.rst +++ b/src/coreComponents/schema/docs/Tasks_other.rst @@ -9,6 +9,7 @@ MultiphasePoromechanicsInitialization node :ref:`DATASTRUC PVTDriver node :ref:`DATASTRUCTURE_PVTDriver` PackCollection node :ref:`DATASTRUCTURE_PackCollection` ReactiveFluidDriver node :ref:`DATASTRUCTURE_ReactiveFluidDriver` +RegionMultiphaseStatistics node :ref:`DATASTRUCTURE_RegionMultiphaseStatistics` RelpermDriver node :ref:`DATASTRUCTURE_RelpermDriver` SinglePhasePoromechanicsInitialization node :ref:`DATASTRUCTURE_SinglePhasePoromechanicsInitialization` SinglePhaseReservoirPoromechanicsInitialization node :ref:`DATASTRUCTURE_SinglePhaseReservoirPoromechanicsInitialization` diff --git a/src/coreComponents/schema/schema.xsd b/src/coreComponents/schema/schema.xsd index bb9ab1cbbdf..937ee56113d 100644 --- a/src/coreComponents/schema/schema.xsd +++ b/src/coreComponents/schema/schema.xsd @@ -483,6 +483,10 @@ + + + + @@ -3544,6 +3548,7 @@ Local - Add stabilization only to interiors of macro elements.--> + @@ -3647,6 +3652,22 @@ Local - Add stabilization only to interiors of macro elements.--> + + + + + + + + + + + + + + + + diff --git a/src/coreComponents/schema/schema.xsd.other b/src/coreComponents/schema/schema.xsd.other index c81efef7141..f3b6e20022d 100644 --- a/src/coreComponents/schema/schema.xsd.other +++ b/src/coreComponents/schema/schema.xsd.other @@ -1285,6 +1285,7 @@ + @@ -1300,6 +1301,7 @@ + diff --git a/src/docs/sphinx/CompleteXMLSchema.rst b/src/docs/sphinx/CompleteXMLSchema.rst index 73517d48697..a8b8449c494 100644 --- a/src/docs/sphinx/CompleteXMLSchema.rst +++ b/src/docs/sphinx/CompleteXMLSchema.rst @@ -997,6 +997,13 @@ Element: Rectangle .. include:: ../../coreComponents/schema/docs/Rectangle.rst +.. _XML_RegionMultiphaseStatistics: + +Element: RegionMultiphaseStatistics +=================================== +.. include:: ../../coreComponents/schema/docs/RegionMultiphaseStatistics.rst + + .. _XML_RelpermDriver: Element: RelpermDriver @@ -2422,6 +2429,13 @@ Datastructure: Rectangle .. include:: ../../coreComponents/schema/docs/Rectangle_other.rst +.. _DATASTRUCTURE_RegionMultiphaseStatistics: + +Datastructure: RegionMultiphaseStatistics +========================================= +.. include:: ../../coreComponents/schema/docs/RegionMultiphaseStatistics_other.rst + + .. _DATASTRUCTURE_RelpermDriver: Datastructure: RelpermDriver From b91f59c2953deec29cf94e0ac9cb136d697dca9e Mon Sep 17 00:00:00 2001 From: dkachuma Date: Mon, 11 Mar 2024 16:38:29 -0500 Subject: [PATCH 2/7] Move kernel to implementation file --- ...rmalCompositionalMultiphaseBaseKernels.hpp | 234 ++++++- .../fluidFlow/RegionMultiphaseStatistics.cpp | 575 +++++++++++++----- .../fluidFlow/RegionMultiphaseStatistics.hpp | 54 +- 3 files changed, 676 insertions(+), 187 deletions(-) diff --git a/src/coreComponents/physicsSolvers/fluidFlow/IsothermalCompositionalMultiphaseBaseKernels.hpp b/src/coreComponents/physicsSolvers/fluidFlow/IsothermalCompositionalMultiphaseBaseKernels.hpp index 87f6a48c75d..6e972ff7ecd 100644 --- a/src/coreComponents/physicsSolvers/fluidFlow/IsothermalCompositionalMultiphaseBaseKernels.hpp +++ b/src/coreComponents/physicsSolvers/fluidFlow/IsothermalCompositionalMultiphaseBaseKernels.hpp @@ -2089,9 +2089,123 @@ struct StatisticsKernel struct RegionStatisticsKernel { + struct RegionStatistics + { + /// Indices of minimum, maximum and total property values + static constexpr integer MINIMUM = 0; + static constexpr integer MAXIMUM = 1; + static constexpr integer TOTAL = 2; + static constexpr integer WEIGHTED_TOTAL = 3; + static constexpr integer TYPE_END = 4; + + /// Indices for property type + static constexpr integer PORE_VOLUME = 0; + static constexpr integer PRESSURE = 1; + static constexpr integer TEMPERATURE = 2; + static constexpr integer PHASE_PORE_VOLUME = 3; + static constexpr integer PHASE_MASS = 4; + static constexpr integer PHASE_COMP_MASS = 5; + static constexpr integer PROP_END = 6; + + RegionStatistics( integer regionCount, integer phaseCount, integer componentCount ) + : m_phaseCount( phaseCount ), + m_componentCount( componentCount ), + m_properties( TYPE_END, regionCount, getSize( phaseCount, componentCount )), + m_elementCounts( regionCount ) + { + LvArray::forValuesInSlice( m_properties.toSlice(), []( real64 & v ){ v = 0.0; } ); + LvArray::forValuesInSlice( m_elementCounts.toSlice(), []( localIndex & i ){ i = 0; } ); + + static constexpr real64 minValue = -LvArray::NumericLimits< real64 >::max; + static constexpr real64 maxValue = LvArray::NumericLimits< real64 >::max; + LvArray::forValuesInSlice( m_properties[MINIMUM], []( real64 & v ){ v = maxValue; } ); + LvArray::forValuesInSlice( m_properties[MAXIMUM], []( real64 & v ){ v = minValue; } ); + } + + real64 & getProperty( integer region, integer type, integer property, integer phase = 0, integer component = 0 ) + { + return m_properties( type, region, calculateIndex( property, phase, component )); + } + + real64 getProperty( integer region, integer type, integer property, integer phase = 0, integer component = 0 ) const + { + return m_properties( type, region, calculateIndex( property, phase, component )); + } + + static real64 * getProperty( arrayView3d< real64 > props, + integer phaseCount, + integer componentCount, + integer region, + integer type, + integer property, + integer phase = 0, + integer component = 0 ) + { + return &props( type, region, calculateIndex( phaseCount, componentCount, property, phase, component )); + } + + arrayView3d< real64 > toView() const + { + return m_properties.toView(); + } + + arrayView1d< localIndex > getElementCounts() const + { + return m_elementCounts.toView(); + } + +private: + static integer getSize( integer phaseCount, integer componentCount ) + { + // 1 pore volume, 1 pressure, 1 temperature, np phase volume, np phase mass, np*nc phase comp mass + return 3 + 2*phaseCount + phaseCount*componentCount; + } + + integer calculateIndex( integer property, integer phase, integer component ) const + { + return calculateIndex( m_phaseCount, m_componentCount, property, phase, component ); + } + + static integer calculateIndex( integer phaseCount, + integer componentCount, + integer property, + integer phase, + integer component ) + { + if( property == PORE_VOLUME || property == PRESSURE || property == TEMPERATURE ) + { + return property; + } + else if( property == PHASE_PORE_VOLUME ) + { + return PHASE_PORE_VOLUME + phase; + } + else if( property == PHASE_MASS ) + { + return PHASE_PORE_VOLUME + phaseCount; + } + else if( property == PHASE_COMP_MASS ) + { + return PHASE_PORE_VOLUME + 2*phaseCount + phase*componentCount + component; + } + return 0; + } + + integer const m_phaseCount{0}; + integer const m_componentCount{0}; + array3d< real64 > m_properties; + array1d< localIndex > m_elementCounts; + }; + template< typename POLICY > static void launch( localIndex const size, + integer const phaseCount, + integer const componentCount, + arrayView3d< real64 > const regionStatistics, + arrayView1d< localIndex > const & elementCounts, + arrayView1d< real64 const > const & regionIndices, + arrayView1d< real64 const > const & regionMarkers, arrayView1d< integer const > const & elementGhostRank, arrayView1d< real64 const > const & elementVolume, arrayView1d< real64 const > const & pressure, @@ -2104,17 +2218,115 @@ struct RegionStatisticsKernel arrayView3d< real64 const, relperm::USD_RELPERM > const & phaseTrappedVolumeFraction, arrayView3d< real64 const, relperm::USD_RELPERM > const & phaseRelativePermeability ) { -GEOS_UNUSED_VAR(elementGhostRank); -GEOS_UNUSED_VAR(elementVolume); -GEOS_UNUSED_VAR(pressure); -GEOS_UNUSED_VAR(temperature); -GEOS_UNUSED_VAR(referencePorosity); -GEOS_UNUSED_VAR(porosity); -GEOS_UNUSED_VAR(phaseDensity); -GEOS_UNUSED_VAR(phaseComponentFraction); -GEOS_UNUSED_VAR(phaseVolumeFraction); -GEOS_UNUSED_VAR(phaseTrappedVolumeFraction); -GEOS_UNUSED_VAR(phaseRelativePermeability); + GEOS_UNUSED_VAR( regionStatistics ); + GEOS_UNUSED_VAR( regionMarkers ); + GEOS_UNUSED_VAR( elementVolume ); + GEOS_UNUSED_VAR( pressure ); + GEOS_UNUSED_VAR( temperature ); + GEOS_UNUSED_VAR( referencePorosity ); + GEOS_UNUSED_VAR( porosity ); + GEOS_UNUSED_VAR( phaseDensity ); + GEOS_UNUSED_VAR( phaseComponentFraction ); + GEOS_UNUSED_VAR( phaseVolumeFraction ); + GEOS_UNUSED_VAR( phaseTrappedVolumeFraction ); + GEOS_UNUSED_VAR( phaseRelativePermeability ); + + integer const regionCount = regionIndices.size(); + + auto populateMinMaxTotal = [regionStatistics, + phaseCount, + componentCount] GEOS_HOST_DEVICE ( real64 const value, + real64 const weight, + integer region, + integer property, + integer phase = 0, + integer component = 0 ) + { + real64 * targetValue = nullptr; + targetValue = RegionStatistics::getProperty( regionStatistics, region, phaseCount, componentCount, RegionStatistics::MINIMUM, property, phase, component ); + RAJA::atomicMin< AtomicPolicy< POLICY > >( targetValue, value ); + targetValue = RegionStatistics::getProperty( regionStatistics, region, phaseCount, componentCount, RegionStatistics::MAXIMUM, property, phase, component ); + RAJA::atomicMax< AtomicPolicy< POLICY > >( targetValue, value ); + targetValue = RegionStatistics::getProperty( regionStatistics, region, phaseCount, componentCount, RegionStatistics::TOTAL, property, phase, component ); + RAJA::atomicAdd< AtomicPolicy< POLICY > >( targetValue, value ); + targetValue = RegionStatistics::getProperty( regionStatistics, region, phaseCount, componentCount, RegionStatistics::WEIGHTED_TOTAL, property, phase, component ); + RAJA::atomicAdd< AtomicPolicy< POLICY > >( targetValue, weight * value ); + return value; + }; + + forAll< POLICY >( size, [phaseCount, + componentCount, + regionCount, + regionStatistics, + regionIndices, + regionMarkers, + elementCounts, + elementGhostRank, + elementVolume, + pressure, + temperature, + referencePorosity, + porosity, + phaseDensity, + phaseComponentFraction, + phaseVolumeFraction, + phaseTrappedVolumeFraction, + phaseRelativePermeability, + &populateMinMaxTotal] GEOS_HOST_DEVICE ( localIndex const ei ) + { + if( elementGhostRank[ei] >= 0 ) + { + return; + } + + // Find the appropriate region + integer region = 0; + real64 const regionMarker = regionMarkers[ei]; + for(; region < regionCount; ++region ) + { + if( LvArray::math::abs( regionIndices[region] - regionMarker ) < 0.1 ) + { + break; + } + } + + if( regionCount <= region ) + { + return; + } + + RAJA::atomicAdd< AtomicPolicy< POLICY > >( &elementCounts[region], 1 ); + + // Uncompacted pore volume which is also used as the weight + real64 const weight = elementVolume[ei] * referencePorosity[ei]; + + populateMinMaxTotal( weight, 1.0, region, RegionStatistics::PORE_VOLUME ); + + populateMinMaxTotal( pressure[ei], weight, region, RegionStatistics::PRESSURE ); + populateMinMaxTotal( temperature[ei], weight, region, RegionStatistics::TEMPERATURE ); + + for( integer phaseIndex = 0; phaseIndex < phaseCount; ++phaseIndex ) + { + real64 const elementPhaseVolume = elementVolume[ei] * porosity[ei][0] * phaseVolumeFraction[ei][phaseIndex]; + real64 const elementPhaseDensity = phaseDensity[ei][0][phaseIndex]; + real64 const elementPhaseMass = elementPhaseVolume * elementPhaseDensity; + + populateMinMaxTotal( elementPhaseVolume, weight, region, RegionStatistics::PHASE_PORE_VOLUME, phaseIndex ); + populateMinMaxTotal( elementPhaseMass, weight, region, RegionStatistics::PHASE_MASS, phaseIndex ); + + for( integer compIndex = 0; compIndex < componentCount; ++compIndex ) + { + real64 const elementComponentPhaseMass = elementPhaseMass * phaseComponentFraction[ei][0][phaseIndex][compIndex]; + populateMinMaxTotal( elementComponentPhaseMass, weight, region, RegionStatistics::PHASE_COMP_MASS, phaseIndex, compIndex ); + } + } + } ); + + // Dummy loop to bring data back to the CPU + forAll< serialPolicy >( 1, [regionStatistics, elementCounts] ( localIndex const ) + { + GEOS_UNUSED_VAR( regionStatistics, elementCounts ); + } ); } }; diff --git a/src/coreComponents/physicsSolvers/fluidFlow/RegionMultiphaseStatistics.cpp b/src/coreComponents/physicsSolvers/fluidFlow/RegionMultiphaseStatistics.cpp index 6774f32beaa..ff55952bd7f 100644 --- a/src/coreComponents/physicsSolvers/fluidFlow/RegionMultiphaseStatistics.cpp +++ b/src/coreComponents/physicsSolvers/fluidFlow/RegionMultiphaseStatistics.cpp @@ -1,136 +1,439 @@ -/* - * ------------------------------------------------------------------------------------------------------------ - * SPDX-License-Identifier: LGPL-2.1-only - * - * Copyright (c) 2018-2020 Lawrence Livermore National Security LLC - * Copyright (c) 2018-2020 The Board of Trustees of the Leland Stanford Junior University - * Copyright (c) 2018-2020 TotalEnergies - * Copyright (c) 2019- GEOSX Contributors - * All rights reserved - * - * See top level LICENSE, COPYRIGHT, CONTRIBUTORS, NOTICE, and ACKNOWLEDGEMENTS files for details. - * ------------------------------------------------------------------------------------------------------------ - */ - -/** - * @file RegionMultiphaseStatistics.cpp - */ - -#include "RegionMultiphaseStatistics.hpp" -//#include "constitutive/fluid/multifluid/MultiFluidBase.hpp" -//#include "constitutive/relativePermeability/RelativePermeabilityBase.hpp" -//#include "constitutive/solid/CoupledSolidBase.hpp" -//#include "finiteVolume/FiniteVolumeManager.hpp" -//#include "finiteVolume/FluxApproximationBase.hpp" -//#include "mainInterface/ProblemManager.hpp" -//#include "physicsSolvers/PhysicsSolverManager.hpp" -#include "physicsSolvers/fluidFlow/CompositionalMultiphaseBase.hpp" -#include "mesh/MeshFields.hpp" -//#include "physicsSolvers/fluidFlow/CompositionalMultiphaseBaseFields.hpp" -//#include "physicsSolvers/fluidFlow/CompositionalMultiphaseHybridFVM.hpp" -//#include "physicsSolvers/fluidFlow/FlowSolverBaseFields.hpp" -//#include "physicsSolvers/fluidFlow/IsothermalCompositionalMultiphaseBaseKernels.hpp" -//#include "physicsSolvers/fluidFlow/IsothermalCompositionalMultiphaseFVMKernels.hpp" - -namespace geos -{ - -using namespace constitutive; -using namespace dataRepository; - -RegionMultiphaseStatistics::RegionMultiphaseStatistics( const string & name, - Group * const parent ): - Base( name, parent ) -{ - registerWrapper( viewKeyStruct::regionNamesString(), &m_regionNames ). - setInputFlag( InputFlags::REQUIRED ). - setDescription( "The list of names of the regions" ); - - registerWrapper( viewKeyStruct::regionIdentifiersString(), &m_regionIdentifiers ). - setInputFlag( InputFlags::REQUIRED ). - setDescription( "The list of integer identifiers of the regions" ); - - registerWrapper( viewKeyStruct::propertyNamesString(), &m_propertyNames ). - setInputFlag( InputFlags::OPTIONAL ). - setDescription( "The list of names of properties to report" ); -} - -void RegionMultiphaseStatistics::postProcessInput() -{ - Base::postProcessInput(); - - // Make sure we have the same number of region names and region indices - GEOS_ERROR_IF_NE(m_regionNames.size(), m_regionIdentifiers.size(), - GEOS_FMT( "The number of region names {} is different from the number of region identifiers {}.", - m_regionNames.size(), m_regionIdentifiers.size() ) ); - - // Convert the property names to integers - m_propertyNameTypes.clear(); - localIndex const propSize = m_propertyNames.size(); - if( propSize == 0 ) - { -m_propertyNameTypes.emplace_back(PropertyNameType::Pressure); -m_propertyNameTypes.emplace_back(PropertyNameType::PoreVolume); -m_propertyNameTypes.emplace_back(PropertyNameType::Saturation); - } - else - { - m_propertyNameTypes.resize( propSize ); - for( integer i = 0; i < propSize; ++i ) - { - m_propertyNameTypes[i] = EnumStrings< PropertyNameType >::fromString( m_propertyNames[i] ); - } - } -} -void RegionMultiphaseStatistics::registerDataOnMesh( Group & meshBodies ) -{ - Base::registerDataOnMesh( 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; - } - - string const fieldName = GEOS_FMT( "{}{}", getName(), viewKeyStruct::fieldNameString() ); - - m_solver->forDiscretizationOnMeshTargets( meshBodies, [&] ( string const &, - MeshLevel & mesh, - arrayView1d< string const > const & regionNames ) - { - ElementRegionManager & elemManager = mesh.getElemManager(); -elemManager.forElementSubRegions( regionNames, [&]( localIndex const, - ElementSubRegionBase & subRegion ) - { - subRegion.registerWrapper>( fieldName ). - setPlotLevel( PlotLevel::NOPLOT ). - setRestartFlags( RestartFlags::NO_WRITE ). - setDescription( fieldName ). - setRegisteringObjects( getName() ); - }); - }); -} - -bool RegionMultiphaseStatistics::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 ) -{ - GEOS_UNUSED_VAR( time_n ); - GEOS_UNUSED_VAR( dt ); - GEOS_UNUSED_VAR( domain ); - std::cout << "*** Executing " << getName() << "\n"; - return false; -} - -REGISTER_CATALOG_ENTRY( TaskBase, - RegionMultiphaseStatistics, - string const &, dataRepository::Group * const ) - -} /* namespace geos */ - +/* + * ------------------------------------------------------------------------------------------------------------ + * SPDX-License-Identifier: LGPL-2.1-only + * + * Copyright (c) 2018-2020 Lawrence Livermore National Security LLC + * Copyright (c) 2018-2020 The Board of Trustees of the Leland Stanford Junior University + * Copyright (c) 2018-2020 TotalEnergies + * Copyright (c) 2019- GEOSX Contributors + * All rights reserved + * + * See top level LICENSE, COPYRIGHT, CONTRIBUTORS, NOTICE, and ACKNOWLEDGEMENTS files for details. + * ------------------------------------------------------------------------------------------------------------ + */ + +/** + * @file RegionMultiphaseStatistics.cpp + */ + +#include "RegionMultiphaseStatistics.hpp" +#include "constitutive/fluid/multifluid/MultiFluidBase.hpp" +#include "constitutive/relativePermeability/RelativePermeabilityBase.hpp" +#include "constitutive/solid/CoupledSolidBase.hpp" +//#include "finiteVolume/FiniteVolumeManager.hpp" +//#include "finiteVolume/FluxApproximationBase.hpp" +//#include "mainInterface/ProblemManager.hpp" +//#include "physicsSolvers/PhysicsSolverManager.hpp" +#include "physicsSolvers/fluidFlow/CompositionalMultiphaseBase.hpp" +#include "mesh/MeshFields.hpp" +#include "physicsSolvers/fluidFlow/FlowSolverBaseFields.hpp" +#include "physicsSolvers/fluidFlow/CompositionalMultiphaseBaseFields.hpp" +//#include "physicsSolvers/fluidFlow/CompositionalMultiphaseHybridFVM.hpp" +//#include "physicsSolvers/fluidFlow/IsothermalCompositionalMultiphaseBaseKernels.hpp" +//#include "physicsSolvers/fluidFlow/IsothermalCompositionalMultiphaseFVMKernels.hpp" + +namespace geos +{ + +using namespace constitutive; +using namespace dataRepository; + +// Container for region statistics +struct RegionMultiphaseStatistics::RegionStatistics +{ + /// Indices of minimum, maximum and total property values + static constexpr integer MINIMUM = 0; + static constexpr integer MAXIMUM = 1; + static constexpr integer TOTAL = 2; + static constexpr integer WEIGHTED_TOTAL = 3; + static constexpr integer TYPE_END = 4; + + /// Indices for property type + static constexpr integer PORE_VOLUME = 0; + static constexpr integer PRESSURE = 1; + static constexpr integer TEMPERATURE = 2; + static constexpr integer PHASE_PORE_VOLUME = 3; + static constexpr integer PHASE_MASS = 4; + static constexpr integer PHASE_COMP_MASS = 5; + static constexpr integer PROP_END = 6; + + static std::unique_ptr< array3d< real64 > > allocate( integer regionCount, integer phaseCount, integer componentCount ) + { + std::unique_ptr< array3d< real64 > > properties = std::make_unique< array3d< real64 > >( + TYPE_END, + regionCount, + getSize( phaseCount, componentCount )); + LvArray::forValuesInSlice( properties->toSlice(), []( real64 & v ){ v = 0.0; } ); + + static constexpr real64 minValue = -LvArray::NumericLimits< real64 >::max; + static constexpr real64 maxValue = LvArray::NumericLimits< real64 >::max; + LvArray::forValuesInSlice( (*properties)[MINIMUM], []( real64 & v ){ v = maxValue; } ); + LvArray::forValuesInSlice( (*properties)[MAXIMUM], []( real64 & v ){ v = minValue; } ); + + return std::move( properties ); + } + + static std::unique_ptr< array1d< localIndex > > allocateCounts( integer regionCount ) + { + std::unique_ptr< array1d< localIndex > > elementCounts = std::make_unique< array1d< localIndex > >( regionCount ); + LvArray::forValuesInSlice( elementCounts->toSlice(), []( localIndex & v ){ v = 0; } ); + return std::move( elementCounts ); + } + + static real64 * getProperty( arrayView3d< real64 > props, + integer phaseCount, + integer componentCount, + integer region, + integer type, + integer property, + integer phase = 0, + integer component = 0 ) + { + return &props( type, region, calculateIndex( phaseCount, componentCount, property, phase, component )); + } + +private: + static integer getSize( integer phaseCount, integer componentCount ) + { + // 1 pore volume, 1 pressure, 1 temperature, np phase volume, np phase mass, np*nc phase comp mass + return 3 + 2*phaseCount + phaseCount*componentCount; + } + + static integer calculateIndex( integer phaseCount, + integer componentCount, + integer property, + integer phase, + integer component ) + { + if( property == PORE_VOLUME || property == PRESSURE || property == TEMPERATURE ) + { + return property; + } + else if( property == PHASE_PORE_VOLUME ) + { + return PHASE_PORE_VOLUME + phase; + } + else if( property == PHASE_MASS ) + { + return PHASE_PORE_VOLUME + phaseCount; + } + else if( property == PHASE_COMP_MASS ) + { + return PHASE_PORE_VOLUME + 2*phaseCount + phase*componentCount + component; + } + return 0; + } +}; + +// The region statistics kernel +struct RegionMultiphaseStatistics::RegionStatisticsKernel +{ + template< typename POLICY > + static void launch( localIndex const size, + integer const phaseCount, + integer const componentCount, + arrayView3d< real64 > const regionStatistics, + arrayView1d< localIndex > const & elementCounts, + arrayView1d< real64 const > const & regionIndices, + arrayView1d< real64 const > const & regionMarkers, + arrayView1d< integer const > const & elementGhostRank, + arrayView1d< real64 const > const & elementVolume, + arrayView1d< real64 const > const & pressure, + arrayView1d< real64 const > const & temperature, + arrayView1d< real64 const > const & referencePorosity, + arrayView2d< real64 const > const & porosity, + arrayView3d< real64 const, multifluid::USD_PHASE > const & phaseDensity, + arrayView4d< real64 const, multifluid::USD_PHASE_COMP > const & phaseComponentFraction, + arrayView2d< real64 const, compflow::USD_PHASE > const & phaseVolumeFraction, + arrayView3d< real64 const, relperm::USD_RELPERM > const & phaseTrappedVolumeFraction, + arrayView3d< real64 const, relperm::USD_RELPERM > const & phaseRelativePermeability ) + { + GEOS_UNUSED_VAR( regionStatistics ); + GEOS_UNUSED_VAR( regionMarkers ); + GEOS_UNUSED_VAR( elementVolume ); + GEOS_UNUSED_VAR( pressure ); + GEOS_UNUSED_VAR( temperature ); + GEOS_UNUSED_VAR( referencePorosity ); + GEOS_UNUSED_VAR( porosity ); + GEOS_UNUSED_VAR( phaseDensity ); + GEOS_UNUSED_VAR( phaseComponentFraction ); + GEOS_UNUSED_VAR( phaseVolumeFraction ); + GEOS_UNUSED_VAR( phaseTrappedVolumeFraction ); + GEOS_UNUSED_VAR( phaseRelativePermeability ); + + integer const regionCount = regionIndices.size(); + + auto populateMinMaxTotal = [regionStatistics, + phaseCount, + componentCount] GEOS_HOST_DEVICE ( real64 const value, + real64 const weight, + integer region, + integer property, + integer phase = 0, + integer component = 0 ) + { + real64 * targetValue = nullptr; + targetValue = RegionStatistics::getProperty( regionStatistics, region, phaseCount, componentCount, RegionStatistics::MINIMUM, property, phase, component ); + RAJA::atomicMin< AtomicPolicy< POLICY > >( targetValue, value ); + targetValue = RegionStatistics::getProperty( regionStatistics, region, phaseCount, componentCount, RegionStatistics::MAXIMUM, property, phase, component ); + RAJA::atomicMax< AtomicPolicy< POLICY > >( targetValue, value ); + targetValue = RegionStatistics::getProperty( regionStatistics, region, phaseCount, componentCount, RegionStatistics::TOTAL, property, phase, component ); + RAJA::atomicAdd< AtomicPolicy< POLICY > >( targetValue, value ); + targetValue = RegionStatistics::getProperty( regionStatistics, region, phaseCount, componentCount, RegionStatistics::WEIGHTED_TOTAL, property, phase, component ); + RAJA::atomicAdd< AtomicPolicy< POLICY > >( targetValue, weight * value ); + return value; + }; + + forAll< POLICY >( size, [phaseCount, + componentCount, + regionCount, + regionStatistics, + regionIndices, + regionMarkers, + elementCounts, + elementGhostRank, + elementVolume, + pressure, + temperature, + referencePorosity, + porosity, + phaseDensity, + phaseComponentFraction, + phaseVolumeFraction, + phaseTrappedVolumeFraction, + phaseRelativePermeability, + &populateMinMaxTotal] GEOS_HOST_DEVICE ( localIndex const ei ) + { + if( elementGhostRank[ei] >= 0 ) + { + return; + } + + // Find the appropriate region + integer region = 0; + real64 const regionMarker = regionMarkers[ei]; + for(; region < regionCount; ++region ) + { + if( LvArray::math::abs( regionIndices[region] - regionMarker ) < 0.1 ) + { + break; + } + } + + if( regionCount <= region ) + { + return; + } + + RAJA::atomicAdd< AtomicPolicy< POLICY > >( &elementCounts[region], 1 ); + + // Uncompacted pore volume which is also used as the weight + real64 const weight = elementVolume[ei] * referencePorosity[ei]; + + populateMinMaxTotal( weight, 1.0, region, RegionStatistics::PORE_VOLUME ); + + populateMinMaxTotal( pressure[ei], weight, region, RegionStatistics::PRESSURE ); + populateMinMaxTotal( temperature[ei], weight, region, RegionStatistics::TEMPERATURE ); + + for( integer phaseIndex = 0; phaseIndex < phaseCount; ++phaseIndex ) + { + real64 const elementPhaseVolume = elementVolume[ei] * porosity[ei][0] * phaseVolumeFraction[ei][phaseIndex]; + real64 const elementPhaseDensity = phaseDensity[ei][0][phaseIndex]; + real64 const elementPhaseMass = elementPhaseVolume * elementPhaseDensity; + + populateMinMaxTotal( elementPhaseVolume, weight, region, RegionStatistics::PHASE_PORE_VOLUME, phaseIndex ); + populateMinMaxTotal( elementPhaseMass, weight, region, RegionStatistics::PHASE_MASS, phaseIndex ); + + for( integer compIndex = 0; compIndex < componentCount; ++compIndex ) + { + real64 const elementComponentPhaseMass = elementPhaseMass * phaseComponentFraction[ei][0][phaseIndex][compIndex]; + populateMinMaxTotal( elementComponentPhaseMass, weight, region, RegionStatistics::PHASE_COMP_MASS, phaseIndex, compIndex ); + } + } + } ); + + // Dummy loop to bring data back to the CPU + forAll< serialPolicy >( 1, [regionStatistics, elementCounts] ( localIndex const ) + { + GEOS_UNUSED_VAR( regionStatistics, elementCounts ); + } ); + } +}; + +RegionMultiphaseStatistics::RegionMultiphaseStatistics( const string & name, + Group * const parent ): + Base( name, parent ) +{ + registerWrapper( viewKeyStruct::regionNamesString(), &m_regionNames ). + setInputFlag( InputFlags::REQUIRED ). + setDescription( "The list of names of the regions" ); + + registerWrapper( viewKeyStruct::regionIdentifiersString(), &m_regionIdentifiers ). + setInputFlag( InputFlags::REQUIRED ). + setDescription( "The list of integer identifiers of the regions" ); + + registerWrapper( viewKeyStruct::propertyNamesString(), &m_propertyNames ). + setInputFlag( InputFlags::OPTIONAL ). + setDescription( "The list of names of properties to report" ); +} + +void RegionMultiphaseStatistics::postProcessInput() +{ + Base::postProcessInput(); + + // Make sure we have the same number of region names and region indices + GEOS_ERROR_IF_NE_MSG( m_regionNames.size(), m_regionIdentifiers.size(), + GEOS_FMT( "The number of region names ({}) is different from the number of region identifiers ({}).", + m_regionNames.size(), m_regionIdentifiers.size() ) ); + + // Convert the property names to integers + m_propertyNameTypes.clear(); + localIndex const propSize = m_propertyNames.size(); + if( propSize == 0 ) + { + m_propertyNameTypes.emplace_back( PropertyNameType::Pressure ); + m_propertyNameTypes.emplace_back( PropertyNameType::PoreVolume ); + m_propertyNameTypes.emplace_back( PropertyNameType::Saturation ); + } + else + { + m_propertyNameTypes.resize( propSize ); + for( integer i = 0; i < propSize; ++i ) + { + m_propertyNameTypes[i] = EnumStrings< PropertyNameType >::fromString( m_propertyNames[i] ); + } + } +} + +void RegionMultiphaseStatistics::registerDataOnMesh( Group & meshBodies ) +{ + Base::registerDataOnMesh( 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; + } + + string const fieldName = GEOS_FMT( "{}{}", getName(), viewKeyStruct::fieldNameString() ); + + m_solver->forDiscretizationOnMeshTargets( meshBodies, [&] ( string const &, + MeshLevel & mesh, + arrayView1d< string const > const & regionNames ) + { + ElementRegionManager & elemManager = mesh.getElemManager(); + elemManager.forElementSubRegions( regionNames, [&]( localIndex const, + ElementSubRegionBase & subRegion ) + { + subRegion.registerWrapper< array1d< real64 > >( fieldName ). + setPlotLevel( PlotLevel::NOPLOT ). + setRestartFlags( RestartFlags::NO_WRITE ). + setDescription( fieldName ). + setRegisteringObjects( getName() ); + } ); + } ); +} + +bool RegionMultiphaseStatistics::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, + arrayView1d< string const > const & regionNames ) + { + real64 const currentTime = time_n + dt; + computeRegionStatistics( currentTime, mesh, regionNames ); + } ); + + return false; +} + +void RegionMultiphaseStatistics::computeRegionStatistics( real64 const time, + MeshLevel & mesh, + arrayView1d< string const > const & regionNames ) const +{ + GEOS_MARK_FUNCTION; + + GEOS_UNUSED_VAR( time ); + + integer const numPhases = m_solver->numFluidPhases(); + integer const numComps = m_solver->numFluidComponents(); + + string const fieldName = GEOS_FMT( "{}{}", getName(), viewKeyStruct::fieldNameString() ); + + // Step 1: initialize the average/min/max quantities + // All properties will be stored in one large array + std::unique_ptr< array3d< real64 > > regionStatistics = RegionStatistics::allocate( m_regionNames.size(), numPhases, numComps ); + std::unique_ptr< array1d< localIndex > > elementCounts = RegionStatistics::allocateCounts( m_regionNames.size() ); + + // Step 2: increment the average/min/max quantities for all the subRegions + mesh.getElemManager().forElementSubRegions( regionNames, [&]( localIndex const, + ElementSubRegionBase & subRegion ) + { + arrayView1d< integer const > const elementGhostRank = subRegion.ghostRank(); + arrayView1d< real64 const > const elementVolume = subRegion.getElementVolume(); + array1d< real64 > const regionMarkers = subRegion.getWrapper< array1d< real64 > >( fieldName ).reference(); + arrayView1d< real64 const > const pressure = subRegion.getField< fields::flow::pressure >(); + arrayView1d< real64 const > const temperature = subRegion.getField< fields::flow::temperature >(); + arrayView2d< real64 const, compflow::USD_PHASE > const phaseVolumeFraction = subRegion.getField< fields::flow::phaseVolumeFraction >(); + + 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 referencePorosity = 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 phaseComponentFraction = fluid.phaseCompFraction(); + + string const & relpermName = subRegion.getReference< string >( CompositionalMultiphaseBase::viewKeyStruct::relPermNamesString() ); + RelativePermeabilityBase const & relperm = constitutiveModels.getGroup< RelativePermeabilityBase >( relpermName ); + arrayView3d< real64 const, relperm::USD_RELPERM > const phaseTrappedVolumeFraction = relperm.phaseTrappedVolFraction(); + arrayView3d< real64 const, relperm::USD_RELPERM > const phaseRelativePermeability = relperm.phaseRelPerm(); + + RegionStatisticsKernel::launch< parallelDevicePolicy<> >( subRegion.size(), + numPhases, + numComps, + regionStatistics->toView(), + elementCounts->toView(), + m_regionIdentifiers.toViewConst(), + regionMarkers.toViewConst(), + elementGhostRank, + elementVolume, + pressure, + temperature, + referencePorosity, + porosity, + phaseDensity, + phaseComponentFraction, + phaseVolumeFraction, + phaseTrappedVolumeFraction, + phaseRelativePermeability ); + } ); + + // Step 3: synchronize the results over the MPI ranks + for( integer regionIndex = 0; regionIndex < m_regionIdentifiers.size(); ++regionIndex ) + {} + + for( integer i = 0; i < elementCounts->size(); ++i ) + { + std::cout << "Region " << i + 1 << ": " << (*elementCounts)[i] << "\n"; + } +} + +REGISTER_CATALOG_ENTRY( TaskBase, + RegionMultiphaseStatistics, + string const &, + dataRepository::Group * const ) + +} /* namespace geos */ diff --git a/src/coreComponents/physicsSolvers/fluidFlow/RegionMultiphaseStatistics.hpp b/src/coreComponents/physicsSolvers/fluidFlow/RegionMultiphaseStatistics.hpp index 49e1bf72571..9c20e7096d3 100644 --- a/src/coreComponents/physicsSolvers/fluidFlow/RegionMultiphaseStatistics.hpp +++ b/src/coreComponents/physicsSolvers/fluidFlow/RegionMultiphaseStatistics.hpp @@ -91,54 +91,28 @@ class RegionMultiphaseStatistics : public FieldStatisticsBase< CompositionalMult constexpr static char const * fieldNameString() { return "_region"; } }; - struct RegionStatistics - { - explicit RegionStatistics(integer regionCount); - /// 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; - }; - void postProcessInput() override; void registerDataOnMesh( Group & meshBodies ) override; + /** + * @brief Compute some statistics on the regions + * @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, + arrayView1d< string const > const & regionNames ) const; + + struct RegionStatistics; + struct RegionStatisticsKernel; + // The names of regions string_array m_regionNames; // The ids of regions - array1d m_regionIdentifiers; + array1d< real64 > m_regionIdentifiers; // The names of properties to output string_array m_propertyNames; From 62149b563d7ba240fad11a18f971e9ecfb548dac Mon Sep 17 00:00:00 2001 From: dkachuma Date: Mon, 11 Mar 2024 16:50:46 -0500 Subject: [PATCH 3/7] Restore kernel file --- ...rmalCompositionalMultiphaseBaseKernels.hpp | 246 ------------------ .../fluidFlow/RegionMultiphaseStatistics.cpp | 18 +- 2 files changed, 2 insertions(+), 262 deletions(-) diff --git a/src/coreComponents/physicsSolvers/fluidFlow/IsothermalCompositionalMultiphaseBaseKernels.hpp b/src/coreComponents/physicsSolvers/fluidFlow/IsothermalCompositionalMultiphaseBaseKernels.hpp index 6e972ff7ecd..8adafb2945a 100644 --- a/src/coreComponents/physicsSolvers/fluidFlow/IsothermalCompositionalMultiphaseBaseKernels.hpp +++ b/src/coreComponents/physicsSolvers/fluidFlow/IsothermalCompositionalMultiphaseBaseKernels.hpp @@ -2085,252 +2085,6 @@ struct StatisticsKernel } }; -/******************************** RegionStatisticsKernel ********************************/ - -struct RegionStatisticsKernel -{ - struct RegionStatistics - { - /// Indices of minimum, maximum and total property values - static constexpr integer MINIMUM = 0; - static constexpr integer MAXIMUM = 1; - static constexpr integer TOTAL = 2; - static constexpr integer WEIGHTED_TOTAL = 3; - static constexpr integer TYPE_END = 4; - - /// Indices for property type - static constexpr integer PORE_VOLUME = 0; - static constexpr integer PRESSURE = 1; - static constexpr integer TEMPERATURE = 2; - static constexpr integer PHASE_PORE_VOLUME = 3; - static constexpr integer PHASE_MASS = 4; - static constexpr integer PHASE_COMP_MASS = 5; - static constexpr integer PROP_END = 6; - - RegionStatistics( integer regionCount, integer phaseCount, integer componentCount ) - : m_phaseCount( phaseCount ), - m_componentCount( componentCount ), - m_properties( TYPE_END, regionCount, getSize( phaseCount, componentCount )), - m_elementCounts( regionCount ) - { - LvArray::forValuesInSlice( m_properties.toSlice(), []( real64 & v ){ v = 0.0; } ); - LvArray::forValuesInSlice( m_elementCounts.toSlice(), []( localIndex & i ){ i = 0; } ); - - static constexpr real64 minValue = -LvArray::NumericLimits< real64 >::max; - static constexpr real64 maxValue = LvArray::NumericLimits< real64 >::max; - LvArray::forValuesInSlice( m_properties[MINIMUM], []( real64 & v ){ v = maxValue; } ); - LvArray::forValuesInSlice( m_properties[MAXIMUM], []( real64 & v ){ v = minValue; } ); - } - - real64 & getProperty( integer region, integer type, integer property, integer phase = 0, integer component = 0 ) - { - return m_properties( type, region, calculateIndex( property, phase, component )); - } - - real64 getProperty( integer region, integer type, integer property, integer phase = 0, integer component = 0 ) const - { - return m_properties( type, region, calculateIndex( property, phase, component )); - } - - static real64 * getProperty( arrayView3d< real64 > props, - integer phaseCount, - integer componentCount, - integer region, - integer type, - integer property, - integer phase = 0, - integer component = 0 ) - { - return &props( type, region, calculateIndex( phaseCount, componentCount, property, phase, component )); - } - - arrayView3d< real64 > toView() const - { - return m_properties.toView(); - } - - arrayView1d< localIndex > getElementCounts() const - { - return m_elementCounts.toView(); - } - -private: - static integer getSize( integer phaseCount, integer componentCount ) - { - // 1 pore volume, 1 pressure, 1 temperature, np phase volume, np phase mass, np*nc phase comp mass - return 3 + 2*phaseCount + phaseCount*componentCount; - } - - integer calculateIndex( integer property, integer phase, integer component ) const - { - return calculateIndex( m_phaseCount, m_componentCount, property, phase, component ); - } - - static integer calculateIndex( integer phaseCount, - integer componentCount, - integer property, - integer phase, - integer component ) - { - if( property == PORE_VOLUME || property == PRESSURE || property == TEMPERATURE ) - { - return property; - } - else if( property == PHASE_PORE_VOLUME ) - { - return PHASE_PORE_VOLUME + phase; - } - else if( property == PHASE_MASS ) - { - return PHASE_PORE_VOLUME + phaseCount; - } - else if( property == PHASE_COMP_MASS ) - { - return PHASE_PORE_VOLUME + 2*phaseCount + phase*componentCount + component; - } - return 0; - } - - integer const m_phaseCount{0}; - integer const m_componentCount{0}; - array3d< real64 > m_properties; - array1d< localIndex > m_elementCounts; - }; - - template< typename POLICY > - static void - launch( localIndex const size, - integer const phaseCount, - integer const componentCount, - arrayView3d< real64 > const regionStatistics, - arrayView1d< localIndex > const & elementCounts, - arrayView1d< real64 const > const & regionIndices, - arrayView1d< real64 const > const & regionMarkers, - arrayView1d< integer const > const & elementGhostRank, - arrayView1d< real64 const > const & elementVolume, - arrayView1d< real64 const > const & pressure, - arrayView1d< real64 const > const & temperature, - arrayView1d< real64 const > const & referencePorosity, - arrayView2d< real64 const > const & porosity, - arrayView3d< real64 const, multifluid::USD_PHASE > const & phaseDensity, - arrayView4d< real64 const, multifluid::USD_PHASE_COMP > const & phaseComponentFraction, - arrayView2d< real64 const, compflow::USD_PHASE > const & phaseVolumeFraction, - arrayView3d< real64 const, relperm::USD_RELPERM > const & phaseTrappedVolumeFraction, - arrayView3d< real64 const, relperm::USD_RELPERM > const & phaseRelativePermeability ) - { - GEOS_UNUSED_VAR( regionStatistics ); - GEOS_UNUSED_VAR( regionMarkers ); - GEOS_UNUSED_VAR( elementVolume ); - GEOS_UNUSED_VAR( pressure ); - GEOS_UNUSED_VAR( temperature ); - GEOS_UNUSED_VAR( referencePorosity ); - GEOS_UNUSED_VAR( porosity ); - GEOS_UNUSED_VAR( phaseDensity ); - GEOS_UNUSED_VAR( phaseComponentFraction ); - GEOS_UNUSED_VAR( phaseVolumeFraction ); - GEOS_UNUSED_VAR( phaseTrappedVolumeFraction ); - GEOS_UNUSED_VAR( phaseRelativePermeability ); - - integer const regionCount = regionIndices.size(); - - auto populateMinMaxTotal = [regionStatistics, - phaseCount, - componentCount] GEOS_HOST_DEVICE ( real64 const value, - real64 const weight, - integer region, - integer property, - integer phase = 0, - integer component = 0 ) - { - real64 * targetValue = nullptr; - targetValue = RegionStatistics::getProperty( regionStatistics, region, phaseCount, componentCount, RegionStatistics::MINIMUM, property, phase, component ); - RAJA::atomicMin< AtomicPolicy< POLICY > >( targetValue, value ); - targetValue = RegionStatistics::getProperty( regionStatistics, region, phaseCount, componentCount, RegionStatistics::MAXIMUM, property, phase, component ); - RAJA::atomicMax< AtomicPolicy< POLICY > >( targetValue, value ); - targetValue = RegionStatistics::getProperty( regionStatistics, region, phaseCount, componentCount, RegionStatistics::TOTAL, property, phase, component ); - RAJA::atomicAdd< AtomicPolicy< POLICY > >( targetValue, value ); - targetValue = RegionStatistics::getProperty( regionStatistics, region, phaseCount, componentCount, RegionStatistics::WEIGHTED_TOTAL, property, phase, component ); - RAJA::atomicAdd< AtomicPolicy< POLICY > >( targetValue, weight * value ); - return value; - }; - - forAll< POLICY >( size, [phaseCount, - componentCount, - regionCount, - regionStatistics, - regionIndices, - regionMarkers, - elementCounts, - elementGhostRank, - elementVolume, - pressure, - temperature, - referencePorosity, - porosity, - phaseDensity, - phaseComponentFraction, - phaseVolumeFraction, - phaseTrappedVolumeFraction, - phaseRelativePermeability, - &populateMinMaxTotal] GEOS_HOST_DEVICE ( localIndex const ei ) - { - if( elementGhostRank[ei] >= 0 ) - { - return; - } - - // Find the appropriate region - integer region = 0; - real64 const regionMarker = regionMarkers[ei]; - for(; region < regionCount; ++region ) - { - if( LvArray::math::abs( regionIndices[region] - regionMarker ) < 0.1 ) - { - break; - } - } - - if( regionCount <= region ) - { - return; - } - - RAJA::atomicAdd< AtomicPolicy< POLICY > >( &elementCounts[region], 1 ); - - // Uncompacted pore volume which is also used as the weight - real64 const weight = elementVolume[ei] * referencePorosity[ei]; - - populateMinMaxTotal( weight, 1.0, region, RegionStatistics::PORE_VOLUME ); - - populateMinMaxTotal( pressure[ei], weight, region, RegionStatistics::PRESSURE ); - populateMinMaxTotal( temperature[ei], weight, region, RegionStatistics::TEMPERATURE ); - - for( integer phaseIndex = 0; phaseIndex < phaseCount; ++phaseIndex ) - { - real64 const elementPhaseVolume = elementVolume[ei] * porosity[ei][0] * phaseVolumeFraction[ei][phaseIndex]; - real64 const elementPhaseDensity = phaseDensity[ei][0][phaseIndex]; - real64 const elementPhaseMass = elementPhaseVolume * elementPhaseDensity; - - populateMinMaxTotal( elementPhaseVolume, weight, region, RegionStatistics::PHASE_PORE_VOLUME, phaseIndex ); - populateMinMaxTotal( elementPhaseMass, weight, region, RegionStatistics::PHASE_MASS, phaseIndex ); - - for( integer compIndex = 0; compIndex < componentCount; ++compIndex ) - { - real64 const elementComponentPhaseMass = elementPhaseMass * phaseComponentFraction[ei][0][phaseIndex][compIndex]; - populateMinMaxTotal( elementComponentPhaseMass, weight, region, RegionStatistics::PHASE_COMP_MASS, phaseIndex, compIndex ); - } - } - } ); - - // Dummy loop to bring data back to the CPU - forAll< serialPolicy >( 1, [regionStatistics, elementCounts] ( localIndex const ) - { - GEOS_UNUSED_VAR( regionStatistics, elementCounts ); - } ); - } -}; - - /******************************** HydrostaticPressureKernel ********************************/ struct HydrostaticPressureKernel diff --git a/src/coreComponents/physicsSolvers/fluidFlow/RegionMultiphaseStatistics.cpp b/src/coreComponents/physicsSolvers/fluidFlow/RegionMultiphaseStatistics.cpp index ff55952bd7f..1424da5fa7d 100644 --- a/src/coreComponents/physicsSolvers/fluidFlow/RegionMultiphaseStatistics.cpp +++ b/src/coreComponents/physicsSolvers/fluidFlow/RegionMultiphaseStatistics.cpp @@ -148,19 +148,6 @@ struct RegionMultiphaseStatistics::RegionStatisticsKernel arrayView3d< real64 const, relperm::USD_RELPERM > const & phaseTrappedVolumeFraction, arrayView3d< real64 const, relperm::USD_RELPERM > const & phaseRelativePermeability ) { - GEOS_UNUSED_VAR( regionStatistics ); - GEOS_UNUSED_VAR( regionMarkers ); - GEOS_UNUSED_VAR( elementVolume ); - GEOS_UNUSED_VAR( pressure ); - GEOS_UNUSED_VAR( temperature ); - GEOS_UNUSED_VAR( referencePorosity ); - GEOS_UNUSED_VAR( porosity ); - GEOS_UNUSED_VAR( phaseDensity ); - GEOS_UNUSED_VAR( phaseComponentFraction ); - GEOS_UNUSED_VAR( phaseVolumeFraction ); - GEOS_UNUSED_VAR( phaseTrappedVolumeFraction ); - GEOS_UNUSED_VAR( phaseRelativePermeability ); - integer const regionCount = regionIndices.size(); auto populateMinMaxTotal = [regionStatistics, @@ -324,9 +311,8 @@ void RegionMultiphaseStatistics::registerDataOnMesh( Group & meshBodies ) MeshLevel & mesh, arrayView1d< string const > const & regionNames ) { - ElementRegionManager & elemManager = mesh.getElemManager(); - elemManager.forElementSubRegions( regionNames, [&]( localIndex const, - ElementSubRegionBase & subRegion ) + mesh.getElemManager().forElementSubRegions( regionNames, [&]( localIndex const, + ElementSubRegionBase & subRegion ) { subRegion.registerWrapper< array1d< real64 > >( fieldName ). setPlotLevel( PlotLevel::NOPLOT ). From 1f8cbe9886de07b726b5cf5e4bfd7978cb81e64e Mon Sep 17 00:00:00 2001 From: dkachuma Date: Tue, 12 Mar 2024 16:04:34 -0500 Subject: [PATCH 4/7] Implement averaging --- .../fluidFlow/RegionMultiphaseStatistics.cpp | 336 +++++++++++++++--- .../fluidFlow/RegionMultiphaseStatistics.hpp | 15 +- .../docs/RegionMultiphaseStatistics.rst | 22 +- src/coreComponents/schema/schema.xsd | 2 +- 4 files changed, 309 insertions(+), 66 deletions(-) diff --git a/src/coreComponents/physicsSolvers/fluidFlow/RegionMultiphaseStatistics.cpp b/src/coreComponents/physicsSolvers/fluidFlow/RegionMultiphaseStatistics.cpp index 1424da5fa7d..bad58bcd56b 100644 --- a/src/coreComponents/physicsSolvers/fluidFlow/RegionMultiphaseStatistics.cpp +++ b/src/coreComponents/physicsSolvers/fluidFlow/RegionMultiphaseStatistics.cpp @@ -17,20 +17,14 @@ */ #include "RegionMultiphaseStatistics.hpp" +#include "constitutive/ConstitutiveManager.hpp" #include "constitutive/fluid/multifluid/MultiFluidBase.hpp" #include "constitutive/relativePermeability/RelativePermeabilityBase.hpp" #include "constitutive/solid/CoupledSolidBase.hpp" -//#include "finiteVolume/FiniteVolumeManager.hpp" -//#include "finiteVolume/FluxApproximationBase.hpp" -//#include "mainInterface/ProblemManager.hpp" -//#include "physicsSolvers/PhysicsSolverManager.hpp" #include "physicsSolvers/fluidFlow/CompositionalMultiphaseBase.hpp" #include "mesh/MeshFields.hpp" #include "physicsSolvers/fluidFlow/FlowSolverBaseFields.hpp" #include "physicsSolvers/fluidFlow/CompositionalMultiphaseBaseFields.hpp" -//#include "physicsSolvers/fluidFlow/CompositionalMultiphaseHybridFVM.hpp" -//#include "physicsSolvers/fluidFlow/IsothermalCompositionalMultiphaseBaseKernels.hpp" -//#include "physicsSolvers/fluidFlow/IsothermalCompositionalMultiphaseFVMKernels.hpp" namespace geos { @@ -45,8 +39,7 @@ struct RegionMultiphaseStatistics::RegionStatistics static constexpr integer MINIMUM = 0; static constexpr integer MAXIMUM = 1; static constexpr integer TOTAL = 2; - static constexpr integer WEIGHTED_TOTAL = 3; - static constexpr integer TYPE_END = 4; + static constexpr integer TYPE_END = 3; /// Indices for property type static constexpr integer PORE_VOLUME = 0; @@ -55,41 +48,40 @@ struct RegionMultiphaseStatistics::RegionStatistics static constexpr integer PHASE_PORE_VOLUME = 3; static constexpr integer PHASE_MASS = 4; static constexpr integer PHASE_COMP_MASS = 5; - static constexpr integer PROP_END = 6; - static std::unique_ptr< array3d< real64 > > allocate( integer regionCount, integer phaseCount, integer componentCount ) + static localIndex allocate( array3d< real64 > & properties, integer regionCount, integer phaseCount, integer componentCount ) { - std::unique_ptr< array3d< real64 > > properties = std::make_unique< array3d< real64 > >( - TYPE_END, - regionCount, - getSize( phaseCount, componentCount )); - LvArray::forValuesInSlice( properties->toSlice(), []( real64 & v ){ v = 0.0; } ); + integer const sizes[3] = { TYPE_END, regionCount, getSize( phaseCount, componentCount ) }; + properties.resize( 3, sizes ); + + LvArray::forValuesInSlice( properties.toSlice(), []( real64 & v ){ v = 0.0; } ); static constexpr real64 minValue = -LvArray::NumericLimits< real64 >::max; static constexpr real64 maxValue = LvArray::NumericLimits< real64 >::max; - LvArray::forValuesInSlice( (*properties)[MINIMUM], []( real64 & v ){ v = maxValue; } ); - LvArray::forValuesInSlice( (*properties)[MAXIMUM], []( real64 & v ){ v = minValue; } ); + LvArray::forValuesInSlice( properties[MINIMUM], []( real64 & v ){ v = maxValue; } ); + LvArray::forValuesInSlice( properties[MAXIMUM], []( real64 & v ){ v = minValue; } ); - return std::move( properties ); + return properties.size(); } - static std::unique_ptr< array1d< localIndex > > allocateCounts( integer regionCount ) + static localIndex allocateCounts( array1d< localIndex > & elementCounts, integer regionCount ) { - std::unique_ptr< array1d< localIndex > > elementCounts = std::make_unique< array1d< localIndex > >( regionCount ); - LvArray::forValuesInSlice( elementCounts->toSlice(), []( localIndex & v ){ v = 0; } ); - return std::move( elementCounts ); + elementCounts.resize( regionCount ); + LvArray::forValuesInSlice( elementCounts.toSlice(), []( localIndex & v ){ v = 0; } ); + return elementCounts.size(); } static real64 * getProperty( arrayView3d< real64 > props, + integer region, integer phaseCount, integer componentCount, - integer region, integer type, integer property, integer phase = 0, integer component = 0 ) { - return &props( type, region, calculateIndex( phaseCount, componentCount, property, phase, component )); + integer const ei = calculateIndex( phaseCount, componentCount, property, phase, component ); + return &props( type, region, ei ); } private: @@ -115,7 +107,7 @@ struct RegionMultiphaseStatistics::RegionStatistics } else if( property == PHASE_MASS ) { - return PHASE_PORE_VOLUME + phaseCount; + return PHASE_PORE_VOLUME + phaseCount + phase; } else if( property == PHASE_COMP_MASS ) { @@ -152,21 +144,21 @@ struct RegionMultiphaseStatistics::RegionStatisticsKernel auto populateMinMaxTotal = [regionStatistics, phaseCount, - componentCount] GEOS_HOST_DEVICE ( real64 const value, - real64 const weight, - integer region, - integer property, - integer phase = 0, - integer component = 0 ) + componentCount] + GEOS_HOST_DEVICE ( real64 const value, + real64 const weight, + integer region, + integer property, + integer phase = 0, + integer component = 0 ) { real64 * targetValue = nullptr; + // These atomics really make the whole point of threading quite useless here targetValue = RegionStatistics::getProperty( regionStatistics, region, phaseCount, componentCount, RegionStatistics::MINIMUM, property, phase, component ); RAJA::atomicMin< AtomicPolicy< POLICY > >( targetValue, value ); targetValue = RegionStatistics::getProperty( regionStatistics, region, phaseCount, componentCount, RegionStatistics::MAXIMUM, property, phase, component ); RAJA::atomicMax< AtomicPolicy< POLICY > >( targetValue, value ); targetValue = RegionStatistics::getProperty( regionStatistics, region, phaseCount, componentCount, RegionStatistics::TOTAL, property, phase, component ); - RAJA::atomicAdd< AtomicPolicy< POLICY > >( targetValue, value ); - targetValue = RegionStatistics::getProperty( regionStatistics, region, phaseCount, componentCount, RegionStatistics::WEIGHTED_TOTAL, property, phase, component ); RAJA::atomicAdd< AtomicPolicy< POLICY > >( targetValue, weight * value ); return value; }; @@ -228,13 +220,13 @@ struct RegionMultiphaseStatistics::RegionStatisticsKernel real64 const elementPhaseDensity = phaseDensity[ei][0][phaseIndex]; real64 const elementPhaseMass = elementPhaseVolume * elementPhaseDensity; - populateMinMaxTotal( elementPhaseVolume, weight, region, RegionStatistics::PHASE_PORE_VOLUME, phaseIndex ); - populateMinMaxTotal( elementPhaseMass, weight, region, RegionStatistics::PHASE_MASS, phaseIndex ); + populateMinMaxTotal( elementPhaseVolume, 1.0, region, RegionStatistics::PHASE_PORE_VOLUME, phaseIndex ); + populateMinMaxTotal( elementPhaseMass, 1.0, region, RegionStatistics::PHASE_MASS, phaseIndex ); for( integer compIndex = 0; compIndex < componentCount; ++compIndex ) { real64 const elementComponentPhaseMass = elementPhaseMass * phaseComponentFraction[ei][0][phaseIndex][compIndex]; - populateMinMaxTotal( elementComponentPhaseMass, weight, region, RegionStatistics::PHASE_COMP_MASS, phaseIndex, compIndex ); + populateMinMaxTotal( elementComponentPhaseMass, 1.0, region, RegionStatistics::PHASE_COMP_MASS, phaseIndex, compIndex ); } } } ); @@ -279,8 +271,10 @@ void RegionMultiphaseStatistics::postProcessInput() if( propSize == 0 ) { m_propertyNameTypes.emplace_back( PropertyNameType::Pressure ); + m_propertyNameTypes.emplace_back( PropertyNameType::Temperature ); m_propertyNameTypes.emplace_back( PropertyNameType::PoreVolume ); - m_propertyNameTypes.emplace_back( PropertyNameType::Saturation ); + m_propertyNameTypes.emplace_back( PropertyNameType::VolumeFraction ); + m_propertyNameTypes.emplace_back( PropertyNameType::Mass ); } else { @@ -292,6 +286,22 @@ void RegionMultiphaseStatistics::postProcessInput() } } +template< typename ARRAY > +std::ostream & RegionMultiphaseStatistics::writeArray( std::ostream & os, ARRAY const & array ) const +{ + integer size = array.size(); + if( 0 < size ) + { + os << array[0]; + for( integer i = 1; i < array.size(); ++i ) + { + os << "," << array[i]; + } + os << "\n"; + } + return os; +} + void RegionMultiphaseStatistics::registerDataOnMesh( Group & meshBodies ) { Base::registerDataOnMesh( meshBodies ); @@ -321,6 +331,12 @@ void RegionMultiphaseStatistics::registerDataOnMesh( Group & meshBodies ) setRegisteringObjects( getName() ); } ); } ); + + // Create the csv file if requires + if( 0 < getLogLevel() && 0 < m_writeCSV && MpiWrapper::commRank() == 0 ) + { + initializeFile(); + } } bool RegionMultiphaseStatistics::execute( real64 const time_n, @@ -330,6 +346,12 @@ bool RegionMultiphaseStatistics::execute( real64 const time_n, real64 const GEOS_UNUSED_PARAM( eventProgress ), DomainPartition & domain ) { + // If the logLevel is low and we don't need the csv then there's nothing to do + if( getLogLevel() <= 0 && m_writeCSV == 0 ) + { + return false; + } + m_solver->forDiscretizationOnMeshTargets( domain.getMeshBodies(), [&] ( string const &, MeshLevel & mesh, arrayView1d< string const > const & regionNames ) @@ -347,17 +369,17 @@ void RegionMultiphaseStatistics::computeRegionStatistics( real64 const time, { GEOS_MARK_FUNCTION; - GEOS_UNUSED_VAR( time ); - - integer const numPhases = m_solver->numFluidPhases(); - integer const numComps = m_solver->numFluidComponents(); + integer const phaseCount = m_solver->numFluidPhases(); + integer const componentCount = m_solver->numFluidComponents(); string const fieldName = GEOS_FMT( "{}{}", getName(), viewKeyStruct::fieldNameString() ); // Step 1: initialize the average/min/max quantities // All properties will be stored in one large array - std::unique_ptr< array3d< real64 > > regionStatistics = RegionStatistics::allocate( m_regionNames.size(), numPhases, numComps ); - std::unique_ptr< array1d< localIndex > > elementCounts = RegionStatistics::allocateCounts( m_regionNames.size() ); + array3d< real64 > regionStatistics; + array1d< localIndex > elementCounts; + RegionStatistics::allocate( regionStatistics, m_regionNames.size(), phaseCount, componentCount ); + RegionStatistics::allocateCounts( elementCounts, m_regionNames.size() ); // Step 2: increment the average/min/max quantities for all the subRegions mesh.getElemManager().forElementSubRegions( regionNames, [&]( localIndex const, @@ -388,10 +410,10 @@ void RegionMultiphaseStatistics::computeRegionStatistics( real64 const time, arrayView3d< real64 const, relperm::USD_RELPERM > const phaseRelativePermeability = relperm.phaseRelPerm(); RegionStatisticsKernel::launch< parallelDevicePolicy<> >( subRegion.size(), - numPhases, - numComps, - regionStatistics->toView(), - elementCounts->toView(), + phaseCount, + componentCount, + regionStatistics.toView(), + elementCounts.toView(), m_regionIdentifiers.toViewConst(), regionMarkers.toViewConst(), elementGhostRank, @@ -408,13 +430,223 @@ void RegionMultiphaseStatistics::computeRegionStatistics( real64 const time, } ); // Step 3: synchronize the results over the MPI ranks - for( integer regionIndex = 0; regionIndex < m_regionIdentifiers.size(); ++regionIndex ) - {} + integer const regionCount = m_regionIdentifiers.size(); + for( integer region = 0; region < regionCount; ++region ) + { + elementCounts[region] = MpiWrapper::sum( elementCounts[region] ); + } - for( integer i = 0; i < elementCounts->size(); ++i ) + localIndex const size1 = regionStatistics.size( 1 ); + localIndex const size2 = regionStatistics.size( 2 ); + for( localIndex i1 = 0; i1 < size1; ++i1 ) { - std::cout << "Region " << i + 1 << ": " << (*elementCounts)[i] << "\n"; + for( localIndex i2 = 0; i2 < size2; ++i2 ) + { + regionStatistics( RegionStatistics::MINIMUM, i1, i2 ) = MpiWrapper::min( regionStatistics( RegionStatistics::MINIMUM, i1, i2 ) ); + regionStatistics( RegionStatistics::MAXIMUM, i1, i2 ) = MpiWrapper::max( regionStatistics( RegionStatistics::MAXIMUM, i1, i2 ) ); + regionStatistics( RegionStatistics::TOTAL, i1, i2 ) = MpiWrapper::sum( regionStatistics( RegionStatistics::TOTAL, i1, i2 ) ); + } } + + array1d< real64 > phaseVolumes( phaseCount ); + real64 * propValue = nullptr; + + array1d< real64 > propColumns; + propColumns.emplace_back( time ); + for( integer region = 0; region < regionCount; ++region ) + { + // Extract the uncompacted region pore volume as a weight + propValue = RegionStatistics::getProperty( regionStatistics, region, phaseCount, componentCount, RegionStatistics::TOTAL, RegionStatistics::PORE_VOLUME ); + real64 const regionPoreVolume = *propValue; + real64 const inverseRegionPoreVolume = regionPoreVolume < LvArray::NumericLimits< real64 >::epsilon ? 0.0 : 1.0 / regionPoreVolume; + + for( auto const & prop : m_propertyNameTypes ) + { + if( prop == PropertyNameType::Pressure ) + { + propValue = RegionStatistics::getProperty( regionStatistics, region, phaseCount, componentCount, RegionStatistics::MINIMUM, RegionStatistics::PRESSURE ); + propColumns.emplace_back( *propValue ); + + propValue = RegionStatistics::getProperty( regionStatistics, region, phaseCount, componentCount, RegionStatistics::TOTAL, RegionStatistics::PRESSURE ); + real64 const avgPressure = (*propValue) * inverseRegionPoreVolume; + propColumns.emplace_back( avgPressure ); + + propValue = RegionStatistics::getProperty( regionStatistics, region, phaseCount, componentCount, RegionStatistics::MAXIMUM, RegionStatistics::PRESSURE ); + propColumns.emplace_back( *propValue ); + } + + if( prop == PropertyNameType::Temperature ) + { + propValue = RegionStatistics::getProperty( regionStatistics, region, phaseCount, componentCount, RegionStatistics::MINIMUM, RegionStatistics::TEMPERATURE ); + propColumns.emplace_back( *propValue ); + + propValue = RegionStatistics::getProperty( regionStatistics, region, phaseCount, componentCount, RegionStatistics::TOTAL, RegionStatistics::TEMPERATURE ); + real64 const avgPressure = (*propValue) * inverseRegionPoreVolume; + propColumns.emplace_back( avgPressure ); + + propValue = RegionStatistics::getProperty( regionStatistics, region, phaseCount, componentCount, RegionStatistics::MAXIMUM, RegionStatistics::TEMPERATURE ); + propColumns.emplace_back( *propValue ); + } + + if( prop == PropertyNameType::PoreVolume ) + { + for( integer phaseIndex = 0; phaseIndex < phaseCount; ++phaseIndex ) + { + propValue = RegionStatistics::getProperty( regionStatistics, region, phaseCount, componentCount, RegionStatistics::TOTAL, RegionStatistics::PHASE_PORE_VOLUME, phaseIndex ); + propColumns.emplace_back( *propValue ); + } + } + + if( prop == PropertyNameType::Mass ) + { + for( integer phaseIndex = 0; phaseIndex < phaseCount; ++phaseIndex ) + { + propValue = RegionStatistics::getProperty( regionStatistics, region, phaseCount, componentCount, RegionStatistics::TOTAL, RegionStatistics::PHASE_MASS, phaseIndex ); + propColumns.emplace_back( *propValue ); + } + } + + if( prop == PropertyNameType::VolumeFraction ) + { + real64 totalPhaseVolume = 0.0; + for( integer phaseIndex = 0; phaseIndex < phaseCount; ++phaseIndex ) + { + propValue = RegionStatistics::getProperty( regionStatistics, region, phaseCount, componentCount, RegionStatistics::TOTAL, RegionStatistics::PHASE_PORE_VOLUME, phaseIndex ); + phaseVolumes[phaseIndex] = *propValue; + totalPhaseVolume += phaseVolumes[phaseIndex]; + } + real64 const inverseTotalPhaseVolume = totalPhaseVolume < LvArray::NumericLimits< real64 >::epsilon ? 0.0 : 1.0 / totalPhaseVolume; + for( integer phaseIndex = 0; phaseIndex < phaseCount; ++phaseIndex ) + { + propColumns.emplace_back( phaseVolumes[phaseIndex] * inverseTotalPhaseVolume ); + } + } + } + } + + if( m_writeCSV > 0 && MpiWrapper::commRank() == 0 ) + { + std::ofstream outputFile( m_outputDir + "/" + viewKeyStruct::fileNameString() + ".csv", std::ios_base::app ); + writeArray( outputFile, propColumns ); + outputFile.close(); + } +} + +void RegionMultiphaseStatistics::initializeFile() const +{ + string const fluidModelName = m_solver->referenceFluidModelName(); + + DomainPartition const & domain = getGroupByPath< DomainPartition >( "/Problem/domain" ); + ConstitutiveManager const & constitutiveManager = domain.getConstitutiveManager(); + MultiFluidBase const & fluidModel = constitutiveManager.getConstitutiveRelation< MultiFluidBase >( fluidModelName ); + auto fluidComponentNames = fluidModel.componentNames(); + auto fluidPhaseNames = fluidModel.phaseNames(); + + integer const phaseCount = m_solver->numFluidPhases(); + integer const regionCount = m_regionIdentifiers.size(); + auto regionNames = m_regionNames.toView(); + + string_array propNames; + string_array regNames; + string_array propUnits; + string_array phaseNames; + string_array compNames; + + auto addColumn = [&]( string const propName, + string const propUnit, + integer const regionIndex=-1, + integer const phaseIndex=-1, + integer const compIndex=-1 ){ + propNames.emplace_back( propName ); + propUnits.emplace_back( propUnit ); + if( 0 <= regionIndex ) + { + regNames.emplace_back( regionNames[regionIndex] ); + } + else + { + regNames.emplace_back( "" ); + } + if( 0 <= phaseIndex ) + { + phaseNames.emplace_back( fluidPhaseNames[phaseIndex] ); + } + else + { + phaseNames.emplace_back( "" ); + } + if( 0 <= compIndex ) + { + compNames.emplace_back( fluidComponentNames[compIndex] ); + } + else + { + compNames.emplace_back( "" ); + } + }; + + integer const useMass = m_solver->getReference< integer >( CompositionalMultiphaseBase::viewKeyStruct::useMassFlagString() ); + string const massUnit = useMass ? "kg" : "mol"; + + addColumn( "Time", "s" ); + for( integer region = 0; region < regionCount; ++region ) + { + for( auto const & prop : m_propertyNameTypes ) + { + string const propName = EnumStrings< PropertyNameType >::toString( prop ); + string propUnit = ""; + if( prop == PropertyNameType::Pressure ) + { + propUnit = "Pa"; + addColumn( GEOS_FMT( "Min {}", propName ), propUnit, region ); + addColumn( GEOS_FMT( "Average {}", propName ), propUnit, region ); + addColumn( GEOS_FMT( "Max {}", propName ), propUnit, region ); + } + + if( prop == PropertyNameType::Temperature ) + { + propUnit = "K"; + addColumn( GEOS_FMT( "Min {}", propName ), propUnit, region ); + addColumn( GEOS_FMT( "Average {}", propName ), propUnit, region ); + addColumn( GEOS_FMT( "Max {}", propName ), propUnit, region ); + } + + if( prop == PropertyNameType::PoreVolume ) + { + propUnit = "rm^3"; + for( integer phaseIndex = 0; phaseIndex < phaseCount; ++phaseIndex ) + { + addColumn( propName, propUnit, region, phaseIndex ); + } + } + + if( prop == PropertyNameType::Mass ) + { + for( integer phaseIndex = 0; phaseIndex < phaseCount; ++phaseIndex ) + { + addColumn( propName, massUnit, region, phaseIndex ); + } + } + + if( prop == PropertyNameType::VolumeFraction ) + { + for( integer phaseIndex = 0; phaseIndex < phaseCount; ++phaseIndex ) + { + addColumn( propName, propUnit, region, phaseIndex ); + } + } + } + } + + std::ofstream outputFile( m_outputDir + "/" + viewKeyStruct::fileNameString() + ".csv" ); + + writeArray( outputFile, propNames ); + writeArray( outputFile, propUnits ); + writeArray( outputFile, regNames ); + writeArray( outputFile, phaseNames ); + writeArray( outputFile, compNames ); + + outputFile.close(); } REGISTER_CATALOG_ENTRY( TaskBase, diff --git a/src/coreComponents/physicsSolvers/fluidFlow/RegionMultiphaseStatistics.hpp b/src/coreComponents/physicsSolvers/fluidFlow/RegionMultiphaseStatistics.hpp index 9c20e7096d3..60322442ed7 100644 --- a/src/coreComponents/physicsSolvers/fluidFlow/RegionMultiphaseStatistics.hpp +++ b/src/coreComponents/physicsSolvers/fluidFlow/RegionMultiphaseStatistics.hpp @@ -39,8 +39,10 @@ class RegionMultiphaseStatistics : public FieldStatisticsBase< CompositionalMult enum PropertyNameType : integer { Pressure, + Temperature, PoreVolume, - Saturation + VolumeFraction, + Mass }; public: @@ -89,6 +91,8 @@ class RegionMultiphaseStatistics : public FieldStatisticsBase< CompositionalMult constexpr static char const * propertyNamesString() { return "propertyNames"; } /// String for the suffix of the field name constexpr static char const * fieldNameString() { return "_region"; } + /// String for the file name + constexpr static char const * fileNameString() { return "values"; } }; void postProcessInput() override; @@ -105,6 +109,11 @@ class RegionMultiphaseStatistics : public FieldStatisticsBase< CompositionalMult MeshLevel & mesh, arrayView1d< string const > const & regionNames ) const; + void initializeFile() const; + + template< typename ARRAY > + std::ostream & writeArray( std::ostream & os, ARRAY const & array ) const; + struct RegionStatistics; struct RegionStatisticsKernel; @@ -123,8 +132,10 @@ class RegionMultiphaseStatistics : public FieldStatisticsBase< CompositionalMult ENUM_STRINGS( RegionMultiphaseStatistics::PropertyNameType, "pressure", + "temperature", "poreVolume", - "saturation" ); + "volumeFraction", + "mass" ); } /* namespace geos */ diff --git a/src/coreComponents/schema/docs/RegionMultiphaseStatistics.rst b/src/coreComponents/schema/docs/RegionMultiphaseStatistics.rst index 245e94393dc..e4cd46cb3f2 100644 --- a/src/coreComponents/schema/docs/RegionMultiphaseStatistics.rst +++ b/src/coreComponents/schema/docs/RegionMultiphaseStatistics.rst @@ -1,15 +1,15 @@ -================= ============= ======== ============================================== -Name Type Default Description -================= ============= ======== ============================================== -flowSolverName groupNameRef required Name of the flow solver -logLevel integer 0 Log level -name groupName required A name is required for any non-unique nodes -propertyNames string_array {} The list of names of properties to report -regionIdentifiers integer_array required The list of integer identifiers of the regions -regionNames string_array required The list of names of the regions -writeCSV integer 0 Write statistics into a CSV file -================= ============= ======== ============================================== +================= ============ ======== ============================================== +Name Type Default Description +================= ============ ======== ============================================== +flowSolverName groupNameRef required Name of the flow solver +logLevel integer 0 Log level +name groupName required A name is required for any non-unique nodes +propertyNames string_array {} The list of names of properties to report +regionIdentifiers real64_array required The list of integer identifiers of the regions +regionNames string_array required The list of names of the regions +writeCSV integer 0 Write statistics into a CSV file +================= ============ ======== ============================================== diff --git a/src/coreComponents/schema/schema.xsd b/src/coreComponents/schema/schema.xsd index 937ee56113d..b8d04b7f0e3 100644 --- a/src/coreComponents/schema/schema.xsd +++ b/src/coreComponents/schema/schema.xsd @@ -3660,7 +3660,7 @@ Local - Add stabilization only to interiors of macro elements.--> - + From a69de3279253fa2cb4d88c9dfbfeac381937114b Mon Sep 17 00:00:00 2001 From: dkachuma Date: Thu, 14 Mar 2024 16:53:50 -0500 Subject: [PATCH 5/7] Use element-wise array reduction --- .../fluidFlow/RegionMultiphaseStatistics.cpp | 752 ++++++++++++------ .../fluidFlow/RegionMultiphaseStatistics.hpp | 41 +- 2 files changed, 533 insertions(+), 260 deletions(-) diff --git a/src/coreComponents/physicsSolvers/fluidFlow/RegionMultiphaseStatistics.cpp b/src/coreComponents/physicsSolvers/fluidFlow/RegionMultiphaseStatistics.cpp index bad58bcd56b..cf6c2390380 100644 --- a/src/coreComponents/physicsSolvers/fluidFlow/RegionMultiphaseStatistics.cpp +++ b/src/coreComponents/physicsSolvers/fluidFlow/RegionMultiphaseStatistics.cpp @@ -29,94 +29,124 @@ namespace geos { -using namespace constitutive; -using namespace dataRepository; - -// Container for region statistics -struct RegionMultiphaseStatistics::RegionStatistics +struct RegionMultiphaseStatistics::Statistics { - /// Indices of minimum, maximum and total property values - static constexpr integer MINIMUM = 0; - static constexpr integer MAXIMUM = 1; - static constexpr integer TOTAL = 2; - static constexpr integer TYPE_END = 3; - - /// Indices for property type - static constexpr integer PORE_VOLUME = 0; - static constexpr integer PRESSURE = 1; - static constexpr integer TEMPERATURE = 2; - static constexpr integer PHASE_PORE_VOLUME = 3; - static constexpr integer PHASE_MASS = 4; - static constexpr integer PHASE_COMP_MASS = 5; - - static localIndex allocate( array3d< real64 > & properties, integer regionCount, integer phaseCount, integer componentCount ) + static constexpr integer MAX_NUM_PHASE = 3; + static constexpr integer MAX_NUM_COMPS = 5; + /// Indices for properties + static constexpr integer STATIC_PORE_VOLUME = 0; + static constexpr integer PORE_VOLUME = 1; + + static constexpr integer MINIMUM_PRESSURE = 2; + static constexpr integer MAXIMUM_PRESSURE = 3; + static constexpr integer AVERAGE_PRESSURE = 4; + static constexpr integer MINIMUM_TEMPERATURE = 5; + static constexpr integer MAXIMUM_TEMPERATURE = 6; + static constexpr integer AVERAGE_TEMPERATURE = 7; + // These are starting positions + static constexpr integer PHASE_PORE_VOLUME = 8; + static constexpr integer PHASE_DENSITY = PHASE_PORE_VOLUME + MAX_NUM_PHASE; + static constexpr integer PHASE_VISCOSITY = PHASE_DENSITY + MAX_NUM_PHASE; + static constexpr integer PHASE_MASS = PHASE_VISCOSITY + MAX_NUM_PHASE; + static constexpr integer PHASE_COMP_MASS = PHASE_MASS + MAX_NUM_PHASE; + // End marker + static constexpr integer END = PHASE_COMP_MASS + MAX_NUM_COMPS*MAX_NUM_PHASE; + + GEOS_HOST_DEVICE Statistics( real64 const val ) { - integer const sizes[3] = { TYPE_END, regionCount, getSize( phaseCount, componentCount ) }; - properties.resize( 3, sizes ); - - LvArray::forValuesInSlice( properties.toSlice(), []( real64 & v ){ v = 0.0; } ); - - static constexpr real64 minValue = -LvArray::NumericLimits< real64 >::max; - static constexpr real64 maxValue = LvArray::NumericLimits< real64 >::max; - LvArray::forValuesInSlice( properties[MINIMUM], []( real64 & v ){ v = maxValue; } ); - LvArray::forValuesInSlice( properties[MAXIMUM], []( real64 & v ){ v = minValue; } ); - - return properties.size(); + for( integer i = 0; i < END; i++ ) + { + m_data[i] = val; + } } - static localIndex allocateCounts( array1d< localIndex > & elementCounts, integer regionCount ) + GEOS_HOST_DEVICE Statistics( Statistics const & rhs ) { - elementCounts.resize( regionCount ); - LvArray::forValuesInSlice( elementCounts.toSlice(), []( localIndex & v ){ v = 0; } ); - return elementCounts.size(); + for( integer i = 0; i < END; i++ ) + { + m_data[i] = rhs.m_data[i]; + } } - static real64 * getProperty( arrayView3d< real64 > props, - integer region, - integer phaseCount, - integer componentCount, - integer type, - integer property, - integer phase = 0, - integer component = 0 ) + GEOS_HOST_DEVICE Statistics const & operator=( Statistics const & rhs ) { - integer const ei = calculateIndex( phaseCount, componentCount, property, phase, component ); - return &props( type, region, ei ); + for( integer i = 0; i < END; i++ ) + { + m_data[i] = rhs.m_data[i]; + } + return *this; } -private: - static integer getSize( integer phaseCount, integer componentCount ) + Statistics( Statistics && ) = delete; + ~Statistics() = default; + + GEOS_HOST_DEVICE bool operator!=( Statistics const & ) const { - // 1 pore volume, 1 pressure, 1 temperature, np phase volume, np phase mass, np*nc phase comp mass - return 3 + 2*phaseCount + phaseCount*componentCount; + return false; } - static integer calculateIndex( integer phaseCount, - integer componentCount, - integer property, - integer phase, - integer component ) + real64 m_data[END]; +}; + +} + +namespace RAJA +{ +namespace operators +{ +template< > +struct plus< geos::RegionMultiphaseStatistics::Statistics > +{ + using DataType = geos::RegionMultiphaseStatistics::Statistics; + static constexpr int END = geos::RegionMultiphaseStatistics::Statistics::END; + + GEOS_HOST_DEVICE DataType operator()( const DataType & lhs, const DataType & rhs ) const { - if( property == PORE_VOLUME || property == PRESSURE || property == TEMPERATURE ) + DataType res( lhs ); + for( int i = 0; i < END; ++i ) { - return property; + res.m_data[i] += rhs.m_data[i]; } - else if( property == PHASE_PORE_VOLUME ) - { - return PHASE_PORE_VOLUME + phase; - } - else if( property == PHASE_MASS ) - { - return PHASE_PORE_VOLUME + phaseCount + phase; - } - else if( property == PHASE_COMP_MASS ) + return res; + } + + GEOS_HOST_DEVICE static DataType identity() + { + return DataType( 0.0 ); + } +}; + +template<> +struct minimum< geos::RegionMultiphaseStatistics::Statistics > +{ + using DataType = geos::RegionMultiphaseStatistics::Statistics; + static constexpr int END = geos::RegionMultiphaseStatistics::Statistics::END; + + GEOS_HOST_DEVICE DataType operator()( const DataType & lhs, const DataType & rhs ) const + { + DataType res( lhs ); + for( int i = 0; i < END; ++i ) { - return PHASE_PORE_VOLUME + 2*phaseCount + phase*componentCount + component; + res.m_data[i] = LvArray::math::min( res.m_data[i], rhs.m_data[i] ); } - return 0; + return res; + } + + GEOS_HOST_DEVICE static DataType identity() + { + return DataType( LvArray::NumericLimits< geos::real64 >::max ); } }; +} +} + +namespace geos +{ + +using namespace constitutive; +using namespace dataRepository; + // The region statistics kernel struct RegionMultiphaseStatistics::RegionStatisticsKernel { @@ -124,118 +154,107 @@ struct RegionMultiphaseStatistics::RegionStatisticsKernel static void launch( localIndex const size, integer const phaseCount, integer const componentCount, - arrayView3d< real64 > const regionStatistics, - arrayView1d< localIndex > const & elementCounts, - arrayView1d< real64 const > const & regionIndices, - arrayView1d< real64 const > const & regionMarkers, - arrayView1d< integer const > const & elementGhostRank, + integer const region, + arrayView2d< real64 > const regionStatistics, + arrayView1d< localIndex const > const & targetSet, arrayView1d< real64 const > const & elementVolume, arrayView1d< real64 const > const & pressure, arrayView1d< real64 const > const & temperature, arrayView1d< real64 const > const & referencePorosity, arrayView2d< real64 const > const & porosity, arrayView3d< real64 const, multifluid::USD_PHASE > const & phaseDensity, + arrayView3d< real64 const, multifluid::USD_PHASE > const & phaseViscosity, arrayView4d< real64 const, multifluid::USD_PHASE_COMP > const & phaseComponentFraction, arrayView2d< real64 const, compflow::USD_PHASE > const & phaseVolumeFraction, arrayView3d< real64 const, relperm::USD_RELPERM > const & phaseTrappedVolumeFraction, arrayView3d< real64 const, relperm::USD_RELPERM > const & phaseRelativePermeability ) { - integer const regionCount = regionIndices.size(); - - auto populateMinMaxTotal = [regionStatistics, - phaseCount, - componentCount] - GEOS_HOST_DEVICE ( real64 const value, - real64 const weight, - integer region, - integer property, - integer phase = 0, - integer component = 0 ) - { - real64 * targetValue = nullptr; - // These atomics really make the whole point of threading quite useless here - targetValue = RegionStatistics::getProperty( regionStatistics, region, phaseCount, componentCount, RegionStatistics::MINIMUM, property, phase, component ); - RAJA::atomicMin< AtomicPolicy< POLICY > >( targetValue, value ); - targetValue = RegionStatistics::getProperty( regionStatistics, region, phaseCount, componentCount, RegionStatistics::MAXIMUM, property, phase, component ); - RAJA::atomicMax< AtomicPolicy< POLICY > >( targetValue, value ); - targetValue = RegionStatistics::getProperty( regionStatistics, region, phaseCount, componentCount, RegionStatistics::TOTAL, property, phase, component ); - RAJA::atomicAdd< AtomicPolicy< POLICY > >( targetValue, weight * value ); - return value; - }; + + GEOS_UNUSED_VAR( region ); + + RAJA::ReduceSum< ReducePolicy< POLICY >, Statistics > sumStats( RAJA::operators::plus< Statistics >::identity() ); + RAJA::ReduceMin< ReducePolicy< POLICY >, Statistics > minStats( RAJA::operators::minimum< Statistics >::identity() ); forAll< POLICY >( size, [phaseCount, componentCount, - regionCount, - regionStatistics, - regionIndices, - regionMarkers, - elementCounts, - elementGhostRank, + targetSet, elementVolume, pressure, temperature, referencePorosity, porosity, phaseDensity, + phaseViscosity, phaseComponentFraction, phaseVolumeFraction, phaseTrappedVolumeFraction, phaseRelativePermeability, - &populateMinMaxTotal] GEOS_HOST_DEVICE ( localIndex const ei ) + sumStats, + minStats] GEOS_HOST_DEVICE ( localIndex const i ) { - if( elementGhostRank[ei] >= 0 ) - { - return; - } + localIndex const ei = targetSet[i]; - // Find the appropriate region - integer region = 0; - real64 const regionMarker = regionMarkers[ei]; - for(; region < regionCount; ++region ) - { - if( LvArray::math::abs( regionIndices[region] - regionMarker ) < 0.1 ) - { - break; - } - } - - if( regionCount <= region ) - { - return; - } - - RAJA::atomicAdd< AtomicPolicy< POLICY > >( &elementCounts[region], 1 ); + RegionMultiphaseStatistics::Statistics stats( RAJA::operators::plus< Statistics >::identity() ); // Uncompacted pore volume which is also used as the weight - real64 const weight = elementVolume[ei] * referencePorosity[ei]; + real64 const staticPoreVolume = elementVolume[ei] * referencePorosity[ei]; + real64 const poreVolume = elementVolume[ei] * porosity[ei][0]; - populateMinMaxTotal( weight, 1.0, region, RegionStatistics::PORE_VOLUME ); + stats.m_data[Statistics::STATIC_PORE_VOLUME] = staticPoreVolume; + stats.m_data[Statistics::PORE_VOLUME] = poreVolume; - populateMinMaxTotal( pressure[ei], weight, region, RegionStatistics::PRESSURE ); - populateMinMaxTotal( temperature[ei], weight, region, RegionStatistics::TEMPERATURE ); + stats.m_data[Statistics::MINIMUM_PRESSURE] = pressure[ei]; + stats.m_data[Statistics::MAXIMUM_PRESSURE] = -pressure[ei]; + stats.m_data[Statistics::AVERAGE_PRESSURE] = staticPoreVolume * pressure[ei]; + + stats.m_data[Statistics::MINIMUM_TEMPERATURE] = temperature[ei]; + stats.m_data[Statistics::MAXIMUM_TEMPERATURE] = -temperature[ei]; + stats.m_data[Statistics::AVERAGE_TEMPERATURE] = staticPoreVolume * temperature[ei]; for( integer phaseIndex = 0; phaseIndex < phaseCount; ++phaseIndex ) { - real64 const elementPhaseVolume = elementVolume[ei] * porosity[ei][0] * phaseVolumeFraction[ei][phaseIndex]; + real64 const elementPhaseVolume = poreVolume * phaseVolumeFraction[ei][phaseIndex]; real64 const elementPhaseDensity = phaseDensity[ei][0][phaseIndex]; + real64 const elementPhaseViscosity = phaseViscosity[ei][0][phaseIndex]; real64 const elementPhaseMass = elementPhaseVolume * elementPhaseDensity; - populateMinMaxTotal( elementPhaseVolume, 1.0, region, RegionStatistics::PHASE_PORE_VOLUME, phaseIndex ); - populateMinMaxTotal( elementPhaseMass, 1.0, region, RegionStatistics::PHASE_MASS, phaseIndex ); + stats.m_data[Statistics::PHASE_PORE_VOLUME + phaseIndex] = elementPhaseVolume; + stats.m_data[Statistics::PHASE_MASS + phaseIndex] = elementPhaseMass; + + stats.m_data[Statistics::PHASE_DENSITY + phaseIndex] = elementPhaseVolume * elementPhaseDensity; + stats.m_data[Statistics::PHASE_VISCOSITY + phaseIndex] = elementPhaseVolume * elementPhaseViscosity; + integer const phaseMassIndex = Statistics::PHASE_COMP_MASS + phaseIndex*Statistics::MAX_NUM_COMPS; for( integer compIndex = 0; compIndex < componentCount; ++compIndex ) { real64 const elementComponentPhaseMass = elementPhaseMass * phaseComponentFraction[ei][0][phaseIndex][compIndex]; - populateMinMaxTotal( elementComponentPhaseMass, 1.0, region, RegionStatistics::PHASE_COMP_MASS, phaseIndex, compIndex ); + stats.m_data[phaseMassIndex + compIndex] = elementComponentPhaseMass; } } + + // Atomic operations + sumStats += stats; + minStats.min( stats ); } ); - // Dummy loop to bring data back to the CPU - forAll< serialPolicy >( 1, [regionStatistics, elementCounts] ( localIndex const ) + auto const sumData = sumStats.get().m_data; + auto const minData = minStats.get().m_data; + + arraySlice1d< real64 > regionData = regionStatistics[region]; + + regionData[Statistics::MINIMUM_PRESSURE] = LvArray::math::min( regionData[Statistics::MINIMUM_PRESSURE], minData[Statistics::MINIMUM_PRESSURE] ); + regionData[Statistics::MAXIMUM_PRESSURE] = LvArray::math::max( regionData[Statistics::MAXIMUM_PRESSURE], -minData[Statistics::MAXIMUM_PRESSURE] ); + regionData[Statistics::MINIMUM_TEMPERATURE] = LvArray::math::min( regionData[Statistics::MINIMUM_TEMPERATURE], minData[Statistics::MINIMUM_TEMPERATURE] ); + regionData[Statistics::MAXIMUM_TEMPERATURE] = LvArray::math::max( regionData[Statistics::MAXIMUM_TEMPERATURE], -minData[Statistics::MAXIMUM_TEMPERATURE] ); + + regionData[Statistics::STATIC_PORE_VOLUME] += sumData[Statistics::STATIC_PORE_VOLUME]; + regionData[Statistics::PORE_VOLUME] += sumData[Statistics::PORE_VOLUME]; + regionData[Statistics::AVERAGE_PRESSURE] += sumData[Statistics::AVERAGE_PRESSURE]; + regionData[Statistics::AVERAGE_TEMPERATURE] += sumData[Statistics::AVERAGE_TEMPERATURE]; + for( integer propIndex = Statistics::PHASE_PORE_VOLUME; propIndex < Statistics::END; ++propIndex ) { - GEOS_UNUSED_VAR( regionStatistics, elementCounts ); - } ); + regionData[propIndex] += sumData[propIndex]; + } } }; @@ -272,9 +291,15 @@ void RegionMultiphaseStatistics::postProcessInput() { m_propertyNameTypes.emplace_back( PropertyNameType::Pressure ); m_propertyNameTypes.emplace_back( PropertyNameType::Temperature ); + m_propertyNameTypes.emplace_back( PropertyNameType::StaticPoreVolume ); m_propertyNameTypes.emplace_back( PropertyNameType::PoreVolume ); - m_propertyNameTypes.emplace_back( PropertyNameType::VolumeFraction ); - m_propertyNameTypes.emplace_back( PropertyNameType::Mass ); + m_propertyNameTypes.emplace_back( PropertyNameType::PhaseVolumeFraction ); + m_propertyNameTypes.emplace_back( PropertyNameType::PhasePoreVolume ); + m_propertyNameTypes.emplace_back( PropertyNameType::PhaseMass ); + m_propertyNameTypes.emplace_back( PropertyNameType::PhaseDensity ); + m_propertyNameTypes.emplace_back( PropertyNameType::PhaseViscosity ); + m_propertyNameTypes.emplace_back( PropertyNameType::PhaseComponentMass ); + m_propertyNameTypes.emplace_back( PropertyNameType::ComponentMass ); } else { @@ -286,6 +311,22 @@ void RegionMultiphaseStatistics::postProcessInput() } } +template< typename LAMBDA > +void RegionMultiphaseStatistics::forRegions( Group & meshBodies, LAMBDA && lambda ) +{ + m_solver->forDiscretizationOnMeshTargets( meshBodies, [&] ( string const &, + MeshLevel & mesh, + arrayView1d< string const > const & regionNames ) + { + mesh.getElemManager().forElementSubRegions( regionNames, [&]( localIndex const regionIndex, + ElementSubRegionBase & subRegion ) + { + string const regionName = regionNames[regionIndex]; + lambda( regionIndex, regionName, subRegion ); + } ); + } ); +} + template< typename ARRAY > std::ostream & RegionMultiphaseStatistics::writeArray( std::ostream & os, ARRAY const & array ) const { @@ -317,19 +358,26 @@ void RegionMultiphaseStatistics::registerDataOnMesh( Group & meshBodies ) string const fieldName = GEOS_FMT( "{}{}", getName(), viewKeyStruct::fieldNameString() ); - m_solver->forDiscretizationOnMeshTargets( meshBodies, [&] ( string const &, - MeshLevel & mesh, - arrayView1d< string const > const & regionNames ) + forRegions( meshBodies, [&]( localIndex const, + string const & regionName, + ElementSubRegionBase & subRegion ) { - mesh.getElemManager().forElementSubRegions( regionNames, [&]( localIndex const, - ElementSubRegionBase & subRegion ) + subRegion.registerWrapper< array1d< real64 > >( fieldName ). + setPlotLevel( PlotLevel::NOPLOT ). + setRestartFlags( RestartFlags::NO_WRITE ). + setDescription( fieldName ). + setRegisteringObjects( getName() ); + + for( string const name : m_regionNames ) { - subRegion.registerWrapper< array1d< real64 > >( fieldName ). + string const wrapperName = GEOS_FMT( "{}_{}_{}_{}", getName(), name, regionName, subRegion.getName() ); + + registerWrapper< array1d< localIndex > >( wrapperName ). setPlotLevel( PlotLevel::NOPLOT ). setRestartFlags( RestartFlags::NO_WRITE ). - setDescription( fieldName ). - setRegisteringObjects( getName() ); - } ); + setDescription( wrapperName ). + setSizedFromParent( 0 ); + } } ); // Create the csv file if requires @@ -339,6 +387,71 @@ void RegionMultiphaseStatistics::registerDataOnMesh( Group & meshBodies ) } } +void RegionMultiphaseStatistics::initializePostInitialConditionsPreSubGroups() +{ + string const fieldName = GEOS_FMT( "{}{}", getName(), viewKeyStruct::fieldNameString() ); + + integer const regionCount = m_regionNames.size(); + array1d< integer > regionIdentifiers( regionCount ); + forAll< serialPolicy >( m_regionIdentifiers.size(), [&] ( localIndex const ei ) + { + regionIdentifiers[ei] = LvArray::math::convert< integer >( m_regionIdentifiers[ei] + 0.5 ); + } ); + + DomainPartition & domain = getGroupByPath< DomainPartition >( "/Problem/domain" ); + + forRegions( domain.getMeshBodies(), [&]( localIndex const, + string const & solverRegionName, + ElementSubRegionBase & subRegion ) + { + arrayView1d< real64 const > const regionMarkers = subRegion.getWrapper< array1d< real64 > >( fieldName ) + .reference().toViewConst(); + arrayView1d< integer const > const elementGhostRank = subRegion.ghostRank(); + + array1d< localIndex > regionElementCounts( regionCount ); + regionElementCounts.zero(); + + array2d< localIndex > regionElements( regionCount, regionMarkers.size() ); + + forAll< serialPolicy >( regionMarkers.size(), [&] ( localIndex const ei ) + { + if( 0 <= elementGhostRank[ei] ) + { + return; + } + integer const marker = LvArray::math::convert< integer >( regionMarkers[ei] + 0.5 ); + for( integer region = 0; region < regionCount; region++ ) + { + if( regionIdentifiers[region] == marker ) + { + regionElements( region, regionElementCounts[region]++ ) = ei; + return; + } + } + } ); + + arrayView2d< localIndex const > const regionElementsView = regionElements.toViewConst(); + + for( integer region = 0; region < regionCount; region++ ) + { + string const wrapperName = GEOS_FMT( "{}_{}_{}_{}", getName(), m_regionNames[region], solverRegionName, subRegion.getName() ); + auto & regionWrapper = getWrapper< array1d< localIndex > >( wrapperName ); + regionWrapper.resize( regionElementCounts[region] ); + + arrayView1d< localIndex > const regionArrayView = regionWrapper.reference().toView(); + + forAll< parallelDevicePolicy<> >( regionElementCounts[region], + [region, + regionArrayView, + regionElementsView] + GEOS_HOST_DEVICE ( localIndex const ei ) + { + regionArrayView[ei] = regionElementsView( region, ei ); + } ); + } + } ); +} + bool RegionMultiphaseStatistics::execute( real64 const time_n, real64 const dt, integer const GEOS_UNUSED_PARAM( cycleNumber ), @@ -369,25 +482,32 @@ void RegionMultiphaseStatistics::computeRegionStatistics( real64 const time, { GEOS_MARK_FUNCTION; + GEOS_UNUSED_VAR( time ); + integer const phaseCount = m_solver->numFluidPhases(); integer const componentCount = m_solver->numFluidComponents(); + integer const regionCount = m_regionNames.size(); string const fieldName = GEOS_FMT( "{}{}", getName(), viewKeyStruct::fieldNameString() ); // Step 1: initialize the average/min/max quantities // All properties will be stored in one large array - array3d< real64 > regionStatistics; - array1d< localIndex > elementCounts; - RegionStatistics::allocate( regionStatistics, m_regionNames.size(), phaseCount, componentCount ); - RegionStatistics::allocateCounts( elementCounts, m_regionNames.size() ); + array2d< real64 > regionStatisticsData( regionCount, Statistics::END ); + regionStatisticsData.zero(); + arrayView2d< real64 > regionStatistics = regionStatisticsData.toView(); + for( integer region = 0; region < regionCount; ++region ) + { + regionStatistics( region, Statistics::MINIMUM_PRESSURE ) = LvArray::NumericLimits< real64 >::max; + regionStatistics( region, Statistics::MAXIMUM_PRESSURE ) = -LvArray::NumericLimits< real64 >::max; + regionStatistics( region, Statistics::MINIMUM_TEMPERATURE ) = LvArray::NumericLimits< real64 >::max; + regionStatistics( region, Statistics::MAXIMUM_TEMPERATURE ) = -LvArray::NumericLimits< real64 >::max; + } // Step 2: increment the average/min/max quantities for all the subRegions - mesh.getElemManager().forElementSubRegions( regionNames, [&]( localIndex const, + mesh.getElemManager().forElementSubRegions( regionNames, [&]( localIndex const solverRegionIndex, ElementSubRegionBase & subRegion ) { - arrayView1d< integer const > const elementGhostRank = subRegion.ghostRank(); arrayView1d< real64 const > const elementVolume = subRegion.getElementVolume(); - array1d< real64 > const regionMarkers = subRegion.getWrapper< array1d< real64 > >( fieldName ).reference(); arrayView1d< real64 const > const pressure = subRegion.getField< fields::flow::pressure >(); arrayView1d< real64 const > const temperature = subRegion.getField< fields::flow::temperature >(); arrayView2d< real64 const, compflow::USD_PHASE > const phaseVolumeFraction = subRegion.getField< fields::flow::phaseVolumeFraction >(); @@ -402,6 +522,7 @@ void RegionMultiphaseStatistics::computeRegionStatistics( real64 const time, 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(); + arrayView3d< real64 const, multifluid::USD_PHASE > const phaseViscosity = fluid.phaseViscosity(); arrayView4d< real64 const, multifluid::USD_PHASE_COMP > const phaseComponentFraction = fluid.phaseCompFraction(); string const & relpermName = subRegion.getReference< string >( CompositionalMultiphaseBase::viewKeyStruct::relPermNamesString() ); @@ -409,123 +530,84 @@ void RegionMultiphaseStatistics::computeRegionStatistics( real64 const time, arrayView3d< real64 const, relperm::USD_RELPERM > const phaseTrappedVolumeFraction = relperm.phaseTrappedVolFraction(); arrayView3d< real64 const, relperm::USD_RELPERM > const phaseRelativePermeability = relperm.phaseRelPerm(); - RegionStatisticsKernel::launch< parallelDevicePolicy<> >( subRegion.size(), - phaseCount, - componentCount, - regionStatistics.toView(), - elementCounts.toView(), - m_regionIdentifiers.toViewConst(), - regionMarkers.toViewConst(), - elementGhostRank, - elementVolume, - pressure, - temperature, - referencePorosity, - porosity, - phaseDensity, - phaseComponentFraction, - phaseVolumeFraction, - phaseTrappedVolumeFraction, - phaseRelativePermeability ); - } ); - - // Step 3: synchronize the results over the MPI ranks - integer const regionCount = m_regionIdentifiers.size(); - for( integer region = 0; region < regionCount; ++region ) - { - elementCounts[region] = MpiWrapper::sum( elementCounts[region] ); - } - - localIndex const size1 = regionStatistics.size( 1 ); - localIndex const size2 = regionStatistics.size( 2 ); - for( localIndex i1 = 0; i1 < size1; ++i1 ) - { - for( localIndex i2 = 0; i2 < size2; ++i2 ) + for( integer region = 0; region < regionCount; ++region ) { - regionStatistics( RegionStatistics::MINIMUM, i1, i2 ) = MpiWrapper::min( regionStatistics( RegionStatistics::MINIMUM, i1, i2 ) ); - regionStatistics( RegionStatistics::MAXIMUM, i1, i2 ) = MpiWrapper::max( regionStatistics( RegionStatistics::MAXIMUM, i1, i2 ) ); - regionStatistics( RegionStatistics::TOTAL, i1, i2 ) = MpiWrapper::sum( regionStatistics( RegionStatistics::TOTAL, i1, i2 ) ); - } - } + string const wrapperName = GEOS_FMT( "{}_{}_{}_{}", getName(), m_regionNames[region], regionNames[solverRegionIndex], subRegion.getName() ); + arrayView1d< localIndex const > const targetSet = getWrapper< array1d< localIndex > >( wrapperName ) + .reference() + .toViewConst(); + if( targetSet.size() == 0 ) + { + continue; + } - array1d< real64 > phaseVolumes( phaseCount ); - real64 * propValue = nullptr; + RegionStatisticsKernel::launch< parallelDevicePolicy<> >( targetSet.size(), + phaseCount, + componentCount, + region, + regionStatistics, + targetSet, + elementVolume, + pressure, + temperature, + referencePorosity, + porosity, + phaseDensity, + phaseViscosity, + phaseComponentFraction, + phaseVolumeFraction, + phaseTrappedVolumeFraction, + phaseRelativePermeability ); + } - array1d< real64 > propColumns; - propColumns.emplace_back( time ); - for( integer region = 0; region < regionCount; ++region ) - { - // Extract the uncompacted region pore volume as a weight - propValue = RegionStatistics::getProperty( regionStatistics, region, phaseCount, componentCount, RegionStatistics::TOTAL, RegionStatistics::PORE_VOLUME ); - real64 const regionPoreVolume = *propValue; - real64 const inverseRegionPoreVolume = regionPoreVolume < LvArray::NumericLimits< real64 >::epsilon ? 0.0 : 1.0 / regionPoreVolume; + } ); - for( auto const & prop : m_propertyNameTypes ) - { - if( prop == PropertyNameType::Pressure ) - { - propValue = RegionStatistics::getProperty( regionStatistics, region, phaseCount, componentCount, RegionStatistics::MINIMUM, RegionStatistics::PRESSURE ); - propColumns.emplace_back( *propValue ); + // Step 3: Force data back to cpu if necessary + forAll< serialPolicy >( 1, [regionStatistics]( localIndex const ){ + GEOS_UNUSED_VAR( regionStatistics ); + } ); - propValue = RegionStatistics::getProperty( regionStatistics, region, phaseCount, componentCount, RegionStatistics::TOTAL, RegionStatistics::PRESSURE ); - real64 const avgPressure = (*propValue) * inverseRegionPoreVolume; - propColumns.emplace_back( avgPressure ); + // Step 4: synchronize the results over the MPI ranks + array1d< real64 > phaseVolumes( phaseCount ); - propValue = RegionStatistics::getProperty( regionStatistics, region, phaseCount, componentCount, RegionStatistics::MAXIMUM, RegionStatistics::PRESSURE ); - propColumns.emplace_back( *propValue ); - } + for( integer region =0; region regData = regionStatistics[region]; - if( prop == PropertyNameType::Temperature ) - { - propValue = RegionStatistics::getProperty( regionStatistics, region, phaseCount, componentCount, RegionStatistics::MINIMUM, RegionStatistics::TEMPERATURE ); - propColumns.emplace_back( *propValue ); + regData[Statistics::STATIC_PORE_VOLUME] = MpiWrapper::sum( regData[Statistics::STATIC_PORE_VOLUME] ); + regData[Statistics::PORE_VOLUME] = MpiWrapper::sum( regData[Statistics::PORE_VOLUME] ); - propValue = RegionStatistics::getProperty( regionStatistics, region, phaseCount, componentCount, RegionStatistics::TOTAL, RegionStatistics::TEMPERATURE ); - real64 const avgPressure = (*propValue) * inverseRegionPoreVolume; - propColumns.emplace_back( avgPressure ); + regData[Statistics::MINIMUM_PRESSURE] = MpiWrapper::min( regData[Statistics::MINIMUM_PRESSURE] ); + regData[Statistics::MAXIMUM_PRESSURE] = MpiWrapper::max( regData[Statistics::MAXIMUM_PRESSURE] ); + regData[Statistics::MINIMUM_TEMPERATURE] = MpiWrapper::min( regData[Statistics::MINIMUM_TEMPERATURE] ); + regData[Statistics::MAXIMUM_TEMPERATURE] = MpiWrapper::max( regData[Statistics::MAXIMUM_TEMPERATURE] ); - propValue = RegionStatistics::getProperty( regionStatistics, region, phaseCount, componentCount, RegionStatistics::MAXIMUM, RegionStatistics::TEMPERATURE ); - propColumns.emplace_back( *propValue ); - } + regData[Statistics::AVERAGE_PRESSURE] = MpiWrapper::sum( regData[Statistics::AVERAGE_PRESSURE] ); + regData[Statistics::AVERAGE_TEMPERATURE] = MpiWrapper::sum( regData[Statistics::AVERAGE_TEMPERATURE] ); - if( prop == PropertyNameType::PoreVolume ) - { - for( integer phaseIndex = 0; phaseIndex < phaseCount; ++phaseIndex ) - { - propValue = RegionStatistics::getProperty( regionStatistics, region, phaseCount, componentCount, RegionStatistics::TOTAL, RegionStatistics::PHASE_PORE_VOLUME, phaseIndex ); - propColumns.emplace_back( *propValue ); - } - } + for( integer propIndex = Statistics::PHASE_PORE_VOLUME; propIndex < Statistics::END; ++propIndex ) + { + regData[propIndex] = MpiWrapper::sum( regData[propIndex] ); + } - if( prop == PropertyNameType::Mass ) - { - for( integer phaseIndex = 0; phaseIndex < phaseCount; ++phaseIndex ) - { - propValue = RegionStatistics::getProperty( regionStatistics, region, phaseCount, componentCount, RegionStatistics::TOTAL, RegionStatistics::PHASE_MASS, phaseIndex ); - propColumns.emplace_back( *propValue ); - } - } + // Actually calculate averages + real64 const invStaticPoreVolume = safeInverse( regData[Statistics::STATIC_PORE_VOLUME] ); + regData[Statistics::AVERAGE_PRESSURE] *= invStaticPoreVolume; + regData[Statistics::AVERAGE_TEMPERATURE] *= invStaticPoreVolume; - if( prop == PropertyNameType::VolumeFraction ) - { - real64 totalPhaseVolume = 0.0; - for( integer phaseIndex = 0; phaseIndex < phaseCount; ++phaseIndex ) - { - propValue = RegionStatistics::getProperty( regionStatistics, region, phaseCount, componentCount, RegionStatistics::TOTAL, RegionStatistics::PHASE_PORE_VOLUME, phaseIndex ); - phaseVolumes[phaseIndex] = *propValue; - totalPhaseVolume += phaseVolumes[phaseIndex]; - } - real64 const inverseTotalPhaseVolume = totalPhaseVolume < LvArray::NumericLimits< real64 >::epsilon ? 0.0 : 1.0 / totalPhaseVolume; - for( integer phaseIndex = 0; phaseIndex < phaseCount; ++phaseIndex ) - { - propColumns.emplace_back( phaseVolumes[phaseIndex] * inverseTotalPhaseVolume ); - } - } + // Calculate average density and viscosity + for( integer phaseIndex = 0; phaseIndex < phaseCount; ++phaseIndex ) + { + real64 const invPhaseVolume = safeInverse( regData[Statistics::PHASE_PORE_VOLUME + phaseIndex] ); + regData[Statistics::PHASE_DENSITY + phaseIndex] *= invPhaseVolume; + regData[Statistics::PHASE_VISCOSITY + phaseIndex] *= invPhaseVolume; } } if( m_writeCSV > 0 && MpiWrapper::commRank() == 0 ) { + array1d< real64 > propColumns; + populateColumns( time, propColumns, regionStatistics ); std::ofstream outputFile( m_outputDir + "/" + viewKeyStruct::fileNameString() + ".csv", std::ios_base::app ); writeArray( outputFile, propColumns ); outputFile.close(); @@ -543,6 +625,7 @@ void RegionMultiphaseStatistics::initializeFile() const auto fluidPhaseNames = fluidModel.phaseNames(); integer const phaseCount = m_solver->numFluidPhases(); + integer const componentCount = m_solver->numFluidComponents(); integer const regionCount = m_regionIdentifiers.size(); auto regionNames = m_regionNames.toView(); @@ -556,7 +639,8 @@ void RegionMultiphaseStatistics::initializeFile() const string const propUnit, integer const regionIndex=-1, integer const phaseIndex=-1, - integer const compIndex=-1 ){ + integer const compIndex=-1 ) + { propNames.emplace_back( propName ); propUnits.emplace_back( propUnit ); if( 0 <= regionIndex ) @@ -587,6 +671,7 @@ void RegionMultiphaseStatistics::initializeFile() const integer const useMass = m_solver->getReference< integer >( CompositionalMultiphaseBase::viewKeyStruct::useMassFlagString() ); string const massUnit = useMass ? "kg" : "mol"; + string const densityUnit = useMass ? "kg/m^3" : "mol/m^3"; addColumn( "Time", "s" ); for( integer region = 0; region < regionCount; ++region ) @@ -611,7 +696,27 @@ void RegionMultiphaseStatistics::initializeFile() const addColumn( GEOS_FMT( "Max {}", propName ), propUnit, region ); } + if( prop == PropertyNameType::StaticPoreVolume ) + { + propUnit = "rm^3"; + addColumn( propName, propUnit, region ); + } + if( prop == PropertyNameType::PoreVolume ) + { + propUnit = "rm^3"; + addColumn( propName, propUnit, region ); + } + + if( prop == PropertyNameType::PhaseVolumeFraction ) + { + for( integer phaseIndex = 0; phaseIndex < phaseCount; ++phaseIndex ) + { + addColumn( propName, propUnit, region, phaseIndex ); + } + } + + if( prop == PropertyNameType::PhasePoreVolume ) { propUnit = "rm^3"; for( integer phaseIndex = 0; phaseIndex < phaseCount; ++phaseIndex ) @@ -620,7 +725,7 @@ void RegionMultiphaseStatistics::initializeFile() const } } - if( prop == PropertyNameType::Mass ) + if( prop == PropertyNameType::PhaseMass ) { for( integer phaseIndex = 0; phaseIndex < phaseCount; ++phaseIndex ) { @@ -628,13 +733,41 @@ void RegionMultiphaseStatistics::initializeFile() const } } - if( prop == PropertyNameType::VolumeFraction ) + if( prop == PropertyNameType::PhaseDensity ) { + for( integer phaseIndex = 0; phaseIndex < phaseCount; ++phaseIndex ) + { + addColumn( propName, densityUnit, region, phaseIndex ); + } + } + + if( prop == PropertyNameType::PhaseViscosity ) + { + propUnit = "Pa s"; for( integer phaseIndex = 0; phaseIndex < phaseCount; ++phaseIndex ) { addColumn( propName, propUnit, region, phaseIndex ); } } + + if( prop == PropertyNameType::PhaseComponentMass ) + { + for( integer phaseIndex = 0; phaseIndex < phaseCount; ++phaseIndex ) + { + for( integer compIndex = 0; compIndex < componentCount; ++compIndex ) + { + addColumn( propName, massUnit, region, phaseIndex, compIndex ); + } + } + } + + if( prop == PropertyNameType::ComponentMass ) + { + for( integer compIndex = 0; compIndex < componentCount; ++compIndex ) + { + addColumn( propName, massUnit, region, -1, compIndex ); + } + } } } @@ -649,6 +782,115 @@ void RegionMultiphaseStatistics::initializeFile() const outputFile.close(); } +void RegionMultiphaseStatistics::populateColumns( real64 const time, + array1d< real64 > & columns, + arrayView2d< real64 const > const & regionStatistics ) const +{ + integer const phaseCount = m_solver->numFluidPhases(); + integer const componentCount = m_solver->numFluidComponents(); + integer const regionCount = m_regionIdentifiers.size(); + + columns.clear(); + + columns.emplace_back( time ); + for( integer region = 0; region < regionCount; ++region ) + { + for( auto const & prop : m_propertyNameTypes ) + { + if( prop == PropertyNameType::Pressure ) + { + columns.emplace_back( regionStatistics( region, Statistics::MINIMUM_PRESSURE ) ); + columns.emplace_back( regionStatistics( region, Statistics::AVERAGE_PRESSURE ) ); + columns.emplace_back( regionStatistics( region, Statistics::MAXIMUM_PRESSURE ) ); + } + + if( prop == PropertyNameType::Temperature ) + { + columns.emplace_back( regionStatistics( region, Statistics::MINIMUM_TEMPERATURE ) ); + columns.emplace_back( regionStatistics( region, Statistics::AVERAGE_TEMPERATURE ) ); + columns.emplace_back( regionStatistics( region, Statistics::MAXIMUM_TEMPERATURE ) ); + } + + if( prop == PropertyNameType::StaticPoreVolume ) + { + columns.emplace_back( regionStatistics( region, Statistics::STATIC_PORE_VOLUME ) ); + } + + if( prop == PropertyNameType::PoreVolume ) + { + columns.emplace_back( regionStatistics( region, Statistics::PORE_VOLUME ) ); + } + + if( prop == PropertyNameType::PhaseVolumeFraction ) + { + real64 const invPoreVolume = safeInverse( regionStatistics( region, Statistics::PORE_VOLUME )); + for( integer phaseIndex = 0; phaseIndex < phaseCount; ++phaseIndex ) + { + columns.emplace_back( invPoreVolume * regionStatistics( region, Statistics::PHASE_PORE_VOLUME + phaseIndex ) ); + } + } + + if( prop == PropertyNameType::PhasePoreVolume ) + { + for( integer phaseIndex = 0; phaseIndex < phaseCount; ++phaseIndex ) + { + columns.emplace_back( regionStatistics( region, Statistics::PHASE_PORE_VOLUME + phaseIndex ) ); + } + } + + if( prop == PropertyNameType::PhaseMass ) + { + for( integer phaseIndex = 0; phaseIndex < phaseCount; ++phaseIndex ) + { + columns.emplace_back( regionStatistics( region, Statistics::PHASE_MASS + phaseIndex ) ); + } + } + + if( prop == PropertyNameType::PhaseDensity ) + { + for( integer phaseIndex = 0; phaseIndex < phaseCount; ++phaseIndex ) + { + columns.emplace_back( regionStatistics( region, Statistics::PHASE_DENSITY + phaseIndex ) ); + } + } + + if( prop == PropertyNameType::PhaseViscosity ) + { + for( integer phaseIndex = 0; phaseIndex < phaseCount; ++phaseIndex ) + { + columns.emplace_back( regionStatistics( region, Statistics::PHASE_VISCOSITY + phaseIndex ) ); + } + } + + if( prop == PropertyNameType::PhaseComponentMass ) + { + for( integer phaseIndex = 0; phaseIndex < phaseCount; ++phaseIndex ) + { + integer const phaseMassIndex = Statistics::PHASE_COMP_MASS + phaseIndex*Statistics::MAX_NUM_COMPS; + for( integer compIndex = 0; compIndex < componentCount; ++compIndex ) + { + columns.emplace_back( regionStatistics( region, phaseMassIndex + compIndex ) ); + } + } + } + + if( prop == PropertyNameType::ComponentMass ) + { + for( integer compIndex = 0; compIndex < componentCount; ++compIndex ) + { + real64 componentMass = 0.0; + for( integer phaseIndex = 0; phaseIndex < phaseCount; ++phaseIndex ) + { + integer const phaseMassIndex = Statistics::PHASE_COMP_MASS + phaseIndex*Statistics::MAX_NUM_COMPS; + componentMass += regionStatistics( region, phaseMassIndex + compIndex ); + } + columns.emplace_back( componentMass ); + } + } + } + } +} + REGISTER_CATALOG_ENTRY( TaskBase, RegionMultiphaseStatistics, string const &, diff --git a/src/coreComponents/physicsSolvers/fluidFlow/RegionMultiphaseStatistics.hpp b/src/coreComponents/physicsSolvers/fluidFlow/RegionMultiphaseStatistics.hpp index 60322442ed7..2c7164cbee0 100644 --- a/src/coreComponents/physicsSolvers/fluidFlow/RegionMultiphaseStatistics.hpp +++ b/src/coreComponents/physicsSolvers/fluidFlow/RegionMultiphaseStatistics.hpp @@ -40,11 +40,20 @@ class RegionMultiphaseStatistics : public FieldStatisticsBase< CompositionalMult { Pressure, Temperature, + StaticPoreVolume, PoreVolume, - VolumeFraction, - Mass + PhaseVolumeFraction, + PhasePoreVolume, + PhaseMass, + PhaseDensity, + PhaseViscosity, + PhaseComponentMass, + ComponentMass, }; +public: + struct Statistics; + public: /** @@ -78,6 +87,8 @@ class RegionMultiphaseStatistics : public FieldStatisticsBase< CompositionalMult using Base = FieldStatisticsBase< CompositionalMultiphaseBase >; + static constexpr real64 hugeValue = LvArray::NumericLimits< real64 >::max; + /** * @struct viewKeyStruct holds char strings and viewKeys for fast lookup */ @@ -99,6 +110,11 @@ class RegionMultiphaseStatistics : public FieldStatisticsBase< CompositionalMult void registerDataOnMesh( Group & meshBodies ) override; + void initializePostInitialConditionsPreSubGroups() override; + + template< typename LAMBDA > + void forRegions( Group & meshBodies, LAMBDA && lambda ); + /** * @brief Compute some statistics on the regions * @param[in] time current time @@ -111,10 +127,19 @@ class RegionMultiphaseStatistics : public FieldStatisticsBase< CompositionalMult void initializeFile() const; + void populateColumns( real64 const time, + array1d< real64 > & columns, + arrayView2d< real64 const > const & regionStatistics ) const; + template< typename ARRAY > std::ostream & writeArray( std::ostream & os, ARRAY const & array ) const; - struct RegionStatistics; + GEOS_HOST_DEVICE + static constexpr real64 safeInverse( real64 const val ) + { + return (val < LvArray::NumericLimits< real64 >::epsilon) ? 0.0 : 1.0 / val; + } + struct RegionStatisticsKernel; // The names of regions @@ -133,9 +158,15 @@ class RegionMultiphaseStatistics : public FieldStatisticsBase< CompositionalMult ENUM_STRINGS( RegionMultiphaseStatistics::PropertyNameType, "pressure", "temperature", + "staticPoreVolume", "poreVolume", - "volumeFraction", - "mass" ); + "phaseVolumeFraction", + "phasePoreVolume", + "phaseMass", + "phaseDensity", + "phaseViscosity", + "phaseComponentMass", + "componentMass" ); } /* namespace geos */ From be2bd30facf411ffd9637db07a1bcdffa2713909 Mon Sep 17 00:00:00 2001 From: dkachuma Date: Thu, 14 Mar 2024 18:07:00 -0500 Subject: [PATCH 6/7] Use cell count to average density and viscosity --- .../fluidFlow/RegionMultiphaseStatistics.cpp | 75 ++++++++++--------- 1 file changed, 39 insertions(+), 36 deletions(-) diff --git a/src/coreComponents/physicsSolvers/fluidFlow/RegionMultiphaseStatistics.cpp b/src/coreComponents/physicsSolvers/fluidFlow/RegionMultiphaseStatistics.cpp index cf6c2390380..0cc2e82ffaa 100644 --- a/src/coreComponents/physicsSolvers/fluidFlow/RegionMultiphaseStatistics.cpp +++ b/src/coreComponents/physicsSolvers/fluidFlow/RegionMultiphaseStatistics.cpp @@ -34,17 +34,18 @@ struct RegionMultiphaseStatistics::Statistics static constexpr integer MAX_NUM_PHASE = 3; static constexpr integer MAX_NUM_COMPS = 5; /// Indices for properties - static constexpr integer STATIC_PORE_VOLUME = 0; - static constexpr integer PORE_VOLUME = 1; - - static constexpr integer MINIMUM_PRESSURE = 2; - static constexpr integer MAXIMUM_PRESSURE = 3; - static constexpr integer AVERAGE_PRESSURE = 4; - static constexpr integer MINIMUM_TEMPERATURE = 5; - static constexpr integer MAXIMUM_TEMPERATURE = 6; - static constexpr integer AVERAGE_TEMPERATURE = 7; + static constexpr integer CELL_COUNT = 0; + static constexpr integer STATIC_PORE_VOLUME = 1; + static constexpr integer PORE_VOLUME = 2; + + static constexpr integer MINIMUM_PRESSURE = 3; + static constexpr integer MAXIMUM_PRESSURE = 4; + static constexpr integer AVERAGE_PRESSURE = 5; + static constexpr integer MINIMUM_TEMPERATURE = 6; + static constexpr integer MAXIMUM_TEMPERATURE = 7; + static constexpr integer AVERAGE_TEMPERATURE = 8; // These are starting positions - static constexpr integer PHASE_PORE_VOLUME = 8; + static constexpr integer PHASE_PORE_VOLUME = 9; static constexpr integer PHASE_DENSITY = PHASE_PORE_VOLUME + MAX_NUM_PHASE; static constexpr integer PHASE_VISCOSITY = PHASE_DENSITY + MAX_NUM_PHASE; static constexpr integer PHASE_MASS = PHASE_VISCOSITY + MAX_NUM_PHASE; @@ -200,6 +201,7 @@ struct RegionMultiphaseStatistics::RegionStatisticsKernel real64 const staticPoreVolume = elementVolume[ei] * referencePorosity[ei]; real64 const poreVolume = elementVolume[ei] * porosity[ei][0]; + stats.m_data[Statistics::CELL_COUNT] = 1.0; stats.m_data[Statistics::STATIC_PORE_VOLUME] = staticPoreVolume; stats.m_data[Statistics::PORE_VOLUME] = poreVolume; @@ -221,8 +223,8 @@ struct RegionMultiphaseStatistics::RegionStatisticsKernel stats.m_data[Statistics::PHASE_PORE_VOLUME + phaseIndex] = elementPhaseVolume; stats.m_data[Statistics::PHASE_MASS + phaseIndex] = elementPhaseMass; - stats.m_data[Statistics::PHASE_DENSITY + phaseIndex] = elementPhaseVolume * elementPhaseDensity; - stats.m_data[Statistics::PHASE_VISCOSITY + phaseIndex] = elementPhaseVolume * elementPhaseViscosity; + stats.m_data[Statistics::PHASE_DENSITY + phaseIndex] = elementPhaseDensity; + stats.m_data[Statistics::PHASE_VISCOSITY + phaseIndex] = elementPhaseViscosity; integer const phaseMassIndex = Statistics::PHASE_COMP_MASS + phaseIndex*Statistics::MAX_NUM_COMPS; for( integer compIndex = 0; compIndex < componentCount; ++compIndex ) @@ -237,23 +239,27 @@ struct RegionMultiphaseStatistics::RegionStatisticsKernel minStats.min( stats ); } ); - auto const sumData = sumStats.get().m_data; - auto const minData = minStats.get().m_data; + auto const sumStatsData = sumStats.get(); + auto const minStatsData = minStats.get(); + auto const & sumData = sumStatsData.m_data; + auto const & minData = minStatsData.m_data; - arraySlice1d< real64 > regionData = regionStatistics[region]; + arraySlice1d< real64 > regData = regionStatistics[region]; + + regData[Statistics::CELL_COUNT] += size; + regData[Statistics::STATIC_PORE_VOLUME] += sumData[Statistics::STATIC_PORE_VOLUME]; + regData[Statistics::PORE_VOLUME] += sumData[Statistics::PORE_VOLUME]; - regionData[Statistics::MINIMUM_PRESSURE] = LvArray::math::min( regionData[Statistics::MINIMUM_PRESSURE], minData[Statistics::MINIMUM_PRESSURE] ); - regionData[Statistics::MAXIMUM_PRESSURE] = LvArray::math::max( regionData[Statistics::MAXIMUM_PRESSURE], -minData[Statistics::MAXIMUM_PRESSURE] ); - regionData[Statistics::MINIMUM_TEMPERATURE] = LvArray::math::min( regionData[Statistics::MINIMUM_TEMPERATURE], minData[Statistics::MINIMUM_TEMPERATURE] ); - regionData[Statistics::MAXIMUM_TEMPERATURE] = LvArray::math::max( regionData[Statistics::MAXIMUM_TEMPERATURE], -minData[Statistics::MAXIMUM_TEMPERATURE] ); + regData[Statistics::MINIMUM_PRESSURE] = LvArray::math::min( regData[Statistics::MINIMUM_PRESSURE], minData[Statistics::MINIMUM_PRESSURE] ); + regData[Statistics::MAXIMUM_PRESSURE] = LvArray::math::max( regData[Statistics::MAXIMUM_PRESSURE], -minData[Statistics::MAXIMUM_PRESSURE] ); + regData[Statistics::MINIMUM_TEMPERATURE] = LvArray::math::min( regData[Statistics::MINIMUM_TEMPERATURE], minData[Statistics::MINIMUM_TEMPERATURE] ); + regData[Statistics::MAXIMUM_TEMPERATURE] = LvArray::math::max( regData[Statistics::MAXIMUM_TEMPERATURE], -minData[Statistics::MAXIMUM_TEMPERATURE] ); - regionData[Statistics::STATIC_PORE_VOLUME] += sumData[Statistics::STATIC_PORE_VOLUME]; - regionData[Statistics::PORE_VOLUME] += sumData[Statistics::PORE_VOLUME]; - regionData[Statistics::AVERAGE_PRESSURE] += sumData[Statistics::AVERAGE_PRESSURE]; - regionData[Statistics::AVERAGE_TEMPERATURE] += sumData[Statistics::AVERAGE_TEMPERATURE]; + regData[Statistics::AVERAGE_PRESSURE] += sumData[Statistics::AVERAGE_PRESSURE]; + regData[Statistics::AVERAGE_TEMPERATURE] += sumData[Statistics::AVERAGE_TEMPERATURE]; for( integer propIndex = Statistics::PHASE_PORE_VOLUME; propIndex < Statistics::END; ++propIndex ) { - regionData[propIndex] += sumData[propIndex]; + regData[propIndex] += sumData[propIndex]; } } }; @@ -482,8 +488,6 @@ void RegionMultiphaseStatistics::computeRegionStatistics( real64 const time, { GEOS_MARK_FUNCTION; - GEOS_UNUSED_VAR( time ); - integer const phaseCount = m_solver->numFluidPhases(); integer const componentCount = m_solver->numFluidComponents(); integer const regionCount = m_regionNames.size(); @@ -497,10 +501,10 @@ void RegionMultiphaseStatistics::computeRegionStatistics( real64 const time, arrayView2d< real64 > regionStatistics = regionStatisticsData.toView(); for( integer region = 0; region < regionCount; ++region ) { - regionStatistics( region, Statistics::MINIMUM_PRESSURE ) = LvArray::NumericLimits< real64 >::max; - regionStatistics( region, Statistics::MAXIMUM_PRESSURE ) = -LvArray::NumericLimits< real64 >::max; - regionStatistics( region, Statistics::MINIMUM_TEMPERATURE ) = LvArray::NumericLimits< real64 >::max; - regionStatistics( region, Statistics::MAXIMUM_TEMPERATURE ) = -LvArray::NumericLimits< real64 >::max; + regionStatistics( region, Statistics::MINIMUM_PRESSURE ) = hugeValue; + regionStatistics( region, Statistics::MAXIMUM_PRESSURE ) = -hugeValue; + regionStatistics( region, Statistics::MINIMUM_TEMPERATURE ) = hugeValue; + regionStatistics( region, Statistics::MAXIMUM_TEMPERATURE ) = -hugeValue; } // Step 2: increment the average/min/max quantities for all the subRegions @@ -568,12 +572,11 @@ void RegionMultiphaseStatistics::computeRegionStatistics( real64 const time, } ); // Step 4: synchronize the results over the MPI ranks - array1d< real64 > phaseVolumes( phaseCount ); - for( integer region =0; region regData = regionStatistics[region]; + regData[Statistics::CELL_COUNT] = MpiWrapper::sum( regData[Statistics::CELL_COUNT] ); regData[Statistics::STATIC_PORE_VOLUME] = MpiWrapper::sum( regData[Statistics::STATIC_PORE_VOLUME] ); regData[Statistics::PORE_VOLUME] = MpiWrapper::sum( regData[Statistics::PORE_VOLUME] ); @@ -595,12 +598,12 @@ void RegionMultiphaseStatistics::computeRegionStatistics( real64 const time, regData[Statistics::AVERAGE_PRESSURE] *= invStaticPoreVolume; regData[Statistics::AVERAGE_TEMPERATURE] *= invStaticPoreVolume; - // Calculate average density and viscosity + // Calculate average density and viscosity using cell count as weight for( integer phaseIndex = 0; phaseIndex < phaseCount; ++phaseIndex ) { - real64 const invPhaseVolume = safeInverse( regData[Statistics::PHASE_PORE_VOLUME + phaseIndex] ); - regData[Statistics::PHASE_DENSITY + phaseIndex] *= invPhaseVolume; - regData[Statistics::PHASE_VISCOSITY + phaseIndex] *= invPhaseVolume; + real64 const invCellCount = safeInverse( regData[Statistics::CELL_COUNT] ); + regData[Statistics::PHASE_DENSITY + phaseIndex] *= invCellCount; + regData[Statistics::PHASE_VISCOSITY + phaseIndex] *= invCellCount; } } From 4e553bfe6f892a9034acd447f1ba543b8f9d5c16 Mon Sep 17 00:00:00 2001 From: dkachuma Date: Thu, 14 Mar 2024 18:52:22 -0500 Subject: [PATCH 7/7] Implement move constructor --- .../fluidFlow/RegionMultiphaseStatistics.cpp | 6 +++++- .../fluidFlow/RegionMultiphaseStatistics.hpp | 7 +++---- 2 files changed, 8 insertions(+), 5 deletions(-) diff --git a/src/coreComponents/physicsSolvers/fluidFlow/RegionMultiphaseStatistics.cpp b/src/coreComponents/physicsSolvers/fluidFlow/RegionMultiphaseStatistics.cpp index 0cc2e82ffaa..2b3e072f674 100644 --- a/src/coreComponents/physicsSolvers/fluidFlow/RegionMultiphaseStatistics.cpp +++ b/src/coreComponents/physicsSolvers/fluidFlow/RegionMultiphaseStatistics.cpp @@ -78,7 +78,11 @@ struct RegionMultiphaseStatistics::Statistics return *this; } - Statistics( Statistics && ) = delete; + Statistics( Statistics && rhs ) + { + std::swap( m_data, rhs.m_data ); + } + ~Statistics() = default; GEOS_HOST_DEVICE bool operator!=( Statistics const & ) const diff --git a/src/coreComponents/physicsSolvers/fluidFlow/RegionMultiphaseStatistics.hpp b/src/coreComponents/physicsSolvers/fluidFlow/RegionMultiphaseStatistics.hpp index 2c7164cbee0..fa99e593fe0 100644 --- a/src/coreComponents/physicsSolvers/fluidFlow/RegionMultiphaseStatistics.hpp +++ b/src/coreComponents/physicsSolvers/fluidFlow/RegionMultiphaseStatistics.hpp @@ -53,6 +53,7 @@ class RegionMultiphaseStatistics : public FieldStatisticsBase< CompositionalMult public: struct Statistics; + struct RegionStatisticsKernel; public: @@ -83,6 +84,8 @@ class RegionMultiphaseStatistics : public FieldStatisticsBase< CompositionalMult /**@}*/ + void initializePostInitialConditionsPreSubGroups() override; + private: using Base = FieldStatisticsBase< CompositionalMultiphaseBase >; @@ -110,8 +113,6 @@ class RegionMultiphaseStatistics : public FieldStatisticsBase< CompositionalMult void registerDataOnMesh( Group & meshBodies ) override; - void initializePostInitialConditionsPreSubGroups() override; - template< typename LAMBDA > void forRegions( Group & meshBodies, LAMBDA && lambda ); @@ -140,8 +141,6 @@ class RegionMultiphaseStatistics : public FieldStatisticsBase< CompositionalMult return (val < LvArray::NumericLimits< real64 >::epsilon) ? 0.0 : 1.0 / val; } - struct RegionStatisticsKernel; - // The names of regions string_array m_regionNames;