diff --git a/src/coreComponents/physicsSolvers/fluidFlow/CompositionalMultiphaseBase.cpp b/src/coreComponents/physicsSolvers/fluidFlow/CompositionalMultiphaseBase.cpp index efbb0f9ff6c..f49c477ad05 100644 --- a/src/coreComponents/physicsSolvers/fluidFlow/CompositionalMultiphaseBase.cpp +++ b/src/coreComponents/physicsSolvers/fluidFlow/CompositionalMultiphaseBase.cpp @@ -1877,6 +1877,7 @@ void CompositionalMultiphaseBase::applySourceFluxBC( real64 const time, localIndex const rankOffset = dofManager.rankOffset(); RAJA::ReduceSum< parallelDeviceReduce, real64 > massProd( 0.0 ); + RAJA::ReduceSum< parallelDeviceReduce, localIndex > elementCount( 0 ); // note that the dofArray will not be used after this step (simpler to use dofNumber instead) FieldSpecificationImpl::computeRhsContribution< FieldSpecificationAdd, @@ -1913,7 +1914,8 @@ void CompositionalMultiphaseBase::applySourceFluxBC( real64 const time, dofNumber, rhsContributionArrayView, localRhs, - massProd] GEOS_HOST_DEVICE ( localIndex const a ) + massProd, + elementCount] GEOS_HOST_DEVICE ( localIndex const a ) { // we need to filter out ghosts here, because targetSet may contain them localIndex const ei = targetSet[a]; @@ -1924,6 +1926,7 @@ void CompositionalMultiphaseBase::applySourceFluxBC( real64 const time, real64 const rhsValue = rhsContributionArrayView[a] / sizeScalingFactor; // scale the contribution by the sizeScalingFactor here! massProd += rhsValue; + elementCount += 1; if( useTotalMassEquation > 0 ) { // for all "fluid components", we add the value to the total mass balance equation @@ -1948,7 +1951,7 @@ void CompositionalMultiphaseBase::applySourceFluxBC( real64 const time, // set the new sub-region statistics for this timestep array1d< real64 > massProdArr{ m_numComponents }; massProdArr[fluidComponentId] = massProd.get(); - wrapper.gatherTimeStepStats( time, dt, massProdArr.toViewConst(), targetSet.size() ); + wrapper.gatherTimeStepStats( time, dt, massProdArr.toViewConst(), elementCount.get() ); } ); } ); } ); diff --git a/src/coreComponents/physicsSolvers/fluidFlow/SinglePhaseBase.cpp b/src/coreComponents/physicsSolvers/fluidFlow/SinglePhaseBase.cpp index 562cfae6df4..5b04db3c016 100644 --- a/src/coreComponents/physicsSolvers/fluidFlow/SinglePhaseBase.cpp +++ b/src/coreComponents/physicsSolvers/fluidFlow/SinglePhaseBase.cpp @@ -1055,6 +1055,7 @@ void SinglePhaseBase::applySourceFluxBC( real64 const time_n, localIndex const rankOffset = dofManager.rankOffset(); RAJA::ReduceSum< parallelDeviceReduce, real64 > massProd( 0.0 ); + RAJA::ReduceSum< parallelDeviceReduce, localIndex > elementCount( 0 ); // note that the dofArray will not be used after this step (simpler to use dofNumber instead) FieldSpecificationImpl::computeRhsContribution< FieldSpecificationAdd, @@ -1096,7 +1097,8 @@ void SinglePhaseBase::applySourceFluxBC( real64 const time_n, rhsContributionArrayView, localRhs, localMatrix, - massProd] GEOS_HOST_DEVICE ( localIndex const a ) + massProd, + elementCount] GEOS_HOST_DEVICE ( localIndex const a ) { // we need to filter out ghosts here, because targetSet may contain them localIndex const ei = targetSet[a]; @@ -1111,6 +1113,7 @@ void SinglePhaseBase::applySourceFluxBC( real64 const time_n, real64 const rhsValue = rhsContributionArrayView[a] / sizeScalingFactor; // scale the contribution by the sizeScalingFactor here! localRhs[massRowIndex] += rhsValue; massProd += rhsValue; + elementCount += 1; //add the value to the energy balance equation if the flux is positive (i.e., it's a producer) if( rhsContributionArrayView[a] > 0.0 ) { @@ -1138,7 +1141,8 @@ void SinglePhaseBase::applySourceFluxBC( real64 const time_n, dofNumber, rhsContributionArrayView, localRhs, - massProd] GEOS_HOST_DEVICE ( localIndex const a ) + massProd, + elementCount] GEOS_HOST_DEVICE ( localIndex const a ) { // we need to filter out ghosts here, because targetSet may contain them localIndex const ei = targetSet[a]; @@ -1152,6 +1156,7 @@ void SinglePhaseBase::applySourceFluxBC( real64 const time_n, real64 const rhsValue = rhsContributionArrayView[a] / sizeScalingFactor; localRhs[rowIndex] += rhsValue; massProd += rhsValue; + elementCount += 1; } ); } @@ -1161,7 +1166,7 @@ void SinglePhaseBase::applySourceFluxBC( real64 const time_n, // set the new sub-region statistics for this timestep array1d< real64 > massProdArr{ 1 }; massProdArr[0] = massProd.get(); - wrapper.gatherTimeStepStats( time_n, dt, massProdArr.toViewConst(), targetSet.size() ); + wrapper.gatherTimeStepStats( time_n, dt, massProdArr.toViewConst(), elementCount.get() ); } ); } ); } ); diff --git a/src/coreComponents/physicsSolvers/fluidFlow/SourceFluxStatistics.cpp b/src/coreComponents/physicsSolvers/fluidFlow/SourceFluxStatistics.cpp index 26f734efd64..5ae6cf4152e 100644 --- a/src/coreComponents/physicsSolvers/fluidFlow/SourceFluxStatistics.cpp +++ b/src/coreComponents/physicsSolvers/fluidFlow/SourceFluxStatistics.cpp @@ -139,6 +139,14 @@ void SourceFluxStatsAggregator::registerDataOnMesh( Group & meshBodies ) } ); } } ); + + if( m_writeCSV > 0 && MpiWrapper::commRank() == 0 ) + { + std::ofstream outputFile( m_csvFilename ); + TableCSVFormatter const tableStatFormatter( m_csvLayout ); + outputFile << tableStatFormatter.headerToString(); + outputFile.close(); + } } void SourceFluxStatsAggregator::gatherStatsForLog( bool logLevelActive, @@ -215,9 +223,9 @@ void SourceFluxStatsAggregator::outputStatsToCSV( TableData & csvData ) { if( m_writeCSV > 0 && MpiWrapper::commRank() == 0 ) { - std::ofstream outputFile( m_csvFilename ); + std::ofstream outputFile( m_csvFilename, std::ios::app ); TableCSVFormatter const tableStatFormatter( m_csvLayout ); - outputFile << tableStatFormatter.toString( csvData ); + outputFile << tableStatFormatter.dataToString( csvData ); outputFile.close(); csvData.clear(); } @@ -241,30 +249,31 @@ bool SourceFluxStatsAggregator::execute( real64 const GEOS_UNUSED_PARAM( time_n { TableData logData; TableData csvData; - meshLevelStats.stats() = StatData(); + meshLevelStats.stats().reset(); forAllFluxStatsWrappers( meshLevel, [&] ( MeshLevel &, WrappedStats & fluxStats ) { - fluxStats.stats() = StatData(); + fluxStats.stats().reset(); forAllRegionStatsWrappers( meshLevel, fluxStats.getFluxName(), [&] ( ElementRegionBase & region, WrappedStats & regionStats ) { - regionStats.stats() = StatData(); + regionStats.stats().reset(); forAllSubRegionStatsWrappers( region, regionStats.getFluxName(), [&] ( ElementSubRegionBase &, WrappedStats & subRegionStats ) { + subRegionStats.stats().reset(); subRegionStats.finalizePeriod(); - regionStats.stats().combine( subRegionStats.stats() ); + regionStats.combine( subRegionStats ); } ); - fluxStats.stats().combine( regionStats.stats() ); + fluxStats.combine( regionStats ); gatherStatsForLog( regionsStatsOn, fluxStats.getFluxName(), region.getName(), logData, regionStats ); gatherStatsForCSV( fluxStats.getFluxName(), region.getName(), csvData, regionStats ); } ); - meshLevelStats.stats().combine( fluxStats.stats() ); + meshLevelStats.combine( fluxStats ); gatherStatsForLog( fluxesStatsOn, fluxStats.getFluxName(), allRegionsStr, logData, fluxStats ); @@ -370,7 +379,9 @@ void SourceFluxStatsAggregator::WrappedStats::finalizePeriod() // produce the period stats of this rank m_stats.m_elementCount = m_periodStats.m_elementCount; - m_statsPeriodStart = m_periodStats.m_periodStart; + // MPI ranks without cells in the flux region may bypass gatherTimeStepStats() entirely, + // leaving their m_periodStart behind. Take the max so all ranks agree on the most advanced time. + m_statsPeriodStart = MpiWrapper::max( m_periodStats.m_periodStart ); m_statsPeriodDT = m_periodStats.m_timeStepDeltaTime + m_periodStats.m_periodPendingDeltaTime; real64 const timeDivisor = m_statsPeriodDT > 0.0 ? 1.0 / m_statsPeriodDT : 0.0; @@ -387,6 +398,13 @@ void SourceFluxStatsAggregator::WrappedStats::finalizePeriod() // start a new timestep m_periodStats.reset(); } +void SourceFluxStatsAggregator::WrappedStats::combine( WrappedStats const & other ) +{ + stats().combine( other.stats() ); + // Some fluxes may lag behind (their ranks may have skipped gatherTimeStepStats()), + // so take the most advanced period start time across all combined stats. + m_statsPeriodStart = LvArray::math::max( m_statsPeriodStart, other.m_statsPeriodStart ); +} void SourceFluxStatsAggregator::WrappedStats::PeriodStats::allocate( integer phaseCount ) { if( m_timeStepMass.size() < phaseCount ) diff --git a/src/coreComponents/physicsSolvers/fluidFlow/SourceFluxStatistics.hpp b/src/coreComponents/physicsSolvers/fluidFlow/SourceFluxStatistics.hpp index b777cd3ba07..5b015e22232 100644 --- a/src/coreComponents/physicsSolvers/fluidFlow/SourceFluxStatistics.hpp +++ b/src/coreComponents/physicsSolvers/fluidFlow/SourceFluxStatistics.hpp @@ -112,6 +112,14 @@ class SourceFluxStatsAggregator final : public FieldStatisticsBase< FlowSolverBa */ void finalizePeriod(); + /** + * @brief Aggregate the statistics of the instance with those of another one. + * @details Combines the statistics from the other object into this one and also advances the period + * start if the other object has a later time recorded. + * @param other the other WrappedStats object. + */ + void combine( WrappedStats const & other ); + /** * @return the reference to the wrapped stats data collected over the last period (one timestep or more), computed by finalizePeriod() */ @@ -151,7 +159,7 @@ class SourceFluxStatsAggregator final : public FieldStatisticsBase< FlowSolverBa /// stats data collected over the last period (one timestep or more), computed by finalizePeriod() StatData m_stats; /// the start time of the wrapped stats period (in s) - real64 m_statsPeriodStart; + real64 m_statsPeriodStart{-LvArray::NumericLimits< real64 >::max}; /// the duration of the wrapped stats period (in s) real64 m_statsPeriodDT; @@ -169,7 +177,7 @@ class SourceFluxStatsAggregator final : public FieldStatisticsBase< FlowSolverBa /// time that the current timestep is simulating. real64 m_timeStepDeltaTime = 0.0; /// start time of the current period. - real64 m_periodStart = 0.0; + real64 m_periodStart = -LvArray::NumericLimits< real64 >::max; /// delta time from all previous time-step of the current period. real64 m_periodPendingDeltaTime = 0.0; /// number of cell elements targeted by this instance diff --git a/src/coreComponents/physicsSolvers/fluidFlow/unitTests/CMakeLists.txt b/src/coreComponents/physicsSolvers/fluidFlow/unitTests/CMakeLists.txt index 77d628e7ae2..f59e71ec23b 100644 --- a/src/coreComponents/physicsSolvers/fluidFlow/unitTests/CMakeLists.txt +++ b/src/coreComponents/physicsSolvers/fluidFlow/unitTests/CMakeLists.txt @@ -17,4 +17,11 @@ foreach( test_source ${fluidFlow_gtest_tests} ) FOLDER UnitTests ) geos_add_test( NAME ${test_name} COMMAND ${test_name} ) + + if( ENABLE_MPI AND NOT ENABLE_CUDA AND NOT ENABLE_HIP ) + set( nranks 2 ) + geos_add_test( NAME ${test_name}_mpi + COMMAND ${test_name} -x ${nranks} + NUM_MPI_TASKS ${nranks} ) + endif() endforeach() diff --git a/src/coreComponents/physicsSolvers/fluidFlow/unitTests/testFlowStatistics.cpp b/src/coreComponents/physicsSolvers/fluidFlow/unitTests/testFlowStatistics.cpp index 356661e434e..61d85956613 100644 --- a/src/coreComponents/physicsSolvers/fluidFlow/unitTests/testFlowStatistics.cpp +++ b/src/coreComponents/physicsSolvers/fluidFlow/unitTests/testFlowStatistics.cpp @@ -19,6 +19,7 @@ #include "mainInterface/GeosxState.hpp" #include "physicsSolvers/fluidFlow/SourceFluxStatistics.hpp" #include "physicsSolvers/fluidFlow/SinglePhaseStatistics.hpp" +#include "common/MpiWrapper.hpp" #include @@ -154,15 +155,19 @@ class FlowStatisticsTest : public ::testing::Test void writeTableFiles( stdMap< string, string > const & files ) { - for( auto const & [fileName, content] : files ) + if( MpiWrapper::commRank() == 0 ) { - std::ofstream os( fileName ); - ASSERT_TRUE( os.is_open() ); - os << content; - os.close(); + for( auto const & [fileName, content] : files ) + { + std::ofstream os( fileName ); + ASSERT_TRUE( os.is_open() ); + os << content; + os.close(); - m_tableFileNames.push_back( fileName ); + m_tableFileNames.push_back( fileName ); + } } + MpiWrapper::barrier(); } void TearDown() override