From be0eb0e2e92d933f0f976a2c341d0ab4d9289ebf Mon Sep 17 00:00:00 2001 From: Pavel Tomin Date: Fri, 6 Dec 2024 21:03:14 -0600 Subject: [PATCH 1/2] fix: make new gravity treatment from #3337 an option (#3467) --- .integrated_tests.yaml | 2 +- BASELINE_NOTES.md | 4 + .../fluidFlow/CompositionalMultiphaseBase.cpp | 8 + .../fluidFlow/CompositionalMultiphaseBase.hpp | 4 + .../fluidFlow/CompositionalMultiphaseFVM.cpp | 1 + .../kernels/compositional/C1PPUPhaseFlux.hpp | 3 +- .../kernels/compositional/CFLKernel.cpp | 102 +++++++------ .../kernels/compositional/CFLKernel.hpp | 18 ++- .../DissipationFluxComputeKernel.hpp | 1 + .../compositional/FluxComputeKernel.hpp | 9 +- .../compositional/FluxComputeKernelBase.hpp | 7 +- .../kernels/compositional/IHUPhaseFlux.hpp | 138 +++++++++++------- .../kernels/compositional/PPUPhaseFlux.hpp | 3 +- .../kernels/compositional/PotGrad.hpp | 111 ++++++++------ .../StabilizedFluxComputeKernel.hpp | 1 + .../ThermalFluxComputeKernel.hpp | 3 +- 16 files changed, 259 insertions(+), 156 deletions(-) diff --git a/.integrated_tests.yaml b/.integrated_tests.yaml index 97f13dd9062..0ce10cbe33b 100644 --- a/.integrated_tests.yaml +++ b/.integrated_tests.yaml @@ -1,6 +1,6 @@ baselines: bucket: geosx - baseline: integratedTests/baseline_integratedTests-pr3361-9139-2fc4131 + baseline: integratedTests/baseline_integratedTests-pr3467-9212-976cc3b allow_fail: all: '' streak: '' diff --git a/BASELINE_NOTES.md b/BASELINE_NOTES.md index e5e337912cb..8f344877826 100644 --- a/BASELINE_NOTES.md +++ b/BASELINE_NOTES.md @@ -6,6 +6,10 @@ This file is designed to track changes to the integrated test baselines. Any developer who updates the baseline ID in the .integrated_tests.yaml file is expected to create an entry in this file with the pull request number, date, and their justification for rebaselining. These notes should be in reverse-chronological order, and use the following time format: (YYYY-MM-DD). +PR #3361 (2024-12-03) +===================== +Revert default gravity treatment to old version. Make the way introduced in #3337 optional. + PR #3361 (2024-12-03) ===================== Baseline diffs after reimplementation of wave equation acoustic gradient for velocity and density parameters: new field "partialGradient2" and "pressureForward" field replacing "pressureDoubleDerivative". diff --git a/src/coreComponents/physicsSolvers/fluidFlow/CompositionalMultiphaseBase.cpp b/src/coreComponents/physicsSolvers/fluidFlow/CompositionalMultiphaseBase.cpp index 9c7d0108c7b..9eb10dbab90 100644 --- a/src/coreComponents/physicsSolvers/fluidFlow/CompositionalMultiphaseBase.cpp +++ b/src/coreComponents/physicsSolvers/fluidFlow/CompositionalMultiphaseBase.cpp @@ -78,6 +78,7 @@ CompositionalMultiphaseBase::CompositionalMultiphaseBase( const string & name, m_allowCompDensChopping( 1 ), m_useTotalMassEquation( 1 ), m_useSimpleAccumulation( 1 ), + m_useNewGravity( 0 ), m_minCompDens( isothermalCompositionalMultiphaseBaseKernels::minDensForDivision ) { //START_SPHINX_INCLUDE_00 @@ -164,6 +165,12 @@ CompositionalMultiphaseBase::CompositionalMultiphaseBase( const string & name, setApplyDefaultValue( 1 ). setDescription( "Flag indicating whether simple accumulation form is used" ); + this->registerWrapper( viewKeyStruct::useNewGravityString(), &m_useNewGravity ). + setSizedFromParent( 0 ). + setInputFlag( InputFlags::OPTIONAL ). + setApplyDefaultValue( 0 ). + setDescription( "Flag indicating whether new gravity treatment is used" ); + this->registerWrapper( viewKeyStruct::minCompDensString(), &m_minCompDens ). setSizedFromParent( 0 ). setInputFlag( InputFlags::OPTIONAL ). @@ -2245,6 +2252,7 @@ void CompositionalMultiphaseBase::computeCFLNumbers( geos::DomainPartition & dom isothermalCompositionalMultiphaseBaseKernels::KernelLaunchSelector1 < isothermalCompositionalMultiphaseFVMKernels::CFLFluxKernel >( numComps, numPhases, + m_useNewGravity, dt, stencilWrapper, compFlowAccessors.get( fields::flow::pressure{} ), diff --git a/src/coreComponents/physicsSolvers/fluidFlow/CompositionalMultiphaseBase.hpp b/src/coreComponents/physicsSolvers/fluidFlow/CompositionalMultiphaseBase.hpp index 278d8d6b453..1e657613848 100644 --- a/src/coreComponents/physicsSolvers/fluidFlow/CompositionalMultiphaseBase.hpp +++ b/src/coreComponents/physicsSolvers/fluidFlow/CompositionalMultiphaseBase.hpp @@ -268,6 +268,7 @@ class CompositionalMultiphaseBase : public FlowSolverBase static constexpr char const * allowLocalCompDensChoppingString() { return "allowLocalCompDensityChopping"; } static constexpr char const * useTotalMassEquationString() { return "useTotalMassEquation"; } static constexpr char const * useSimpleAccumulationString() { return "useSimpleAccumulation"; } + static constexpr char const * useNewGravityString() { return "useNewGravity"; } static constexpr char const * minCompDensString() { return "minCompDens"; } static constexpr char const * maxSequentialCompDensChangeString() { return "maxSequentialCompDensChange"; } static constexpr char const * minScalingFactorString() { return "minScalingFactor"; } @@ -486,6 +487,9 @@ class CompositionalMultiphaseBase : public FlowSolverBase /// flag indicating whether simple accumulation form is used integer m_useSimpleAccumulation; + /// flag indicating whether new gravity treatment is used + integer m_useNewGravity; + /// minimum allowed global component density real64 m_minCompDens; diff --git a/src/coreComponents/physicsSolvers/fluidFlow/CompositionalMultiphaseFVM.cpp b/src/coreComponents/physicsSolvers/fluidFlow/CompositionalMultiphaseFVM.cpp index 2e99cfe6f48..8eec4acbef0 100644 --- a/src/coreComponents/physicsSolvers/fluidFlow/CompositionalMultiphaseFVM.cpp +++ b/src/coreComponents/physicsSolvers/fluidFlow/CompositionalMultiphaseFVM.cpp @@ -231,6 +231,7 @@ void CompositionalMultiphaseFVM::assembleFluxTerms( real64 const dt, elemDofKey, m_hasCapPressure, m_useTotalMassEquation, + m_useNewGravity, fluxApprox.upwindingParams(), getName(), mesh.getElemManager(), diff --git a/src/coreComponents/physicsSolvers/fluidFlow/kernels/compositional/C1PPUPhaseFlux.hpp b/src/coreComponents/physicsSolvers/fluidFlow/kernels/compositional/C1PPUPhaseFlux.hpp index a254aa9f97d..07e5d54df99 100644 --- a/src/coreComponents/physicsSolvers/fluidFlow/kernels/compositional/C1PPUPhaseFlux.hpp +++ b/src/coreComponents/physicsSolvers/fluidFlow/kernels/compositional/C1PPUPhaseFlux.hpp @@ -79,6 +79,7 @@ struct C1PPUPhaseFlux compute( integer const numPhase, integer const ip, integer const hasCapPressure, + integer const useNewGravity, localIndex const ( &seri )[numFluxSupportPoints], localIndex const ( &sesri )[numFluxSupportPoints], localIndex const ( &sei )[numFluxSupportPoints], @@ -110,7 +111,7 @@ struct C1PPUPhaseFlux real64 dPresGrad_dC[numFluxSupportPoints][numComp]{}; real64 dGravHead_dP[numFluxSupportPoints]{}; real64 dGravHead_dC[numFluxSupportPoints][numComp]{}; - PotGrad::compute< numComp, numFluxSupportPoints >( numPhase, ip, hasCapPressure, seri, sesri, sei, trans, dTrans_dPres, pres, + PotGrad::compute< numComp, numFluxSupportPoints >( numPhase, ip, hasCapPressure, useNewGravity, seri, sesri, sei, trans, dTrans_dPres, pres, gravCoef, phaseVolFrac, dPhaseVolFrac, dCompFrac_dCompDens, phaseMassDens, dPhaseMassDens, phaseCapPressure, dPhaseCapPressure_dPhaseVolFrac, potGrad, dPresGrad_dP, dPresGrad_dC, dGravHead_dP, dGravHead_dC ); diff --git a/src/coreComponents/physicsSolvers/fluidFlow/kernels/compositional/CFLKernel.cpp b/src/coreComponents/physicsSolvers/fluidFlow/kernels/compositional/CFLKernel.cpp index b416573af0b..e2f32348bcc 100644 --- a/src/coreComponents/physicsSolvers/fluidFlow/kernels/compositional/CFLKernel.cpp +++ b/src/coreComponents/physicsSolvers/fluidFlow/kernels/compositional/CFLKernel.cpp @@ -22,6 +22,7 @@ #include "finiteVolume/SurfaceElementStencil.hpp" #include "finiteVolume/EmbeddedSurfaceToCellStencil.hpp" #include "finiteVolume/FaceElementToCellStencil.hpp" +#include "CFLKernel.hpp" namespace geos { @@ -32,12 +33,13 @@ namespace isothermalCompositionalMultiphaseFVMKernels /******************************** CFLFluxKernel ********************************/ -template< integer NC, localIndex NUM_ELEMS, localIndex maxStencilSize > +template< integer NC > GEOS_HOST_DEVICE inline void CFLFluxKernel:: compute( integer const numPhases, + integer const useNewGravity, localIndex const stencilSize, real64 const dt, arraySlice1d< localIndex const > const seri, @@ -67,27 +69,7 @@ CFLFluxKernel:: real64 gravHead{}; // calculate quantities on primary connected cells - integer denom = 0; - for( localIndex i = 0; i < NUM_ELEMS; ++i ) - { - localIndex const er = seri[i]; - localIndex const esr = sesri[i]; - localIndex const ei = sei[i]; - - bool const phaseExists = (phaseVolFrac[er][esr][ei][ip] > 0); - if( !phaseExists ) - { - continue; - } - - // average density across the face - densMean += phaseMassDens[er][esr][ei][0][ip]; - denom++; - } - if( denom > 1 ) - { - densMean /= denom; - } + calculateMeanDensity( useNewGravity, ip, stencilSize, seri, sesri, sei, phaseVolFrac, phaseMassDens, densMean ); //***** calculation of phase volumetric flux ***** @@ -138,10 +120,43 @@ CFLFluxKernel:: } } -template< integer NC, typename STENCILWRAPPER_TYPE > +GEOS_HOST_DEVICE +inline void -CFLFluxKernel:: +CFLFluxKernel::calculateMeanDensity( integer const useNewGravity, integer const ip, localIndex const stencilSize, + arraySlice1d< localIndex const > const seri, + arraySlice1d< localIndex const > const sesri, + arraySlice1d< localIndex const > const sei, + ElementViewConst< arrayView2d< real64 const, compflow::USD_PHASE > > const & phaseVolFrac, + ElementViewConst< arrayView3d< real64 const, multifluid::USD_PHASE > > const & phaseMassDens, + real64 & densMean ) +{ + integer denom = 0; + for( localIndex i = 0; i < stencilSize; ++i ) + { + localIndex const er = seri[i]; + localIndex const esr = sesri[i]; + localIndex const ei = sei[i]; + + bool const phaseExists = (phaseVolFrac[er][esr][ei][ip] > 0); + if( useNewGravity && !phaseExists ) + { + continue; + } + + // average density across the face + densMean += phaseMassDens[er][esr][ei][0][ip]; + denom++; + } + if( denom > 1 ) + { + densMean /= denom; + } +} +template< integer NC, typename STENCILWRAPPER_TYPE > +void CFLFluxKernel:: launch( integer const numPhases, + integer const useNewGravity, real64 const dt, STENCILWRAPPER_TYPE const & stencilWrapper, ElementViewConst< arrayView1d< real64 const > > const & pres, @@ -161,9 +176,6 @@ CFLFluxKernel:: typename STENCILWRAPPER_TYPE::IndexContainerViewConstType const & sesri = stencilWrapper.getElementSubRegionIndices(); typename STENCILWRAPPER_TYPE::IndexContainerViewConstType const & sei = stencilWrapper.getElementIndices(); - localIndex constexpr numElems = STENCILWRAPPER_TYPE::maxNumPointsInFlux; - localIndex constexpr maxStencilSize = STENCILWRAPPER_TYPE::maxStencilSize; - forAll< parallelDevicePolicy<> >( stencilWrapper.size(), [=] GEOS_HOST_DEVICE ( localIndex const iconn ) { // compute transmissibility @@ -176,23 +188,24 @@ CFLFluxKernel:: transmissibility, dTrans_dPres ); - CFLFluxKernel::compute< NC, numElems, maxStencilSize >( numPhases, - sei[iconn].size(), - dt, - seri[iconn], - sesri[iconn], - sei[iconn], - transmissibility[0], - pres, - gravCoef, - phaseVolFrac, - phaseRelPerm, - phaseVisc, - phaseDens, - phaseMassDens, - phaseCompFrac, - phaseOutflux, - compOutflux ); + CFLFluxKernel::compute< NC >( numPhases, + useNewGravity, + sei[iconn].size(), + dt, + seri[iconn], + sesri[iconn], + sei[iconn], + transmissibility[0], + pres, + gravCoef, + phaseVolFrac, + phaseRelPerm, + phaseVisc, + phaseDens, + phaseMassDens, + phaseCompFrac, + phaseOutflux, + compOutflux ); } ); } @@ -200,6 +213,7 @@ CFLFluxKernel:: template \ void CFLFluxKernel:: \ launch< NC, STENCILWRAPPER_TYPE >( integer const numPhases, \ + integer const useNewGravity, \ real64 const dt, \ STENCILWRAPPER_TYPE const & stencil, \ ElementViewConst< arrayView1d< real64 const > > const & pres, \ diff --git a/src/coreComponents/physicsSolvers/fluidFlow/kernels/compositional/CFLKernel.hpp b/src/coreComponents/physicsSolvers/fluidFlow/kernels/compositional/CFLKernel.hpp index 8449d57e43e..a2242e202fb 100644 --- a/src/coreComponents/physicsSolvers/fluidFlow/kernels/compositional/CFLKernel.hpp +++ b/src/coreComponents/physicsSolvers/fluidFlow/kernels/compositional/CFLKernel.hpp @@ -23,6 +23,7 @@ #include "common/DataLayouts.hpp" #include "common/DataTypes.hpp" #include "common/GEOS_RAJA_Interface.hpp" +#include "constitutive/fluid/multifluid/Layouts.hpp" #include "constitutive/fluid/multifluid/MultiFluidBase.hpp" #include "constitutive/fluid/multifluid/MultiFluidFields.hpp" #include "constitutive/permeability/PermeabilityBase.hpp" @@ -83,11 +84,10 @@ struct CFLFluxKernel using RelPermAccessors = StencilMaterialAccessors< constitutive::RelativePermeabilityBase, fields::relperm::phaseRelPerm >; - template< integer NC, localIndex NUM_ELEMS, localIndex maxStencilSize > - GEOS_HOST_DEVICE - inline - static void + template< integer NC > + GEOS_HOST_DEVICE inline static void compute( integer const numPhases, + integer const useNewGravity, localIndex const stencilSize, real64 const dt, arraySlice1d< localIndex const > const seri, @@ -105,9 +105,19 @@ struct CFLFluxKernel ElementView< arrayView2d< real64, compflow::USD_PHASE > > const & phaseOutflux, ElementView< arrayView2d< real64, compflow::USD_COMP > > const & compOutflux ); + GEOS_HOST_DEVICE inline static void + calculateMeanDensity( integer const useNewGravity, integer const ip, localIndex const stencilSize, + arraySlice1d< localIndex const > const seri, + arraySlice1d< localIndex const > const sesri, + arraySlice1d< localIndex const > const sei, + ElementViewConst< arrayView2d< real64 const, compflow::USD_PHASE > > const & phaseVolFrac, + ElementViewConst< arrayView3d< real64 const, constitutive::multifluid::USD_PHASE > > const & phaseMassDens, + real64 & densMean ); + template< integer NC, typename STENCILWRAPPER_TYPE > static void launch( integer const numPhases, + integer const useNewGravity, real64 const dt, STENCILWRAPPER_TYPE const & stencil, ElementViewConst< arrayView1d< real64 const > > const & pres, diff --git a/src/coreComponents/physicsSolvers/fluidFlow/kernels/compositional/DissipationFluxComputeKernel.hpp b/src/coreComponents/physicsSolvers/fluidFlow/kernels/compositional/DissipationFluxComputeKernel.hpp index 1b4183c8c1a..62df1c74088 100644 --- a/src/coreComponents/physicsSolvers/fluidFlow/kernels/compositional/DissipationFluxComputeKernel.hpp +++ b/src/coreComponents/physicsSolvers/fluidFlow/kernels/compositional/DissipationFluxComputeKernel.hpp @@ -185,6 +185,7 @@ class FluxComputeKernel : public isothermalCompositionalMultiphaseFVMKernels::Fl // // We use the lambda below (called **inside** the phase loop of the base computeFlux) to compute dissipation terms Base::computeFlux( iconn, stack, [&] ( integer const ip, + integer const GEOS_UNUSED_PARAM( useNewGravity ), localIndex const (&k)[2], localIndex const (&seri)[2], localIndex const (&sesri)[2], diff --git a/src/coreComponents/physicsSolvers/fluidFlow/kernels/compositional/FluxComputeKernel.hpp b/src/coreComponents/physicsSolvers/fluidFlow/kernels/compositional/FluxComputeKernel.hpp index d28ecf41b41..acd5c371d7f 100644 --- a/src/coreComponents/physicsSolvers/fluidFlow/kernels/compositional/FluxComputeKernel.hpp +++ b/src/coreComponents/physicsSolvers/fluidFlow/kernels/compositional/FluxComputeKernel.hpp @@ -277,6 +277,7 @@ class FluxComputeKernel : public FluxComputeKernelBase ( m_numPhases, ip, m_kernelFlags.isSet( KernelFlags::CapPressure ), + m_kernelFlags.isSet( KernelFlags::NewGravity ), seri, sesri, sei, trans, dTrans_dPres, @@ -303,6 +304,7 @@ class FluxComputeKernel : public FluxComputeKernelBase ( m_numPhases, ip, m_kernelFlags.isSet( KernelFlags::CapPressure ), + m_kernelFlags.isSet( KernelFlags::NewGravity ), seri, sesri, sei, trans, dTrans_dPres, @@ -329,6 +331,7 @@ class FluxComputeKernel : public FluxComputeKernelBase ( m_numPhases, ip, m_kernelFlags.isSet( KernelFlags::CapPressure ), + m_kernelFlags.isSet( KernelFlags::NewGravity ), seri, sesri, sei, trans, dTrans_dPres, @@ -352,7 +355,8 @@ class FluxComputeKernel : public FluxComputeKernelBase // call the lambda in the phase loop to allow the reuse of the phase fluxes and their derivatives // possible use: assemble the derivatives wrt temperature, and the flux term of the energy equation for this phase - compFluxKernelOp( ip, k, seri, sesri, sei, connectionIndex, + compFluxKernelOp( ip, m_kernelFlags.isSet( KernelFlags::NewGravity ), + k, seri, sesri, sei, connectionIndex, k_up, seri[k_up], sesri[k_up], sei[k_up], potGrad, phaseFlux, dPhaseFlux_dP, dPhaseFlux_dC ); @@ -525,6 +529,7 @@ class FluxComputeKernelFactory string const & dofKey, integer const hasCapPressure, integer const useTotalMassEquation, + integer const useNewGravity, UpwindingParameters upwindingParams, string const & solverName, ElementRegionManager const & elemManager, @@ -547,6 +552,8 @@ class FluxComputeKernelFactory kernelFlags.set( KernelFlags::CapPressure ); if( useTotalMassEquation ) kernelFlags.set( KernelFlags::TotalMassEquation ); + if( useNewGravity ) + kernelFlags.set( KernelFlags::NewGravity ); if( upwindingParams.upwindingScheme == UpwindingScheme::C1PPU && isothermalCompositionalMultiphaseFVMKernelUtilities::epsC1PPU > 0 ) kernelFlags.set( KernelFlags::C1PPU ); diff --git a/src/coreComponents/physicsSolvers/fluidFlow/kernels/compositional/FluxComputeKernelBase.hpp b/src/coreComponents/physicsSolvers/fluidFlow/kernels/compositional/FluxComputeKernelBase.hpp index f36124d6ed6..340d5c7f5f8 100644 --- a/src/coreComponents/physicsSolvers/fluidFlow/kernels/compositional/FluxComputeKernelBase.hpp +++ b/src/coreComponents/physicsSolvers/fluidFlow/kernels/compositional/FluxComputeKernelBase.hpp @@ -47,12 +47,13 @@ enum class KernelFlags CapPressure = 1 << 0, // 1 /// Flag indicating whether total mass equation is formed or not TotalMassEquation = 1 << 1, // 2 + /// Flag indicating whether new gravity treatment is used or not + NewGravity = 1 << 2, // 4 /// Flag indicating whether C1-PPU is used or not - C1PPU = 1 << 2, // 4 + C1PPU = 1 << 3, // 8 /// Flag indicating whether IHU is used or not - IHU = 1 << 3 // 8 + IHU = 1 << 4 // 16 /// Add more flags like that if needed: - // Flag5 = 1 << 4, // 16 // Flag6 = 1 << 5, // 32 // Flag7 = 1 << 6, // 64 // Flag8 = 1 << 7 //128 diff --git a/src/coreComponents/physicsSolvers/fluidFlow/kernels/compositional/IHUPhaseFlux.hpp b/src/coreComponents/physicsSolvers/fluidFlow/kernels/compositional/IHUPhaseFlux.hpp index ec9b3c55783..bce2f8cd261 100644 --- a/src/coreComponents/physicsSolvers/fluidFlow/kernels/compositional/IHUPhaseFlux.hpp +++ b/src/coreComponents/physicsSolvers/fluidFlow/kernels/compositional/IHUPhaseFlux.hpp @@ -65,7 +65,7 @@ upwindMobilityViscous( localIndex const numPhase, ElementViewConst< arrayView3d< real64 const, compflow::USD_PHASE_DC > > const & dPhaseVolFrac, ElementViewConst< arrayView3d< real64 const, constitutive::cappres::USD_CAPPRES > > const & phaseCapPressure, ElementViewConst< arrayView4d< real64 const, constitutive::cappres::USD_CAPPRES_DS > > const & dPhaseCapPressure_dPhaseVolFrac, - integer const capPressureFlag, + integer const hasCapPressure, localIndex & upwindDir, real64 & mobility, real64( &dMobility_dP), @@ -99,7 +99,7 @@ upwindMobilityViscous( localIndex const numPhase, dPhaseVolFrac, phaseCapPressure, dPhaseCapPressure_dPhaseVolFrac, - capPressureFlag, + hasCapPressure, upwindDir ); localIndex const er_up = seri[upwindDir]; @@ -139,7 +139,8 @@ upwindMobilityGravity( localIndex const numPhase, ElementViewConst< arrayView3d< real64 const, compflow::USD_PHASE_DC > > const & dPhaseVolFrac, ElementViewConst< arrayView3d< real64 const, constitutive::cappres::USD_CAPPRES > > const & phaseCapPressure, ElementViewConst< arrayView4d< real64 const, constitutive::cappres::USD_CAPPRES_DS > > const & dPhaseCapPressure_dPhaseVolFrac, - integer const capPressureFlag, + integer const hasCapPressure, + integer const useNewGravity, localIndex & upwindDir, real64 & mobility, real64( &dMobility_dP), @@ -174,7 +175,8 @@ upwindMobilityGravity( localIndex const numPhase, dPhaseVolFrac, phaseCapPressure, dPhaseCapPressure_dPhaseVolFrac, - capPressureFlag, + hasCapPressure, + useNewGravity, upwindDir ); localIndex const er_up = seri[upwindDir]; @@ -213,7 +215,7 @@ upwindMobilityCapillary( localIndex const numPhase, ElementViewConst< arrayView3d< real64 const, compflow::USD_PHASE_DC > > const & dPhaseVolFrac, ElementViewConst< arrayView3d< real64 const, constitutive::cappres::USD_CAPPRES > > const & phaseCapPressure, ElementViewConst< arrayView4d< real64 const, constitutive::cappres::USD_CAPPRES_DS > > const & dPhaseCapPressure_dPhaseVolFrac, - integer const capPressureFlag, + integer const hasCapPressure, localIndex & upwindDir, real64 & mobility, real64( &dMobility_dP), @@ -247,7 +249,7 @@ upwindMobilityCapillary( localIndex const numPhase, dPhaseVolFrac, phaseCapPressure, dPhaseCapPressure_dPhaseVolFrac, - capPressureFlag, + hasCapPressure, upwindDir ); localIndex const er_up = seri[upwindDir]; @@ -290,7 +292,7 @@ computeFractionalFlowViscous( localIndex const numPhase, ElementViewConst< arrayView3d< real64 const, compflow::USD_PHASE_DC > > const & dPhaseVolFrac, ElementViewConst< arrayView3d< real64 const, constitutive::cappres::USD_CAPPRES > > const & phaseCapPressure, ElementViewConst< arrayView4d< real64 const, constitutive::cappres::USD_CAPPRES_DS > > const & dPhaseCapPressure_dPhaseVolFrac, - integer const capPressureFlag, + integer const hasCapPressure, localIndex & k_up_main, real64 & fractionalFlow, real64 ( & dFractionalFlow_dP)[numFluxSupportPoints], @@ -341,7 +343,7 @@ computeFractionalFlowViscous( localIndex const numPhase, dPhaseVolFrac, phaseCapPressure, dPhaseCapPressure_dPhaseVolFrac, - capPressureFlag, + hasCapPressure, k_up, mob, dMob_dP, @@ -405,7 +407,8 @@ computeFractionalFlowGravity( localIndex const numPhase, ElementViewConst< arrayView3d< real64 const, compflow::USD_PHASE_DC > > const & dPhaseVolFrac, ElementViewConst< arrayView3d< real64 const, constitutive::cappres::USD_CAPPRES > > const & phaseCapPressure, ElementViewConst< arrayView4d< real64 const, constitutive::cappres::USD_CAPPRES_DS > > const & dPhaseCapPressure_dPhaseVolFrac, - integer const capPressureFlag, + integer const hasCapPressure, + integer const useNewGravity, localIndex & k_up_main, real64 & fractionalFlow, real64 ( & dFractionalFlow_dP)[numFluxSupportPoints], @@ -455,7 +458,8 @@ computeFractionalFlowGravity( localIndex const numPhase, dPhaseVolFrac, phaseCapPressure, dPhaseCapPressure_dPhaseVolFrac, - capPressureFlag, + hasCapPressure, + useNewGravity, k_up, mob, dMob_dP, @@ -517,7 +521,7 @@ computeFractionalFlowCapillary( localIndex const numPhase, ElementViewConst< arrayView3d< real64 const, compflow::USD_PHASE_DC > > const & dPhaseVolFrac, ElementViewConst< arrayView3d< real64 const, constitutive::cappres::USD_CAPPRES > > const & phaseCapPressure, ElementViewConst< arrayView4d< real64 const, constitutive::cappres::USD_CAPPRES_DS > > const & dPhaseCapPressure_dPhaseVolFrac, - integer const capPressureFlag, + integer const hasCapPressure, localIndex & k_up_main, real64 & fractionalFlow, real64 ( & dFractionalFlow_dP)[numFluxSupportPoints], @@ -565,7 +569,7 @@ computeFractionalFlowCapillary( localIndex const numPhase, dPhaseVolFrac, phaseCapPressure, dPhaseCapPressure_dPhaseVolFrac, - capPressureFlag, + hasCapPressure, k_up, mob, dMob_dP, @@ -658,6 +662,7 @@ struct computePotentialGravity GEOS_HOST_DEVICE static void compute( localIndex const GEOS_UNUSED_PARAM( numPhase ), localIndex const ip, + integer const useNewGravity, localIndex const (&seri)[numFluxSupportPoints], localIndex const (&sesri)[numFluxSupportPoints], localIndex const (&sei)[numFluxSupportPoints], @@ -697,7 +702,46 @@ struct computePotentialGravity } } - //inner loop to get average density + calculateMeanDensity( useNewGravity, ip, seri, sesri, sei, phaseVolFrac, dCompFrac_dCompDens, phaseMassDens, dPhaseMassDens, dProp_dComp, + densMean, dDensMean_dPres, dDensMean_dComp ); + + // compute potential difference MPFA-style + for( localIndex i = 0; i < numFluxSupportPoints; ++i ) + { + localIndex const er = seri[i]; + localIndex const esr = sesri[i]; + localIndex const ei = sei[i]; + + real64 const gravD = transmissibility[i] * gravCoef[er][esr][ei]; + real64 const dGravD_dP = dTrans_dPres[i] * gravCoef[er][esr][ei]; + pot += densMean * gravD; + + // need to add contributions from both cells the mean density depends on + for( localIndex j = 0; j < numFluxSupportPoints; ++j ) + { + dPot_dPres[j] += dDensMean_dPres[j] * gravD + densMean * dGravD_dP; + for( localIndex jc = 0; jc < numComp; ++jc ) + { + dPot_dComp[j][jc] += dDensMean_dComp[j][jc] * gravD; + } + } + } + } + + template< localIndex numComp, localIndex numFluxSupportPoints > + GEOS_HOST_DEVICE + static void calculateMeanDensity( integer const useNewGravity, + localIndex const ip, + localIndex const (&seri)[numFluxSupportPoints], + localIndex const (&sesri)[numFluxSupportPoints], + localIndex const (&sei)[numFluxSupportPoints], + ElementViewConst< arrayView2d< real64 const, compflow::USD_PHASE > > const & phaseVolFrac, + ElementViewConst< arrayView3d< real64 const, compflow::USD_COMP_DC > > const & dCompFrac_dCompDens, + ElementViewConst< arrayView3d< real64 const, constitutive::multifluid::USD_PHASE > > const & phaseMassDens, + ElementViewConst< arrayView4d< real64 const, constitutive::multifluid::USD_PHASE_DC > > const & dPhaseMassDens, + real64 (& dProp_dComp)[numComp], + real64 & densMean, real64 (& dDensMean_dPres)[numFluxSupportPoints], real64 (& dDensMean_dComp)[numFluxSupportPoints][numComp] ) + { integer denom = 0; for( localIndex i = 0; i < numFluxSupportPoints; ++i ) { @@ -706,7 +750,7 @@ struct computePotentialGravity localIndex const ei = sei[i]; bool const phaseExists = (phaseVolFrac[er][esr][ei][ip] > 0); - if( !phaseExists ) + if( useNewGravity && !phaseExists ) { continue; } @@ -742,29 +786,6 @@ struct computePotentialGravity } } } - - // compute potential difference MPFA-style - for( localIndex i = 0; i < numFluxSupportPoints; ++i ) - { - localIndex const er = seri[i]; - localIndex const esr = sesri[i]; - localIndex const ei = sei[i]; - - real64 const gravD = transmissibility[i] * gravCoef[er][esr][ei]; - real64 const dGravD_dP = dTrans_dPres[i] * gravCoef[er][esr][ei]; - pot += densMean * gravD; - - // need to add contributions from both cells the mean density depends on - for( localIndex j = 0; j < numFluxSupportPoints; ++j ) - { - dPot_dPres[j] += dDensMean_dPres[j] * gravD + densMean * dGravD_dP; - for( localIndex jc = 0; jc < numComp; ++jc ) - { - dPot_dComp[j][jc] += dDensMean_dComp[j][jc] * gravD; - } - } - } - } }; @@ -856,7 +877,8 @@ static void computePotentialFluxesGravity( localIndex const numPhase, ElementViewConst< arrayView4d< real64 const, constitutive::multifluid::USD_PHASE_DC > > const & dPhaseMassDens, ElementViewConst< arrayView3d< real64 const, constitutive::cappres::USD_CAPPRES > > const & phaseCapPressure, ElementViewConst< arrayView4d< real64 const, constitutive::cappres::USD_CAPPRES_DS > > const & dPhaseCapPressure_dPhaseVolFrac, - localIndex const capPressureFlag, + localIndex const hasCapPressure, + integer const useNewGravity, localIndex( &k_up), localIndex (&k_up_o), real64 & phaseFlux, @@ -876,6 +898,7 @@ static void computePotentialFluxesGravity( localIndex const numPhase, // UpwindHelpers::computePotentialGravity::compute< numComp, numFluxSupportPoints >( numPhase, ip, + useNewGravity, seri, sesri, sei, @@ -920,7 +943,8 @@ static void computePotentialFluxesGravity( localIndex const numPhase, dPhaseVolFrac, phaseCapPressure, dPhaseCapPressure_dPhaseVolFrac, - capPressureFlag, + hasCapPressure, + useNewGravity, k_up, fflow, dFflow_dP, @@ -940,6 +964,7 @@ static void computePotentialFluxesGravity( localIndex const numPhase, //Fetch pot for phase j!=i defined as \rho_j g dz/dx UpwindHelpers::computePotentialGravity::compute< numComp, numFluxSupportPoints >( numPhase, jp, + useNewGravity, seri, sesri, sei, @@ -986,7 +1011,8 @@ static void computePotentialFluxesGravity( localIndex const numPhase, dPhaseVolFrac, phaseCapPressure, dPhaseCapPressure_dPhaseVolFrac, - capPressureFlag, + hasCapPressure, + useNewGravity, k_up_o, mobOther, dMobOther_dP, @@ -1049,7 +1075,7 @@ static void computePotentialFluxesCapillary( localIndex const numPhase, ElementViewConst< arrayView4d< real64 const, constitutive::multifluid::USD_PHASE_DC > > const & dPhaseMassDens, ElementViewConst< arrayView3d< real64 const, constitutive::cappres::USD_CAPPRES > > const & phaseCapPressure, ElementViewConst< arrayView4d< real64 const, constitutive::cappres::USD_CAPPRES_DS > > const & dPhaseCapPressure_dPhaseVolFrac, - localIndex const capPressureFlag, + localIndex const hasCapPressure, localIndex( &k_up), localIndex (&k_up_o), real64 & phaseFlux, @@ -1110,7 +1136,7 @@ static void computePotentialFluxesCapillary( localIndex const numPhase, dPhaseVolFrac, phaseCapPressure, dPhaseCapPressure_dPhaseVolFrac, - capPressureFlag, + hasCapPressure, k_up, fflow, dFflow_dP, @@ -1174,7 +1200,7 @@ static void computePotentialFluxesCapillary( localIndex const numPhase, dPhaseVolFrac, phaseCapPressure, dPhaseCapPressure_dPhaseVolFrac, - capPressureFlag, + hasCapPressure, k_up_o, mobOther, dMobOther_dP, @@ -1261,7 +1287,7 @@ class UpwindScheme ElementViewConst< arrayView3d< real64 const, compflow::USD_PHASE_DC > > const & dPhaseVolFrac, ElementViewConst< arrayView3d< real64 const, constitutive::cappres::USD_CAPPRES > > const & phaseCapPressure, ElementViewConst< arrayView4d< real64 const, constitutive::cappres::USD_CAPPRES_DS > > const & dPhaseCapPressure_dPhaseVolFrac, - integer const capPressureFlag, + integer const hasCapPressure, localIndex & upwindDir ) { @@ -1286,7 +1312,7 @@ class UpwindScheme dPhaseVolFrac, phaseCapPressure, dPhaseCapPressure_dPhaseVolFrac, - capPressureFlag, + hasCapPressure, pot ); //all definition has been changed to fit pot>0 => first cell is upstream @@ -1314,7 +1340,8 @@ class UpwindScheme ElementViewConst< arrayView3d< real64 const, compflow::USD_PHASE_DC > > const & dPhaseVolFrac, ElementViewConst< arrayView3d< real64 const, constitutive::cappres::USD_CAPPRES > > const & phaseCapPressure, ElementViewConst< arrayView4d< real64 const, constitutive::cappres::USD_CAPPRES_DS > > const & dPhaseCapPressure_dPhaseVolFrac, - integer const capPressureFlag, + integer const hasCapPressure, + integer const useNewGravity, localIndex & upwindDir ) { @@ -1340,7 +1367,8 @@ class UpwindScheme dPhaseVolFrac, phaseCapPressure, dPhaseCapPressure_dPhaseVolFrac, - capPressureFlag, + hasCapPressure, + useNewGravity, pot ); //all definition has been changed to fit pot>0 => first cell is upstream @@ -1368,7 +1396,7 @@ class UpwindScheme ElementViewConst< arrayView3d< real64 const, compflow::USD_PHASE_DC > > const & dPhaseVolFrac, ElementViewConst< arrayView3d< real64 const, constitutive::cappres::USD_CAPPRES > > const & phaseCapPressure, ElementViewConst< arrayView4d< real64 const, constitutive::cappres::USD_CAPPRES_DS > > const & dPhaseCapPressure_dPhaseVolFrac, - integer const capPressureFlag, + integer const hasCapPressure, localIndex & upwindDir ) { @@ -1393,7 +1421,7 @@ class UpwindScheme dPhaseVolFrac, phaseCapPressure, dPhaseCapPressure_dPhaseVolFrac, - capPressureFlag, + hasCapPressure, pot ); //all definition has been changed to fit pot>0 => first cell is upstream @@ -1486,7 +1514,7 @@ class HybridUpwind : public UpwindScheme ElementViewConst< arrayView3d< real64 const, compflow::USD_PHASE_DC > > const & dPhaseVolFrac, ElementViewConst< arrayView3d< real64 const, constitutive::cappres::USD_CAPPRES > > const & phaseCapPressure, ElementViewConst< arrayView4d< real64 const, constitutive::cappres::USD_CAPPRES_DS > > const & dPhaseCapPressure_dPhaseVolFrac, - integer const GEOS_UNUSED_PARAM( capPressureFlag ), + integer const GEOS_UNUSED_PARAM( hasCapPressure ), real64 & potential ) { @@ -1538,7 +1566,8 @@ class HybridUpwind : public UpwindScheme ElementViewConst< arrayView3d< real64 const, compflow::USD_PHASE_DC > > const & dPhaseVolFrac, ElementViewConst< arrayView3d< real64 const, constitutive::cappres::USD_CAPPRES > > const & phaseCapPressure, ElementViewConst< arrayView4d< real64 const, constitutive::cappres::USD_CAPPRES_DS > > const & dPhaseCapPressure_dPhaseVolFrac, - integer const GEOS_UNUSED_PARAM( capPressureFlag ), + integer const GEOS_UNUSED_PARAM( hasCapPressure ), + integer const useNewGravity, real64 & potential ) { @@ -1558,6 +1587,7 @@ class HybridUpwind : public UpwindScheme UpwindHelpers::computePotentialGravity::compute< numComp, numFluxSupportPoints >( numPhase, ipp, + useNewGravity, seri, sesri, sei, @@ -1602,7 +1632,7 @@ class HybridUpwind : public UpwindScheme ElementViewConst< arrayView3d< real64 const, compflow::USD_PHASE_DC > > const & dPhaseVolFrac, ElementViewConst< arrayView3d< real64 const, constitutive::cappres::USD_CAPPRES > > const & phaseCapPressure, ElementViewConst< arrayView4d< real64 const, constitutive::cappres::USD_CAPPRES_DS > > const & dPhaseCapPressure_dPhaseVolFrac, - integer const GEOS_UNUSED_PARAM( capPressureFlag ), + integer const GEOS_UNUSED_PARAM( hasCapPressure ), real64 & potential ) { @@ -1685,6 +1715,7 @@ struct IHUPhaseFlux compute( integer const numPhase, integer const ip, integer const hasCapPressure, + integer const useNewGravity, localIndex const ( &seri )[numFluxSupportPoints], localIndex const ( &sesri )[numFluxSupportPoints], localIndex const ( &sei )[numFluxSupportPoints], @@ -1732,7 +1763,7 @@ struct IHUPhaseFlux for( integer jp = 0; jp < numPhase; ++jp ) { - PPUPhaseFlux::compute( numPhase, jp, hasCapPressure, + PPUPhaseFlux::compute( numPhase, jp, hasCapPressure, useNewGravity, seri, sesri, sei, trans, dTrans_dPres, pres, gravCoef, @@ -1886,6 +1917,7 @@ struct IHUPhaseFlux phaseCapPressure, dPhaseCapPressure_dPhaseVolFrac, hasCapPressure, + useNewGravity, k_up_g, k_up_og, gravitationalPhaseFlux, diff --git a/src/coreComponents/physicsSolvers/fluidFlow/kernels/compositional/PPUPhaseFlux.hpp b/src/coreComponents/physicsSolvers/fluidFlow/kernels/compositional/PPUPhaseFlux.hpp index c8ac5256d98..58ba51a9f94 100644 --- a/src/coreComponents/physicsSolvers/fluidFlow/kernels/compositional/PPUPhaseFlux.hpp +++ b/src/coreComponents/physicsSolvers/fluidFlow/kernels/compositional/PPUPhaseFlux.hpp @@ -75,6 +75,7 @@ struct PPUPhaseFlux compute( integer const numPhase, integer const ip, integer const hasCapPressure, + integer const useNewGravity, localIndex const ( &seri )[numFluxSupportPoints], localIndex const ( &sesri )[numFluxSupportPoints], localIndex const ( &sei )[numFluxSupportPoints], @@ -106,7 +107,7 @@ struct PPUPhaseFlux real64 dPresGrad_dC[numFluxSupportPoints][numComp]{}; real64 dGravHead_dP[numFluxSupportPoints]{}; real64 dGravHead_dC[numFluxSupportPoints][numComp]{}; - PotGrad::compute< numComp, numFluxSupportPoints >( numPhase, ip, hasCapPressure, seri, sesri, sei, trans, dTrans_dPres, pres, + PotGrad::compute< numComp, numFluxSupportPoints >( numPhase, ip, hasCapPressure, useNewGravity, seri, sesri, sei, trans, dTrans_dPres, pres, gravCoef, phaseVolFrac, dPhaseVolFrac, dCompFrac_dCompDens, phaseMassDens, dPhaseMassDens, phaseCapPressure, dPhaseCapPressure_dPhaseVolFrac, potGrad, dPresGrad_dP, dPresGrad_dC, dGravHead_dP, dGravHead_dC ); diff --git a/src/coreComponents/physicsSolvers/fluidFlow/kernels/compositional/PotGrad.hpp b/src/coreComponents/physicsSolvers/fluidFlow/kernels/compositional/PotGrad.hpp index 14f702792db..138175d5c89 100644 --- a/src/coreComponents/physicsSolvers/fluidFlow/kernels/compositional/PotGrad.hpp +++ b/src/coreComponents/physicsSolvers/fluidFlow/kernels/compositional/PotGrad.hpp @@ -46,6 +46,7 @@ struct PotGrad compute ( integer const numPhase, integer const ip, integer const hasCapPressure, + integer const useNewGravity, localIndex const ( &seri )[numFluxSupportPoints], localIndex const ( &sesri )[numFluxSupportPoints], localIndex const ( &sei )[numFluxSupportPoints], @@ -87,53 +88,7 @@ struct PotGrad real64 gravHead = 0.0; real64 dCapPressure_dC[numComp]{}; - real64 dProp_dC[numComp]{}; - - // calculate quantities on primary connected cells - integer denom = 0; - for( integer i = 0; i < numFluxSupportPoints; ++i ) - { - localIndex const er = seri[i]; - localIndex const esr = sesri[i]; - localIndex const ei = sei[i]; - - bool const phaseExists = (phaseVolFrac[er][esr][ei][ip] > 0); - if( !phaseExists ) - { - continue; - } - - // density - real64 const density = phaseMassDens[er][esr][ei][0][ip]; - real64 const dDens_dP = dPhaseMassDens[er][esr][ei][0][ip][Deriv::dP]; - - applyChainRule( numComp, - dCompFrac_dCompDens[er][esr][ei], - dPhaseMassDens[er][esr][ei][0][ip], - dProp_dC, - Deriv::dC ); - - // average density and derivatives - densMean += density; - dDensMean_dP[i] = dDens_dP; - for( integer jc = 0; jc < numComp; ++jc ) - { - dDensMean_dC[i][jc] = dProp_dC[jc]; - } - denom++; - } - if( denom > 1 ) - { - densMean /= denom; - for( integer i = 0; i < numFluxSupportPoints; ++i ) - { - dDensMean_dP[i] /= denom; - for( integer jc = 0; jc < numComp; ++jc ) - { - dDensMean_dC[i][jc] /= denom; - } - } - } + calculateMeanDensity( useNewGravity, ip, seri, sesri, sei, phaseVolFrac, dCompFrac_dCompDens, phaseMassDens, dPhaseMassDens, densMean, dDensMean_dP, dDensMean_dC ); /// compute the TPFA potential difference for( integer i = 0; i < numFluxSupportPoints; i++ ) @@ -199,6 +154,68 @@ struct PotGrad } + template< integer numComp, integer numFluxSupportPoints > + GEOS_HOST_DEVICE + static void + calculateMeanDensity( integer const useNewGravity, + integer const ip, + localIndex const ( &seri )[numFluxSupportPoints], + localIndex const ( &sesri )[numFluxSupportPoints], + localIndex const ( &sei )[numFluxSupportPoints], + ElementViewConst< arrayView2d< real64 const, compflow::USD_PHASE > > const & phaseVolFrac, + ElementViewConst< arrayView3d< real64 const, compflow::USD_COMP_DC > > const & dCompFrac_dCompDens, + ElementViewConst< arrayView3d< real64 const, constitutive::multifluid::USD_PHASE > > const & phaseMassDens, + ElementViewConst< arrayView4d< real64 const, constitutive::multifluid::USD_PHASE_DC > > const & dPhaseMassDens, + real64 & densMean, real64 ( & dDensMean_dP)[numFluxSupportPoints], real64 ( & dDensMean_dC )[numFluxSupportPoints][numComp] ) + { + real64 dDens_dC[numComp]{}; + + integer denom = 0; + for( integer i = 0; i < numFluxSupportPoints; ++i ) + { + localIndex const er = seri[i]; + localIndex const esr = sesri[i]; + localIndex const ei = sei[i]; + + bool const phaseExists = (phaseVolFrac[er][esr][ei][ip] > 0); + if( useNewGravity && !phaseExists ) + { + continue; + } + + // density + real64 const density = phaseMassDens[er][esr][ei][0][ip]; + real64 const dDens_dP = dPhaseMassDens[er][esr][ei][0][ip][Deriv::dP]; + + applyChainRule( numComp, + dCompFrac_dCompDens[er][esr][ei], + dPhaseMassDens[er][esr][ei][0][ip], + dDens_dC, + Deriv::dC ); + + // average density and derivatives + densMean += density; + dDensMean_dP[i] = dDens_dP; + for( integer jc = 0; jc < numComp; ++jc ) + { + dDensMean_dC[i][jc] = dDens_dC[jc]; + } + denom++; + } + if( denom > 1 ) + { + densMean /= denom; + for( integer i = 0; i < numFluxSupportPoints; ++i ) + { + dDensMean_dP[i] /= denom; + for( integer jc = 0; jc < numComp; ++jc ) + { + dDensMean_dC[i][jc] /= denom; + } + } + } + } + }; } // namespace isothermalCompositionalMultiPhaseFVMKernelUtilities diff --git a/src/coreComponents/physicsSolvers/fluidFlow/kernels/compositional/StabilizedFluxComputeKernel.hpp b/src/coreComponents/physicsSolvers/fluidFlow/kernels/compositional/StabilizedFluxComputeKernel.hpp index de9e92c8b80..71801e9ee59 100644 --- a/src/coreComponents/physicsSolvers/fluidFlow/kernels/compositional/StabilizedFluxComputeKernel.hpp +++ b/src/coreComponents/physicsSolvers/fluidFlow/kernels/compositional/StabilizedFluxComputeKernel.hpp @@ -193,6 +193,7 @@ class FluxComputeKernel : public isothermalCompositionalMultiphaseFVMKernels::Fl // // We use the lambda below (called **inside** the phase loop of the base computeFlux) to compute stabilization terms Base::computeFlux( iconn, stack, [&] ( integer const ip, + integer const GEOS_UNUSED_PARAM( useNewGravity ), localIndex const (&k)[2], localIndex const (&seri)[2], localIndex const (&sesri)[2], diff --git a/src/coreComponents/physicsSolvers/fluidFlow/kernels/compositional/ThermalFluxComputeKernel.hpp b/src/coreComponents/physicsSolvers/fluidFlow/kernels/compositional/ThermalFluxComputeKernel.hpp index 17395f62324..8d7190cc70b 100644 --- a/src/coreComponents/physicsSolvers/fluidFlow/kernels/compositional/ThermalFluxComputeKernel.hpp +++ b/src/coreComponents/physicsSolvers/fluidFlow/kernels/compositional/ThermalFluxComputeKernel.hpp @@ -197,6 +197,7 @@ class FluxComputeKernel : public isothermalCompositionalMultiphaseFVMKernels::Fl // such as potGrad, phaseFlux, and the indices of the upwind cell // We use the lambda below (called **inside** the phase loop of the base computeFlux) to access these variables Base::computeFlux( iconn, stack, [&] ( integer const ip, + integer const useNewGravity, localIndex const (&k)[2], localIndex const (&seri)[2], localIndex const (&sesri)[2], @@ -234,7 +235,7 @@ class FluxComputeKernel : public isothermalCompositionalMultiphaseFVMKernels::Fl localIndex const ei = sei[i]; bool const phaseExists = (m_phaseVolFrac[er_up][esr_up][ei_up][ip] > 0); - if( !phaseExists ) + if( useNewGravity && !phaseExists ) { continue; } From eef8de43ccee83c4dfac15b313a4c783ee69199f Mon Sep 17 00:00:00 2001 From: "Victor A. P. Magri" <50467563+victorapm@users.noreply.github.com> Date: Sat, 7 Dec 2024 07:57:15 -0800 Subject: [PATCH 2/2] feat: IO timers (#3480) --- .integrated_tests.yaml | 2 +- BASELINE_NOTES.md | 4 + .../fileIO/Outputs/BlueprintOutput.cpp | 84 +++++++++++-------- .../fileIO/Outputs/BlueprintOutput.hpp | 6 ++ .../fileIO/Outputs/ChomboIO.cpp | 14 ++++ .../fileIO/Outputs/ChomboIO.hpp | 6 ++ .../fileIO/Outputs/OutputBase.cpp | 21 +++++ .../fileIO/Outputs/OutputBase.hpp | 57 ++++++++++++- .../fileIO/Outputs/PythonOutput.cpp | 14 ++++ .../fileIO/Outputs/PythonOutput.hpp | 3 + .../fileIO/Outputs/RestartOutput.cpp | 29 +++++-- .../fileIO/Outputs/RestartOutput.hpp | 6 ++ .../fileIO/Outputs/SiloOutput.cpp | 68 +++++++++------ .../fileIO/Outputs/SiloOutput.hpp | 6 ++ .../fileIO/Outputs/TimeHistoryOutput.cpp | 28 ++++++- .../fileIO/Outputs/TimeHistoryOutput.hpp | 6 ++ .../fileIO/Outputs/VTKOutput.cpp | 38 +++++++-- .../fileIO/Outputs/VTKOutput.hpp | 9 ++ 18 files changed, 319 insertions(+), 82 deletions(-) diff --git a/.integrated_tests.yaml b/.integrated_tests.yaml index 0ce10cbe33b..3aa31a7e4cb 100644 --- a/.integrated_tests.yaml +++ b/.integrated_tests.yaml @@ -1,6 +1,6 @@ baselines: bucket: geosx - baseline: integratedTests/baseline_integratedTests-pr3467-9212-976cc3b + baseline: integratedTests/baseline_integratedTests-pr3480-9217-caaecb8 allow_fail: all: '' streak: '' diff --git a/BASELINE_NOTES.md b/BASELINE_NOTES.md index 8f344877826..acbeb39e382 100644 --- a/BASELINE_NOTES.md +++ b/BASELINE_NOTES.md @@ -6,6 +6,10 @@ This file is designed to track changes to the integrated test baselines. Any developer who updates the baseline ID in the .integrated_tests.yaml file is expected to create an entry in this file with the pull request number, date, and their justification for rebaselining. These notes should be in reverse-chronological order, and use the following time format: (YYYY-MM-DD). +PR #3480 (2024-12-06) +===================== +Add "logLevel" parameter under /Problem/Outputs in baseline files + PR #3361 (2024-12-03) ===================== Revert default gravity treatment to old version. Make the way introduced in #3337 optional. diff --git a/src/coreComponents/fileIO/Outputs/BlueprintOutput.cpp b/src/coreComponents/fileIO/Outputs/BlueprintOutput.cpp index 21457e7796e..cad26fcbb54 100644 --- a/src/coreComponents/fileIO/Outputs/BlueprintOutput.cpp +++ b/src/coreComponents/fileIO/Outputs/BlueprintOutput.cpp @@ -129,53 +129,57 @@ BlueprintOutput::BlueprintOutput( string const & name, } /////////////////////////////////////////////////////////////////////////////////////////////////// -bool BlueprintOutput::execute( real64 const time, - real64 const, - integer const cycle, - integer const, - real64 const, +bool BlueprintOutput::execute( real64 const time_n, + real64 const GEOS_UNUSED_PARAM( dt ), + integer const cycleNumber, + integer const GEOS_UNUSED_PARAM( eventCounter ), + real64 const GEOS_UNUSED_PARAM( eventProgress ), DomainPartition & domain ) { GEOS_MARK_FUNCTION; - MeshLevel const & meshLevel = domain.getMeshBody( 0 ).getBaseDiscretization(); + { + Timer timer( m_outputTimer ); - conduit::Node meshRoot; - conduit::Node & mesh = meshRoot[ "mesh" ]; - conduit::Node & coordset = mesh[ "coordsets/nodes" ]; - conduit::Node & topologies = mesh[ "topologies" ]; + MeshLevel const & meshLevel = domain.getMeshBody( 0 ).getBaseDiscretization(); - mesh[ "state/time" ] = time; - mesh[ "state/cycle" ] = cycle; + conduit::Node meshRoot; + conduit::Node & mesh = meshRoot[ "mesh" ]; + conduit::Node & coordset = mesh[ "coordsets/nodes" ]; + conduit::Node & topologies = mesh[ "topologies" ]; - addNodalData( meshLevel.getNodeManager(), coordset, topologies, mesh[ "fields" ] ); + mesh[ "state/time" ] = time_n; + mesh[ "state/cycle" ] = cycleNumber; - dataRepository::Group averagedElementData( "averagedElementData", this ); - addElementData( meshLevel.getElemManager(), coordset, topologies, mesh[ "fields" ], averagedElementData ); + addNodalData( meshLevel.getNodeManager(), coordset, topologies, mesh[ "fields" ] ); - /// The Blueprint will complain if the fields node is present but empty. - if( mesh[ "fields" ].number_of_children() == 0 ) - { - mesh.remove( "fields" ); - } + dataRepository::Group averagedElementData( "averagedElementData", this ); + addElementData( meshLevel.getElemManager(), coordset, topologies, mesh[ "fields" ], averagedElementData ); + + /// The Blueprint will complain if the fields node is present but empty. + if( mesh[ "fields" ].number_of_children() == 0 ) + { + mesh.remove( "fields" ); + } - /// Verify that the mesh conforms to the Blueprint. - conduit::Node info; - GEOS_ASSERT_MSG( conduit::blueprint::verify( "mesh", meshRoot, info ), info.to_json() ); + /// Verify that the mesh conforms to the Blueprint. + conduit::Node info; + GEOS_ASSERT_MSG( conduit::blueprint::verify( "mesh", meshRoot, info ), info.to_json() ); - /// Generate the Blueprint index. - conduit::Node fileRoot; - conduit::Node & index = fileRoot[ "blueprint_index/mesh" ]; - conduit::blueprint::mesh::generate_index( mesh, "mesh", MpiWrapper::commSize(), index ); + /// Generate the Blueprint index. + conduit::Node fileRoot; + conduit::Node & index = fileRoot[ "blueprint_index/mesh" ]; + conduit::blueprint::mesh::generate_index( mesh, "mesh", MpiWrapper::commSize(), index ); - /// Verify that the index conforms to the Blueprint. - info.reset(); - GEOS_ASSERT_MSG( conduit::blueprint::mesh::index::verify( index, info ), info.to_json() ); + /// Verify that the index conforms to the Blueprint. + info.reset(); + GEOS_ASSERT_MSG( conduit::blueprint::mesh::index::verify( index, info ), info.to_json() ); - /// Write out the root index file, then write out the mesh. - string const completePath = GEOS_FMT( "{}/blueprintFiles/cycle_{:07}", OutputBase::getOutputDirectory(), cycle ); - string const filePathForRank = dataRepository::writeRootFile( fileRoot, completePath ); - conduit::relay::io::save( meshRoot, filePathForRank, "hdf5" ); + /// Write out the root index file, then write out the mesh. + string const completePath = GEOS_FMT( "{}/blueprintFiles/cycle_{:07}", OutputBase::getOutputDirectory(), cycleNumber ); + string const filePathForRank = dataRepository::writeRootFile( fileRoot, completePath ); + conduit::relay::io::save( meshRoot, filePathForRank, "hdf5" ); + } return false; } @@ -307,7 +311,19 @@ void BlueprintOutput::writeOutConstitutiveData( dataRepository::Group const & co } ); } +namespace logInfo +{ +struct BlueprintOutputTimer : public OutputTimerBase +{ + std::string_view getDescription() const override { return "Blueprint output timing"; } +}; +} +logInfo::OutputTimerBase const & BlueprintOutput::getTimerCategory() const +{ + static logInfo::BlueprintOutputTimer timer; + return timer; +} REGISTER_CATALOG_ENTRY( OutputBase, BlueprintOutput, string const &, dataRepository::Group * const ) diff --git a/src/coreComponents/fileIO/Outputs/BlueprintOutput.hpp b/src/coreComponents/fileIO/Outputs/BlueprintOutput.hpp index 7d17570d47d..fd963698c3c 100644 --- a/src/coreComponents/fileIO/Outputs/BlueprintOutput.hpp +++ b/src/coreComponents/fileIO/Outputs/BlueprintOutput.hpp @@ -36,6 +36,12 @@ class ElementRegionManager; */ class BlueprintOutput : public OutputBase { +protected: + /** + * @copydoc OutputBase::getTimerCategory + */ + logInfo::OutputTimerBase const & getTimerCategory() const override; + public: /** diff --git a/src/coreComponents/fileIO/Outputs/ChomboIO.cpp b/src/coreComponents/fileIO/Outputs/ChomboIO.cpp index 674f1e4eace..26f73a9eb43 100644 --- a/src/coreComponents/fileIO/Outputs/ChomboIO.cpp +++ b/src/coreComponents/fileIO/Outputs/ChomboIO.cpp @@ -30,6 +30,20 @@ namespace geos using namespace dataRepository; +namespace logInfo +{ +struct ChomboOutputTimer : public OutputTimerBase +{ + std::string_view getDescription() const override { return "Chombo output timing"; } +}; +} + +logInfo::OutputTimerBase const & ChomboIO::getTimerCategory() const +{ + static logInfo::ChomboOutputTimer timer; + return timer; +} + ChomboIO::ChomboIO( string const & name, Group * const parent ): OutputBase( name, parent ), m_coupler( nullptr ), diff --git a/src/coreComponents/fileIO/Outputs/ChomboIO.hpp b/src/coreComponents/fileIO/Outputs/ChomboIO.hpp index c483f066968..c9e60cf2f79 100644 --- a/src/coreComponents/fileIO/Outputs/ChomboIO.hpp +++ b/src/coreComponents/fileIO/Outputs/ChomboIO.hpp @@ -89,6 +89,12 @@ class ChomboIO final : public OutputBase } viewKeys; /// @endcond +protected: + /** + * @copydoc OutputBase::getTimerCategory + */ + logInfo::OutputTimerBase const & getTimerCategory() const override; + private: ChomboCoupler * m_coupler; string m_outputPath; diff --git a/src/coreComponents/fileIO/Outputs/OutputBase.cpp b/src/coreComponents/fileIO/Outputs/OutputBase.cpp index 0733dd38cc5..54dcab668e6 100644 --- a/src/coreComponents/fileIO/Outputs/OutputBase.cpp +++ b/src/coreComponents/fileIO/Outputs/OutputBase.cpp @@ -29,6 +29,7 @@ using namespace dataRepository; OutputBase::OutputBase( string const & name, Group * const parent ): ExecutableGroup( name, parent ), + m_outputTimer(), m_childDirectory(), m_parallelThreads( 1 ) { @@ -43,6 +44,8 @@ OutputBase::OutputBase( string const & name, setInputFlag( InputFlags::OPTIONAL ). setDescription( "Number of plot files." ); + // Add the Timers log level + addLogLevel< logInfo::OutputTimers >(); } OutputBase::~OutputBase() @@ -106,5 +109,23 @@ void OutputBase::setupDirectoryStructure() } } +void OutputBase::cleanup( real64 const GEOS_UNUSED_PARAM( time_n ), + integer const GEOS_UNUSED_PARAM( cycleNumber ), + integer const GEOS_UNUSED_PARAM( eventCounter ), + real64 const GEOS_UNUSED_PARAM( eventProgress ), + DomainPartition & GEOS_UNUSED_PARAM( domain ) ) +{ + // Report timing statistics + real64 const time = std::chrono::duration< double >( m_outputTimer ).count(); + real64 const minTime = MpiWrapper::min( time ); + real64 const maxTime = MpiWrapper::max( time ); + if( maxTime > 0 ) + { + GEOS_LOG_LEVEL_INFO_RANK_0( logInfo::OutputTimers, + GEOS_FMT( "{}: file writing time = {} s (min), {} s (max)", + getName(), minTime, maxTime ) ); + } +} + } /* namespace geos */ diff --git a/src/coreComponents/fileIO/Outputs/OutputBase.hpp b/src/coreComponents/fileIO/Outputs/OutputBase.hpp index f4129433f7d..ae957e22213 100644 --- a/src/coreComponents/fileIO/Outputs/OutputBase.hpp +++ b/src/coreComponents/fileIO/Outputs/OutputBase.hpp @@ -21,11 +21,50 @@ #include "dataRepository/Group.hpp" #include "dataRepository/ExecutableGroup.hpp" - +#include "dataRepository/LogLevelsInfo.hpp" // For logInfo namespace +#include "common/Timer.hpp" namespace geos { +namespace logInfo +{ +/** + * @brief Base timer category for output operations + * @details Provides configuration for logging output operation timing information + */ +struct OutputTimers +{ + /** + * @brief Get the description of this timer + * @return String view containing the timer description + */ + static std::string_view getDescription() { return "Output timing information"; } + + /** + * @brief Get the minimum log level for this timer + * @return Integer representing the minimum log level + */ + static constexpr int getMinLogLevel() { return 1; } +}; + +/** + * @brief Base interface for specific output type timers + * @details Each output type (VTK, Silo, etc.) implements this interface to provide + * its own timing category. This is used in conjunction with OutputTimers: + * - OutputTimerBase: For polymorphic behavior in derived output classes + * - OutputTimers: For the general output timing logging infrastructure + */ +struct OutputTimerBase +{ + /** + * @brief Get the description of this timer + * @return String view containing the timer description + */ + virtual std::string_view getDescription() const = 0; +}; +} + /** * @class OutputBase * @@ -102,6 +141,22 @@ class OutputBase : public ExecutableGroup **/ virtual void initializePreSubGroups() override; + /// Timer used to track duration of file writing operations for this specific output type + std::chrono::system_clock::duration m_outputTimer; + + /** + * @brief Get the timer category for this output type + * @return Reference to the output timer base for timing statistics + */ + virtual logInfo::OutputTimerBase const & getTimerCategory() const = 0; + + /// @copydoc geos::ExecutableGroup::cleanup + virtual void cleanup( real64 const time_n, + integer const cycleNumber, + integer const eventCounter, + real64 const eventProgress, + DomainPartition & domain ) override; + private: string m_childDirectory; integer m_parallelThreads; diff --git a/src/coreComponents/fileIO/Outputs/PythonOutput.cpp b/src/coreComponents/fileIO/Outputs/PythonOutput.cpp index e78f8835a06..5891fe55e26 100644 --- a/src/coreComponents/fileIO/Outputs/PythonOutput.cpp +++ b/src/coreComponents/fileIO/Outputs/PythonOutput.cpp @@ -18,6 +18,20 @@ namespace geos { +namespace logInfo +{ +struct PythonOutputTimer : public OutputTimerBase +{ + std::string_view getDescription() const override { return "Python output timing"; } +}; +} + +logInfo::OutputTimerBase const & PythonOutput::getTimerCategory() const +{ + static logInfo::PythonOutputTimer timer; + return timer; +} + REGISTER_CATALOG_ENTRY( OutputBase, PythonOutput, string const &, dataRepository::Group * const ) } // namespace geos diff --git a/src/coreComponents/fileIO/Outputs/PythonOutput.hpp b/src/coreComponents/fileIO/Outputs/PythonOutput.hpp index 0e6c3cebb85..a73b6781654 100644 --- a/src/coreComponents/fileIO/Outputs/PythonOutput.hpp +++ b/src/coreComponents/fileIO/Outputs/PythonOutput.hpp @@ -75,6 +75,9 @@ class PythonOutput : public OutputBase GEOS_UNUSED_VAR( domain ); return true; } + +protected: + logInfo::OutputTimerBase const & getTimerCategory() const override; }; diff --git a/src/coreComponents/fileIO/Outputs/RestartOutput.cpp b/src/coreComponents/fileIO/Outputs/RestartOutput.cpp index 6fb7a2c6b19..6eb07281cf6 100644 --- a/src/coreComponents/fileIO/Outputs/RestartOutput.cpp +++ b/src/coreComponents/fileIO/Outputs/RestartOutput.cpp @@ -24,6 +24,14 @@ namespace geos using namespace dataRepository; +namespace logInfo +{ +struct RestartOutputTimer : public OutputTimerBase +{ + std::string_view getDescription() const override { return "Restart output timing"; } +}; +} + RestartOutput::RestartOutput( string const & name, Group * const parent ): OutputBase( name, parent ) @@ -41,19 +49,24 @@ bool RestartOutput::execute( real64 const GEOS_UNUSED_PARAM( time_n ), { GEOS_MARK_FUNCTION; - Group & rootGroup = this->getGroupByPath( "/Problem" ); + { + Timer timer( m_outputTimer ); - // Ignoring the eventProgress indicator for now to be compliant with the integrated test repo - // integer const eventProgressPercent = static_cast(eventProgress * 100.0); - string const fileName = GEOS_FMT( "{}_restart_{:09}", getFileNameRoot(), cycleNumber ); - - rootGroup.prepareToWrite(); - writeTree( joinPath( OutputBase::getOutputDirectory(), fileName ), *(rootGroup.getConduitNode().parent()) ); - rootGroup.finishWriting(); + Group & rootGroup = this->getGroupByPath( "/Problem" ); + string const fileName = GEOS_FMT( "{}_restart_{:09}", getFileNameRoot(), cycleNumber ); + rootGroup.prepareToWrite(); + writeTree( joinPath( OutputBase::getOutputDirectory(), fileName ), *(rootGroup.getConduitNode().parent()) ); + rootGroup.finishWriting(); + } return false; } +logInfo::OutputTimerBase const & RestartOutput::getTimerCategory() const +{ + static logInfo::RestartOutputTimer timer; + return timer; +} REGISTER_CATALOG_ENTRY( OutputBase, RestartOutput, string const &, Group * const ) } /* namespace geos */ diff --git a/src/coreComponents/fileIO/Outputs/RestartOutput.hpp b/src/coreComponents/fileIO/Outputs/RestartOutput.hpp index a2a12be781c..11c8715ed76 100644 --- a/src/coreComponents/fileIO/Outputs/RestartOutput.hpp +++ b/src/coreComponents/fileIO/Outputs/RestartOutput.hpp @@ -78,6 +78,12 @@ class RestartOutput : public OutputBase dataRepository::ViewKey writeFEMFaces = { "writeFEMFaces" }; } viewKeys; /// @endcond + +protected: + /** + * @copydoc OutputBase::getTimerCategory + */ + logInfo::OutputTimerBase const & getTimerCategory() const override; }; diff --git a/src/coreComponents/fileIO/Outputs/SiloOutput.cpp b/src/coreComponents/fileIO/Outputs/SiloOutput.cpp index 756facc2da1..ba03ceead5f 100644 --- a/src/coreComponents/fileIO/Outputs/SiloOutput.cpp +++ b/src/coreComponents/fileIO/Outputs/SiloOutput.cpp @@ -28,6 +28,20 @@ namespace geos using namespace dataRepository; +namespace logInfo +{ +struct SiloOutputTimer : public OutputTimerBase +{ + std::string_view getDescription() const override { return "Silo output timing"; } +}; +} + +logInfo::OutputTimerBase const & SiloOutput::getTimerCategory() const +{ + static logInfo::SiloOutputTimer timer; + return timer; +} + SiloOutput::SiloOutput( string const & name, Group * const parent ): OutputBase( name, parent ), @@ -122,31 +136,35 @@ bool SiloOutput::execute( real64 const time_n, { GEOS_MARK_FUNCTION; - SiloFile silo; - - int const size = MpiWrapper::commSize( MPI_COMM_GEOS ); - int const rank = MpiWrapper::commRank( MPI_COMM_GEOS ); - MpiWrapper::barrier( MPI_COMM_GEOS ); - - integer const numFiles = parallelThreads() == 0 ? size : parallelThreads(); - - // TODO set this during initialization - // silo.setOutputDirectory( getGlobalState().getCommandLineOptions().outputDirectory ), - silo.setOutputDirectory( getOutputDirectory() ), - silo.setPlotLevel( m_plotLevel ); - silo.setWriteEdgeMesh( m_writeEdgeMesh ); - silo.setWriteFaceMesh( m_writeFaceMesh ); - silo.setWriteCellElementMesh( m_writeCellElementMesh ); - silo.setWriteFaceElementMesh( m_writeFaceElementMesh ); - silo.setOnlyPlotSpecifiedFieldNamesFlag( m_onlyPlotSpecifiedFieldNames ); - silo.setFieldNames( m_fieldNames.toViewConst() ); - silo.setPlotFileRoot( m_plotFileRoot ); - silo.initialize( numFiles ); - silo.waitForBatonWrite( rank, cycleNumber, eventCounter, false ); - silo.writeDomainPartition( domain, cycleNumber, time_n + dt * eventProgress, 0 ); - silo.handOffBaton(); - silo.clearEmptiesFromMultiObjects( cycleNumber ); - silo.finish(); + { + Timer timer( m_outputTimer ); + + SiloFile silo; + + int const size = MpiWrapper::commSize( MPI_COMM_GEOS ); + int const rank = MpiWrapper::commRank( MPI_COMM_GEOS ); + MpiWrapper::barrier( MPI_COMM_GEOS ); + + integer const numFiles = parallelThreads() == 0 ? size : parallelThreads(); + + // TODO set this during initialization + // silo.setOutputDirectory( getGlobalState().getCommandLineOptions().outputDirectory ), + silo.setOutputDirectory( getOutputDirectory() ), + silo.setPlotLevel( m_plotLevel ); + silo.setWriteEdgeMesh( m_writeEdgeMesh ); + silo.setWriteFaceMesh( m_writeFaceMesh ); + silo.setWriteCellElementMesh( m_writeCellElementMesh ); + silo.setWriteFaceElementMesh( m_writeFaceElementMesh ); + silo.setOnlyPlotSpecifiedFieldNamesFlag( m_onlyPlotSpecifiedFieldNames ); + silo.setFieldNames( m_fieldNames.toViewConst() ); + silo.setPlotFileRoot( m_plotFileRoot ); + silo.initialize( numFiles ); + silo.waitForBatonWrite( rank, cycleNumber, eventCounter, false ); + silo.writeDomainPartition( domain, cycleNumber, time_n + dt * eventProgress, 0 ); + silo.handOffBaton(); + silo.clearEmptiesFromMultiObjects( cycleNumber ); + silo.finish(); + } return false; } diff --git a/src/coreComponents/fileIO/Outputs/SiloOutput.hpp b/src/coreComponents/fileIO/Outputs/SiloOutput.hpp index bb3574f3745..6b0a371974f 100644 --- a/src/coreComponents/fileIO/Outputs/SiloOutput.hpp +++ b/src/coreComponents/fileIO/Outputs/SiloOutput.hpp @@ -85,6 +85,12 @@ class SiloOutput : public OutputBase } siloOutputViewKeys; /// @endcond +protected: + /** + * @copydoc OutputBase::getTimerCategory + */ + logInfo::OutputTimerBase const & getTimerCategory() const override; + private: void postInputInitialization() override; diff --git a/src/coreComponents/fileIO/Outputs/TimeHistoryOutput.cpp b/src/coreComponents/fileIO/Outputs/TimeHistoryOutput.cpp index d01edac3dcf..7e1299aac57 100644 --- a/src/coreComponents/fileIO/Outputs/TimeHistoryOutput.cpp +++ b/src/coreComponents/fileIO/Outputs/TimeHistoryOutput.cpp @@ -25,6 +25,20 @@ namespace geos { using namespace dataRepository; +namespace logInfo +{ +struct TimeHistoryOutputTimer : public OutputTimerBase +{ + std::string_view getDescription() const override { return "Time history output timing"; } +}; +} + +logInfo::OutputTimerBase const & TimeHistoryOutput::getTimerCategory() const +{ + static logInfo::TimeHistoryOutputTimer timer; + return timer; +} + TimeHistoryOutput::TimeHistoryOutput( string const & name, Group * const parent ): OutputBase( name, parent ), @@ -166,12 +180,18 @@ bool TimeHistoryOutput::execute( real64 const GEOS_UNUSED_PARAM( time_n ), DomainPartition & GEOS_UNUSED_PARAM( domain ) ) { GEOS_MARK_FUNCTION; - localIndex newBuffered = m_io.front()->getBufferedCount( ); - for( auto & th_io : m_io ) + { - th_io->write( ); + Timer timer( m_outputTimer ); + + localIndex newBuffered = m_io.front()->getBufferedCount( ); + for( auto & th_io : m_io ) + { + th_io->write( ); + } + m_recordCount += newBuffered; } - m_recordCount += newBuffered; + return false; } diff --git a/src/coreComponents/fileIO/Outputs/TimeHistoryOutput.hpp b/src/coreComponents/fileIO/Outputs/TimeHistoryOutput.hpp index d7793943d9c..e2846f7c312 100644 --- a/src/coreComponents/fileIO/Outputs/TimeHistoryOutput.hpp +++ b/src/coreComponents/fileIO/Outputs/TimeHistoryOutput.hpp @@ -114,6 +114,12 @@ class TimeHistoryOutput : public OutputBase virtual PyTypeObject * getPythonType() const override; #endif +protected: + /** + * @copydoc OutputBase::getTimerCategory + */ + logInfo::OutputTimerBase const & getTimerCategory() const override; + private: /** diff --git a/src/coreComponents/fileIO/Outputs/VTKOutput.cpp b/src/coreComponents/fileIO/Outputs/VTKOutput.cpp index 15e4b1dac13..b6f4348899e 100644 --- a/src/coreComponents/fileIO/Outputs/VTKOutput.cpp +++ b/src/coreComponents/fileIO/Outputs/VTKOutput.cpp @@ -29,6 +29,20 @@ namespace geos using namespace dataRepository; +namespace logInfo +{ +struct VTKOutputTimer : public OutputTimerBase +{ + std::string_view getDescription() const override { return "VTK output timing"; } +}; +} + +logInfo::OutputTimerBase const & VTKOutput::getTimerCategory() const +{ + static logInfo::VTKOutputTimer timer; + return timer; +} + VTKOutput::VTKOutput( string const & name, Group * const parent ): OutputBase( name, parent ), @@ -80,7 +94,7 @@ VTKOutput::VTKOutput( string const & name, registerWrapper( viewKeysStruct::fieldNames, &m_fieldNames ). setRTTypeName( rtTypes::CustomTypes::groupNameRefArray ). setInputFlag( InputFlags::OPTIONAL ). - setDescription( "Names of the fields to output. If this attribute is specified, GEOSX outputs all the fields specified by the user, regardless of their `plotLevel`" ); + setDescription( "Names of the fields to output. If this attribute is specified, GEOS outputs all the fields specified by the user, regardless of their `plotLevel`" ); registerWrapper( viewKeysStruct::levelNames, &m_levelNames ). setInputFlag( InputFlags::OPTIONAL ). @@ -159,14 +173,20 @@ bool VTKOutput::execute( real64 const time_n, real64 const GEOS_UNUSED_PARAM ( eventProgress ), DomainPartition & domain ) { - GEOS_LOG_LEVEL_RANK_0( 1, GEOS_FMT( "{}: writing {} at time {} s (cycle number {})", getName(), m_fieldNames, time_n + dt, cycleNumber )); - - m_writer.setWriteGhostCells( m_writeGhostCells ); - m_writer.setWriteFaceElementsAs3D ( m_writeFaceElementsAs3D ); - m_writer.setOutputMode( m_writeBinaryData ); - m_writer.setOutputRegionType( m_outputRegionType ); - m_writer.setPlotLevel( m_plotLevel ); - m_writer.write( time_n, cycleNumber, domain ); + GEOS_MARK_FUNCTION; + + GEOS_LOG_LEVEL_RANK_0( 2, GEOS_FMT( "{}: writing {} at time {} s (cycle number {})", getName(), m_fieldNames, time_n + dt, cycleNumber )); + + { + Timer timer( m_outputTimer ); + + m_writer.setWriteGhostCells( m_writeGhostCells ); + m_writer.setWriteFaceElementsAs3D ( m_writeFaceElementsAs3D ); + m_writer.setOutputMode( m_writeBinaryData ); + m_writer.setOutputRegionType( m_outputRegionType ); + m_writer.setPlotLevel( m_plotLevel ); + m_writer.write( time_n, cycleNumber, domain ); + } return false; } diff --git a/src/coreComponents/fileIO/Outputs/VTKOutput.hpp b/src/coreComponents/fileIO/Outputs/VTKOutput.hpp index 5b167dc66d6..bf36e900eaf 100644 --- a/src/coreComponents/fileIO/Outputs/VTKOutput.hpp +++ b/src/coreComponents/fileIO/Outputs/VTKOutput.hpp @@ -77,6 +77,9 @@ class VTKOutput : public OutputBase DomainPartition & domain ) override { execute( time_n, 0, cycleNumber, eventCounter, eventProgress, domain ); + + // Call parent class cleanup to get the timing statistics + OutputBase::cleanup( time_n, cycleNumber, eventCounter, eventProgress, domain ); } /** @@ -109,6 +112,12 @@ class VTKOutput : public OutputBase virtual PyTypeObject * getPythonType() const override; #endif +protected: + /** + * @copydoc OutputBase::getTimerCategory + */ + logInfo::OutputTimerBase const & getTimerCategory() const override; + private: string m_plotFileRoot;