Implement arbitrary region statistics by dkachuma · Pull Request #3037 · GEOS-DEV/GEOS

dkachuma · 2024-03-11T21:45:16Z

The purpose of this PR is to allow the user to report region statistics for arbitrary cell regions (cell collections) that are not linked to the simulation (solvers or constitutive models).

The idea is that the user will create a task similar to CompositionalMultiphaseStatistics. With that they can specify the properties to report and then target this with an event.

    <Tasks>
        <RegionMultiphaseStatistics
            name="REGION.STATISTICS"
            flowSolverName="FLOW.SOLVER"
            regionNames="{ PLATFORM_A, PLATFORM_B1, PLATFORM_B2 }"
            regionIdentifiers="{ 13, 28, 3 }"
            propertyNames="{ pressure, poreVolume, volumeFraction }"
            writeCSV="1"
            logLevel="1" />
    </Tasks>
    <Events>
        <PeriodicEvent
            name="REGION.STATISTICS"
            target="/Tasks/REGION.STATISTICS"
            timeFrequency="7200" />
    </Events>

In this example, there are 3 non-overlapping regions named PLATFORM_A, PLATFORM_B1 and PLATFORM_B2 not linked to the simulation and the user wants statistics for pressure, pore volume and volume fraction in these regions. The 3 regions are identified by an attribute and these are assigned values 13, 28 and 3 respectively. The creation of the task will create a fieldName REGION.STATISTICS_region for each grid cell. The user can use this field to load a property and mark the cells belonging to each region. For example

    <Mesh>
        <VTKMesh
            name="MESH"
            file="mesh.vtu"
            regionAttribute="attribute"
            fieldsToImport="{ report }"
            fieldNamesInGEOSX="{ REGION.STATISTICS_region }"
            logLevel="5" />
    </Mesh>

jafranc · 2024-03-13T16:12:12Z

Great idea. Can we make it work with '' define boxes too ?

MelReyCG

Here is a batch of (non-blocking) comments. we could meet to discuss about them.

MelReyCG · 2024-03-26T16:18:31Z

+  GEOS_HOST_DEVICE
+  static constexpr real64 safeInverse( real64 const val )
+  {
+    return (val < LvArray::NumericLimits< real64 >::epsilon) ? 0.0 : 1.0 / val;
+  }


Could this be exposed in a header that is commonly used for calculations? I think safe-dividing could be useful to everyone. :)

MelReyCG · 2024-03-26T16:20:25Z

+    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] );
+
+    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] );
+
+    regData[Statistics::AVERAGE_PRESSURE] = MpiWrapper::sum( regData[Statistics::AVERAGE_PRESSURE] );
+    regData[Statistics::AVERAGE_TEMPERATURE] = MpiWrapper::sum( regData[Statistics::AVERAGE_TEMPERATURE] );
+
+    for( integer propIndex = Statistics::PHASE_PORE_VOLUME; propIndex < Statistics::END; ++propIndex )
+    {
+      regData[propIndex] = MpiWrapper::sum( regData[propIndex] );
+    }


I don't know the impact, but I'm wondering if this kind of repeated reductions should & could be grouped to be optimized.

MelReyCG · 2024-03-27T10:01:13Z

+
+  registerWrapper( viewKeyStruct::propertyNamesString(), &m_propertyNames ).
+    setInputFlag( InputFlags::OPTIONAL ).
+    setDescription( "The list of names of properties to report" );


Suggested change

setDescription( "The list of names of properties to report" );

setDescription( "The list of names of properties to report. If leaved empty, all available properties are selected." );

MelReyCG · 2024-03-27T10:49:01Z

+    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::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 );


In order to automatically iterate over enum labels :

Suggested change

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::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 );

constexpr auto labelCount = EnumStrings< PropertyNameType >::count;

for( integer i = 0; i < labelCount; ++i )

{

m_propertyNameTypes.emplace_back( i );

}

... we could add a getEnumStringsCount( ENUM const ) in the ENUM_STRINGS() macro :

inline constexpr size_t getEnumStringsCount( ENUM const ) \ { \ constexpr char const * arr[] { __VA_ARGS__ }; \ constexpr auto s = std::end( arr ) - std::begin( arr ); \ return s; \ } \

... and a constexpr count in EnumStrings :

inline static constexpr size_t count = getEnumStringsCount( enum_type{} );

MelReyCG · 2024-03-27T10:54:40Z

+  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();


Me or @arng40 will certainly make a future PR to add table output in the log (using the component from this PR) and to the csv export with centralized components.

MelReyCG · 2024-03-27T11:16:59Z

+  string const massUnit = useMass ? "kg" : "mol";
+  string const densityUnit = useMass ? "kg/m^3" : "mol/m^3";


Suggested change

string const massUnit = useMass ? "kg" : "mol";

string const densityUnit = useMass ? "kg/m^3" : "mol/m^3";

string const massUnit = getSymbol( useMass ? units::Unit::Mass : units::Unit::Mole );

string const densityUnit = getSymbol( useMass ? units::Unit::MassRate : units::Unit::MoleRate );

MelReyCG · 2024-03-27T15:50:32Z

+struct RegionMultiphaseStatistics::Statistics
+{
+  static constexpr integer MAX_NUM_PHASE = 3;
+  static constexpr integer MAX_NUM_COMPS = 5;
+  /// Indices for properties
+  static constexpr integer CELL_COUNT = 0;
+  static constexpr integer STATIC_PORE_VOLUME = 1;
+  static constexpr integer PORE_VOLUME = 2;


Here, what do you think of grouping the indices in a C Enum?

Suggested change

struct RegionMultiphaseStatistics::Statistics

{

static constexpr integer MAX_NUM_PHASE = 3;

static constexpr integer MAX_NUM_COMPS = 5;

/// Indices for properties

static constexpr integer CELL_COUNT = 0;

static constexpr integer STATIC_PORE_VOLUME = 1;

static constexpr integer PORE_VOLUME = 2;

struct RegionMultiphaseStatistics::Statistics

{

static constexpr integer MAX_NUM_PHASE = 3;

static constexpr integer MAX_NUM_COMPS = 5;

enum PropIds : integer {

CELL_COUNT,

STATIC_PORE_VOLUME,

PORE_VOLUME,

// [...]

END

};

MelReyCG · 2024-03-27T16:13:56Z

+  GEOS_HOST_DEVICE Statistics( Statistics const & rhs )
+  {
+    for( integer i = 0; i < END; i++ )
+    {
+      m_data[i] = rhs.m_data[i];
+    }
+  }
+
+  GEOS_HOST_DEVICE Statistics const & operator=( Statistics const & rhs )
+  {
+    for( integer i = 0; i < END; i++ )
+    {
+      m_data[i] = rhs.m_data[i];
+    }
+    return *this;
+  }
+
+  Statistics( Statistics && rhs )
+  {
+    std::swap( m_data, rhs.m_data );
+  }


Aren't they the default?

MelReyCG · 2024-03-27T16:21:58Z

+    RAJA::ReduceSum< ReducePolicy< POLICY >, Statistics > sumStats( RAJA::operators::plus< Statistics >::identity() );
+    RAJA::ReduceMin< ReducePolicy< POLICY >, Statistics > minStats( RAJA::operators::minimum< Statistics >::identity() );


Is it possible to declare a RAJA::ReduceSum< ReducePolicy< POLICY >, real64[] >?

MelReyCG · 2024-03-28T10:47:41Z

+  } );
+
+  // Create the csv file if requires
+  if( 0 < getLogLevel() && 0 < m_writeCSV && MpiWrapper::commRank() == 0 )


Why this condition is different to the one below? (in RegionMultiphaseStatistics::computeRegionStatistics()).

dkachuma · 2026-05-13T22:58:27Z

Superceded by #4005

dkachuma added 3 commits March 11, 2024 12:21

First pass

4d7e205

Move kernel to implementation file

b91f59c

Restore kernel file

62149b5

dkachuma self-assigned this Mar 12, 2024

dkachuma added the spe11 spe11-required-feature label Mar 12, 2024

Implement averaging

1f8cbe9

dkachuma and others added 4 commits March 14, 2024 16:53

Use element-wise array reduction

a69de32

Use cell count to average density and viscosity

be2bd30

Implement move constructor

4e553bf

Merge branch 'develop' into feature/dkachuma/region-statistics

09e5307

dkachuma requested a review from MelReyCG March 25, 2024 20:05

Merge branch 'develop' into feature/dkachuma/region-statistics

74752aa

MelReyCG reviewed Mar 28, 2024

View reviewed changes

dkachuma closed this May 13, 2026

	setDescription( "The list of names of properties to report" );
	setDescription( "The list of names of properties to report. If leaved empty, all available properties are selected." );

		string const massUnit = useMass ? "kg" : "mol";
		string const densityUnit = useMass ? "kg/m^3" : "mol/m^3";

		RAJA::ReduceSum< ReducePolicy< POLICY >, Statistics > sumStats( RAJA::operators::plus< Statistics >::identity() );
		RAJA::ReduceMin< ReducePolicy< POLICY >, Statistics > minStats( RAJA::operators::minimum< Statistics >::identity() );

Conversation

dkachuma commented Mar 11, 2024 • edited Loading Uh oh! There was an error while loading. Please reload this page.

Uh oh!

Uh oh!

jafranc commented Mar 13, 2024

Uh oh!

MelReyCG left a comment • edited Loading Uh oh! There was an error while loading. Please reload this page.

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

MelReyCG Mar 27, 2024 • edited Loading Uh oh! There was an error while loading. Please reload this page.

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

dkachuma commented May 13, 2026

Uh oh!

Reviewers

Assignees

Labels

Projects

Milestone

Development

Uh oh!

3 participants

dkachuma commented Mar 11, 2024 •

edited

Loading

MelReyCG left a comment •

edited

Loading

MelReyCG Mar 27, 2024 •

edited

Loading