From be0eb0e2e92d933f0f976a2c341d0ab4d9289ebf Mon Sep 17 00:00:00 2001 From: Pavel Tomin Date: Fri, 6 Dec 2024 21:03:14 -0600 Subject: [PATCH 1/3] 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/3] 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; From 7e2c33b57ab7617964c1467faaa8cfedbc5bd816 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Vidar=20Stiernstr=C3=B6m?= Date: Mon, 9 Dec 2024 16:23:06 -0800 Subject: [PATCH 3/3] feat: Rate-and-state friction with explicit time integration (#3450) * Add rate-and-state kernel for explicit time stepping and quasi-dynamic solver using Runge-Kutta 3(2) * Update strings for added rate-and-state fields * Change to using a single multi-d array for the Runge-Kutta stage rates * Add PID controller to adaptive time step update * Bracket slip rate after Newton update * Add kernel for EmbeddedRungeKutta time stepping of slip and state and refactor QuasiDynamicEQRK32 solver * Add some comments to QuasiDynamicEQRK32 solver * Add Bogacki-Shampine 3(2) Runge-Kutta method * Add PIDController class for time step control * Add xml input for explicit solver, update subrepo, add docs * Rename constructor arguments to prevent shadowing * rebaseline --------- Co-authored-by: Matteo Cusini <49037133+CusiniM@users.noreply.github.com> --- .integrated_tests.yaml | 2 +- BASELINE_NOTES.md | 4 + .../SpringSliderExplicit_base.xml | 167 ++++++ .../SpringSliderExplicit_smoke.xml | 32 ++ .../inducedSeismicity/SpringSlider_base.xml | 4 +- .../inducedSeismicity/inducedSeismicity.ats | 9 +- .../physicsSolvers/contact/ContactFields.hpp | 16 + .../inducedSeismicity/CMakeLists.txt | 4 +- .../inducedSeismicity/QuasiDynamicEQ.cpp | 8 +- .../inducedSeismicity/QuasiDynamicEQ.hpp | 2 +- .../inducedSeismicity/QuasiDynamicEQRK32.cpp | 478 ++++++++++++++++++ .../inducedSeismicity/QuasiDynamicEQRK32.hpp | 235 +++++++++ .../kernels/RateAndStateKernels.hpp | 409 ++++++++++++++- .../inducedSeismicity/rateAndStateFields.hpp | 54 +- src/coreComponents/schema/schema.xsd | 75 ++- src/coreComponents/schema/schema.xsd.other | 148 +++--- 16 files changed, 1539 insertions(+), 108 deletions(-) create mode 100644 inputFiles/inducedSeismicity/SpringSliderExplicit_base.xml create mode 100644 inputFiles/inducedSeismicity/SpringSliderExplicit_smoke.xml create mode 100644 src/coreComponents/physicsSolvers/inducedSeismicity/QuasiDynamicEQRK32.cpp create mode 100644 src/coreComponents/physicsSolvers/inducedSeismicity/QuasiDynamicEQRK32.hpp diff --git a/.integrated_tests.yaml b/.integrated_tests.yaml index 3aa31a7e4cb..31ce96119a9 100644 --- a/.integrated_tests.yaml +++ b/.integrated_tests.yaml @@ -1,6 +1,6 @@ baselines: bucket: geosx - baseline: integratedTests/baseline_integratedTests-pr3480-9217-caaecb8 + baseline: integratedTests/baseline_integratedTests-pr3450-9221-37d940c allow_fail: all: '' streak: '' diff --git a/BASELINE_NOTES.md b/BASELINE_NOTES.md index acbeb39e382..3b097223255 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 #3450 (2024-12-08) +===================== +Added test for explicit runge kutta sprinslider. + PR #3480 (2024-12-06) ===================== Add "logLevel" parameter under /Problem/Outputs in baseline files diff --git a/inputFiles/inducedSeismicity/SpringSliderExplicit_base.xml b/inputFiles/inducedSeismicity/SpringSliderExplicit_base.xml new file mode 100644 index 00000000000..a0c0381fb23 --- /dev/null +++ b/inputFiles/inducedSeismicity/SpringSliderExplicit_base.xml @@ -0,0 +1,167 @@ + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + \ No newline at end of file diff --git a/inputFiles/inducedSeismicity/SpringSliderExplicit_smoke.xml b/inputFiles/inducedSeismicity/SpringSliderExplicit_smoke.xml new file mode 100644 index 00000000000..2de09b79a03 --- /dev/null +++ b/inputFiles/inducedSeismicity/SpringSliderExplicit_smoke.xml @@ -0,0 +1,32 @@ + + + + + + + + + + + + + + + + + \ No newline at end of file diff --git a/inputFiles/inducedSeismicity/SpringSlider_base.xml b/inputFiles/inducedSeismicity/SpringSlider_base.xml index 9f141031607..5f6aefc94ad 100644 --- a/inputFiles/inducedSeismicity/SpringSlider_base.xml +++ b/inputFiles/inducedSeismicity/SpringSlider_base.xml @@ -123,7 +123,7 @@ initialCondition="1" objectPath="ElementRegions/Fault/FractureSubRegion" component="0" - scale="0." + scale="0.70710678118e-6" setNames="{all}"/> diff --git a/inputFiles/inducedSeismicity/inducedSeismicity.ats b/inputFiles/inducedSeismicity/inducedSeismicity.ats index 08b031177f0..813bdf08b8b 100644 --- a/inputFiles/inducedSeismicity/inducedSeismicity.ats +++ b/inputFiles/inducedSeismicity/inducedSeismicity.ats @@ -28,6 +28,13 @@ decks = [ partitions=((1, 1, 1), ), restart_step=0, check_step=3262, - restartcheck_params=RestartcheckParameters(atol=1e-4, rtol=1e-3)) + restartcheck_params=RestartcheckParameters(atol=1e-4, rtol=1e-3)), + TestDeck( + name="SpringSliderExplicit_smoke", + description="Spring slider 0D system", + partitions=((1, 1, 1), ), + restart_step=0, + check_step=532, + restartcheck_params=RestartcheckParameters(atol=1e-4, rtol=1e-3)) ] generate_geos_tests(decks) diff --git a/src/coreComponents/physicsSolvers/contact/ContactFields.hpp b/src/coreComponents/physicsSolvers/contact/ContactFields.hpp index 59671a5ad6b..3a5005774e9 100644 --- a/src/coreComponents/physicsSolvers/contact/ContactFields.hpp +++ b/src/coreComponents/physicsSolvers/contact/ContactFields.hpp @@ -74,6 +74,14 @@ DECLARE_FIELD( dispJump, WRITE_AND_READ, "Displacement jump vector in the local reference system" ); +DECLARE_FIELD( dispJump_n, + "displacementJump", + array2d< real64 >, + 0, + NOPLOT, + WRITE_AND_READ, + "Displacement jump vector in the local reference system at the current time-step" ); + DECLARE_FIELD( slip, "slip", array1d< real64 >, @@ -114,6 +122,14 @@ DECLARE_FIELD( traction, WRITE_AND_READ, "Fracture traction vector in the local reference system." ); +DECLARE_FIELD( traction_n, + "traction_n", + array2d< real64 >, + 0, + NOPLOT, + WRITE_AND_READ, + "Initial fracture traction vector in the local reference system at this time-step." ); + DECLARE_FIELD( deltaTraction, "deltaTraction", array2d< real64 >, diff --git a/src/coreComponents/physicsSolvers/inducedSeismicity/CMakeLists.txt b/src/coreComponents/physicsSolvers/inducedSeismicity/CMakeLists.txt index 7edf040454d..c8f23c55323 100644 --- a/src/coreComponents/physicsSolvers/inducedSeismicity/CMakeLists.txt +++ b/src/coreComponents/physicsSolvers/inducedSeismicity/CMakeLists.txt @@ -4,6 +4,7 @@ set( physicsSolvers_headers inducedSeismicity/inducedSeismicityFields.hpp inducedSeismicity/rateAndStateFields.hpp inducedSeismicity/QuasiDynamicEQ.hpp + inducedSeismicity/QuasiDynamicEQRK32.hpp inducedSeismicity/SeismicityRate.hpp inducedSeismicity/kernels/RateAndStateKernels.hpp inducedSeismicity/kernels/SeismicityRateKernels.hpp @@ -12,6 +13,7 @@ set( physicsSolvers_headers # Specify solver sources set( physicsSolvers_sources ${physicsSolvers_sources} - inducedSeismicity/QuasiDynamicEQ.cpp + inducedSeismicity/QuasiDynamicEQ.cpp + inducedSeismicity/QuasiDynamicEQRK32.cpp inducedSeismicity/SeismicityRate.cpp PARENT_SCOPE ) diff --git a/src/coreComponents/physicsSolvers/inducedSeismicity/QuasiDynamicEQ.cpp b/src/coreComponents/physicsSolvers/inducedSeismicity/QuasiDynamicEQ.cpp index 059fa9d33e1..d00e64357e0 100644 --- a/src/coreComponents/physicsSolvers/inducedSeismicity/QuasiDynamicEQ.cpp +++ b/src/coreComponents/physicsSolvers/inducedSeismicity/QuasiDynamicEQ.cpp @@ -32,6 +32,7 @@ namespace geos using namespace dataRepository; using namespace fields; using namespace constitutive; +using namespace rateAndStateKernels; QuasiDynamicEQ::QuasiDynamicEQ( const string & name, Group * const parent ): @@ -149,8 +150,8 @@ real64 QuasiDynamicEQ::solverStep( real64 const & time_n, /// 2. Solve for slip rate and state variable and, compute slip GEOS_LOG_LEVEL_RANK_0( 1, "Rate and State solver" ); - - integer const maxNewtonIter = m_nonlinearSolverParameters.m_maxIterNewton; + integer const maxIterNewton = m_nonlinearSolverParameters.m_maxIterNewton; + real64 const newtonTol = m_nonlinearSolverParameters.m_newtonTol; forDiscretizationOnMeshTargets( domain.getMeshBodies(), [&]( string const &, MeshLevel & mesh, arrayView1d< string const > const & regionNames ) @@ -161,7 +162,8 @@ real64 QuasiDynamicEQ::solverStep( real64 const & time_n, SurfaceElementSubRegion & subRegion ) { // solve rate and state equations. - rateAndStateKernels::createAndLaunch< parallelDevicePolicy<> >( subRegion, viewKeyStruct::frictionLawNameString(), m_shearImpedance, maxNewtonIter, time_n, dtStress ); + createAndLaunch< ImplicitFixedStressRateAndStateKernel, parallelDevicePolicy<> >( subRegion, viewKeyStruct::frictionLawNameString(), m_shearImpedance, maxIterNewton, newtonTol, time_n, + dtStress ); // save old state saveOldStateAndUpdateSlip( subRegion, dtStress ); } ); diff --git a/src/coreComponents/physicsSolvers/inducedSeismicity/QuasiDynamicEQ.hpp b/src/coreComponents/physicsSolvers/inducedSeismicity/QuasiDynamicEQ.hpp index edff334c003..97b202ee7e6 100644 --- a/src/coreComponents/physicsSolvers/inducedSeismicity/QuasiDynamicEQ.hpp +++ b/src/coreComponents/physicsSolvers/inducedSeismicity/QuasiDynamicEQ.hpp @@ -47,7 +47,7 @@ class QuasiDynamicEQ : public PhysicsSolverBase struct viewKeyStruct : public PhysicsSolverBase::viewKeyStruct { /// stress solver name - static constexpr char const * stressSolverNameString() { return "stressSolverName"; } + constexpr static char const * stressSolverNameString() { return "stressSolverName"; } /// Friction law name string constexpr static char const * frictionLawNameString() { return "frictionLawName"; } /// Friction law name string diff --git a/src/coreComponents/physicsSolvers/inducedSeismicity/QuasiDynamicEQRK32.cpp b/src/coreComponents/physicsSolvers/inducedSeismicity/QuasiDynamicEQRK32.cpp new file mode 100644 index 00000000000..7b569118934 --- /dev/null +++ b/src/coreComponents/physicsSolvers/inducedSeismicity/QuasiDynamicEQRK32.cpp @@ -0,0 +1,478 @@ +/* + * ------------------------------------------------------------------------------------------------------------ + * SPDX-License-Identifier: LGPL-2.1-only + * + * Copyright (c) 2016-2024 Lawrence Livermore National Security LLC + * Copyright (c) 2018-2024 Total, S.A + * Copyright (c) 2018-2024 The Board of Trustees of the Leland Stanford Junior University + * Copyright (c) 2023-2024 Chevron + * Copyright (c) 2019- GEOS/GEOSX Contributors + * All rights reserved + * + * See top level LICENSE, COPYRIGHT, CONTRIBUTORS, NOTICE, and ACKNOWLEDGEMENTS files for details. + * ------------------------------------------------------------------------------------------------------------ + */ + +/** + * @file QuasiDynamicEQRK32.cpp + */ + +#include "QuasiDynamicEQRK32.hpp" + +#include "dataRepository/InputFlags.hpp" +#include "mesh/DomainPartition.hpp" +#include "kernels/RateAndStateKernels.hpp" +#include "rateAndStateFields.hpp" +#include "physicsSolvers/contact/ContactFields.hpp" +#include "fieldSpecification/FieldSpecificationManager.hpp" + +namespace geos +{ + +using namespace dataRepository; +using namespace fields; +using namespace constitutive; +using namespace rateAndStateKernels; + +QuasiDynamicEQRK32::QuasiDynamicEQRK32( const string & name, + Group * const parent ): + PhysicsSolverBase( name, parent ), + m_stressSolver( nullptr ), + m_stressSolverName( "SpringSlider" ), + m_shearImpedance( 0.0 ), + m_butcherTable( BogackiShampine32Table()), // TODO: The butcher table should be specified in the XML input. + m_successfulStep( false ), + m_controller( PIDController( { 1.0/18.0, 1.0/9.0, 1.0/18.0 }, + 1.0e-6, 1.0e-6, 0.81 )) // TODO: The control parameters should be specified in the XML input +{ + this->registerWrapper( viewKeyStruct::shearImpedanceString(), &m_shearImpedance ). + setInputFlag( InputFlags::REQUIRED ). + setDescription( "Shear impedance." ); + + this->registerWrapper( viewKeyStruct::stressSolverNameString(), &m_stressSolverName ). + setInputFlag( InputFlags::OPTIONAL ). + setDescription( "Name of solver for computing stress. If empty, the spring-slider model is run." ); +} + +void QuasiDynamicEQRK32::postInputInitialization() +{ + + // Initialize member stress solver as specified in XML input + if( !m_stressSolverName.empty() ) + { + m_stressSolver = &this->getParent().getGroup< PhysicsSolverBase >( m_stressSolverName ); + } + + PhysicsSolverBase::postInputInitialization(); +} + +QuasiDynamicEQRK32::~QuasiDynamicEQRK32() +{ + // TODO Auto-generated destructor stub +} + + +void QuasiDynamicEQRK32::registerDataOnMesh( Group & meshBodies ) +{ + PhysicsSolverBase::registerDataOnMesh( meshBodies ); + + forDiscretizationOnMeshTargets( meshBodies, [&] ( string const &, + MeshLevel & mesh, + arrayView1d< string const > const & regionNames ) + { + ElementRegionManager & elemManager = mesh.getElemManager(); + + elemManager.forElementSubRegions< SurfaceElementSubRegion >( regionNames, + [&]( localIndex const, + SurfaceElementSubRegion & subRegion ) + { + // Scalar functions on fault + subRegion.registerField< rateAndState::stateVariable >( getName() ); + subRegion.registerField< rateAndState::stateVariable_n >( getName() ); + subRegion.registerField< rateAndState::slipRate >( getName() ); + + // Tangent (2-component) functions on fault + string const labels2Comp[2] = {"tangent1", "tangent2" }; + subRegion.registerField< rateAndState::slipVelocity >( getName() ). + setDimLabels( 1, labels2Comp ).reference().resizeDimension< 1 >( 2 ); + subRegion.registerField< rateAndState::slipVelocity_n >( getName() ). + setDimLabels( 1, labels2Comp ).reference().resizeDimension< 1 >( 2 ); + subRegion.registerField< rateAndState::deltaSlip >( getName() ). + setDimLabels( 1, labels2Comp ).reference().resizeDimension< 1 >( 2 ); + subRegion.registerField< rateAndState::deltaSlip_n >( getName() ). + setDimLabels( 1, labels2Comp ).reference().resizeDimension< 1 >( 2 ); + + // Runge-Kutta stage rates and error + integer const numRKComponents = 3; + subRegion.registerField< rateAndState::rungeKuttaStageRates >( getName() ).reference().resizeDimension< 1, 2 >( m_butcherTable.numStages, numRKComponents ); + subRegion.registerField< rateAndState::error >( getName() ).reference().resizeDimension< 1 >( numRKComponents ); + + + if( !subRegion.hasWrapper( contact::dispJump::key() )) + { + // 3-component functions on fault + string const labels3Comp[3] = { "normal", "tangent1", "tangent2" }; + subRegion.registerField< contact::dispJump >( getName() ). + setDimLabels( 1, labels3Comp ). + reference().resizeDimension< 1 >( 3 ); + subRegion.registerField< contact::dispJump_n >( getName() ). + setDimLabels( 1, labels3Comp ). + reference().resizeDimension< 1 >( 3 ); + subRegion.registerField< contact::traction >( getName() ). + setDimLabels( 1, labels3Comp ). + reference().resizeDimension< 1 >( 3 ); + subRegion.registerField< contact::traction_n >( getName() ). + setDimLabels( 1, labels3Comp ). + reference().resizeDimension< 1 >( 3 ); + + subRegion.registerWrapper< string >( viewKeyStruct::frictionLawNameString() ). + setPlotLevel( PlotLevel::NOPLOT ). + setRestartFlags( RestartFlags::NO_WRITE ). + setSizedFromParent( 0 ); + + string & frictionLawName = subRegion.getReference< string >( viewKeyStruct::frictionLawNameString() ); + frictionLawName = PhysicsSolverBase::getConstitutiveName< FrictionBase >( subRegion ); + GEOS_ERROR_IF( frictionLawName.empty(), GEOS_FMT( "{}: FrictionBase model not found on subregion {}", + getDataContext(), subRegion.getDataContext() ) ); + } + } ); + } ); +} + +real64 QuasiDynamicEQRK32::solverStep( real64 const & time_n, + real64 const & dt, + int const cycleNumber, + DomainPartition & domain ) +{ + if( cycleNumber == 0 ) + { + /// Apply initial conditions to the Fault + FieldSpecificationManager & fieldSpecificationManager = FieldSpecificationManager::getInstance(); + + forDiscretizationOnMeshTargets( domain.getMeshBodies(), [&]( string const &, + MeshLevel & mesh, + arrayView1d< string const > const & ) + + { + fieldSpecificationManager.applyInitialConditions( mesh ); + + } ); + saveState( domain ); + } + + real64 dtAdaptive = dt; + + GEOS_LOG_LEVEL_RANK_0( 1, "Begin adaptive time step" ); + while( true ) // Adaptive time step loop. Performs a Runge-Kutta time stepping with error control on state and slip + { + real64 dtStress; GEOS_UNUSED_VAR( dtStress ); + + // Initial Runge-Kutta stage + stepRateStateODEInitialSubstage( dtAdaptive, domain ); + real64 dtStage = m_butcherTable.c[1]*dtAdaptive; + dtStress = updateStresses( time_n, dtStage, cycleNumber, domain ); + updateSlipVelocity( time_n, dtStage, domain ); + + // Remaining stages + for( integer stageIndex = 1; stageIndex < m_butcherTable.numStages-1; stageIndex++ ) + { + stepRateStateODESubstage( stageIndex, dtAdaptive, domain ); + dtStage = m_butcherTable.c[stageIndex+1]*dtAdaptive; + dtStress = updateStresses( time_n, dtStage, cycleNumber, domain ); + updateSlipVelocity( time_n, dtStage, domain ); + } + + stepRateStateODEAndComputeError( dtAdaptive, domain ); + // Update timestep based on the time step error + real64 const dtNext = setNextDt( dtAdaptive, domain ); + if( m_successfulStep ) // set in setNextDt + { + // Compute stresses, and slip velocity and save results at updated time, + if( !m_butcherTable.FSAL ) + { + dtStress = updateStresses( time_n, dtAdaptive, cycleNumber, domain ); + updateSlipVelocity( time_n, dtAdaptive, domain ); + } + saveState( domain ); + // update the time step and exit the adaptive time step loop + dtAdaptive = dtNext; + break; + } + else + { + // Retry with updated time step + dtAdaptive = dtNext; + } + } + // return time step size achieved by stress solver + return dtAdaptive; +} + +void QuasiDynamicEQRK32::stepRateStateODEInitialSubstage( real64 const dt, DomainPartition & domain ) const +{ + + forDiscretizationOnMeshTargets( domain.getMeshBodies(), [&]( string const &, + MeshLevel & mesh, + arrayView1d< string const > const & regionNames ) + + { + mesh.getElemManager().forElementSubRegions< SurfaceElementSubRegion >( regionNames, + [&]( localIndex const, + SurfaceElementSubRegion & subRegion ) + { + + string const & fricitonLawName = subRegion.template getReference< string >( viewKeyStruct::frictionLawNameString() ); + RateAndStateFriction const & frictionLaw = getConstitutiveModel< RateAndStateFriction >( subRegion, fricitonLawName ); + rateAndStateKernels::EmbeddedRungeKuttaKernel rkKernel( subRegion, frictionLaw, m_butcherTable ); + arrayView3d< real64 > const rkStageRates = subRegion.getField< rateAndState::rungeKuttaStageRates >(); + + if( m_butcherTable.FSAL && m_successfulStep ) + { + forAll< parallelDevicePolicy<> >( subRegion.size(), [=] GEOS_HOST_DEVICE ( localIndex const k ) + { + rkKernel.updateStageRatesFSAL( k ); + rkKernel.updateStageValues( k, 1, dt ); + } ); + } + else + { + forAll< parallelDevicePolicy<> >( subRegion.size(), [=] GEOS_HOST_DEVICE ( localIndex const k ) + { + rkKernel.initialize( k ); + rkKernel.updateStageRates( k, 0 ); + rkKernel.updateStageValues( k, 1, dt ); + } ); + } + } ); + } ); +} + +void QuasiDynamicEQRK32::stepRateStateODESubstage( integer const stageIndex, + real64 const dt, + DomainPartition & domain ) const +{ + + forDiscretizationOnMeshTargets( domain.getMeshBodies(), [&]( string const &, + MeshLevel & mesh, + arrayView1d< string const > const & regionNames ) + + { + mesh.getElemManager().forElementSubRegions< SurfaceElementSubRegion >( regionNames, + [&]( localIndex const, + SurfaceElementSubRegion & subRegion ) + { + + string const & fricitonLawName = subRegion.template getReference< string >( viewKeyStruct::frictionLawNameString() ); + RateAndStateFriction const & frictionLaw = getConstitutiveModel< RateAndStateFriction >( subRegion, fricitonLawName ); + rateAndStateKernels::EmbeddedRungeKuttaKernel rkKernel( subRegion, frictionLaw, m_butcherTable ); + arrayView3d< real64 > const rkStageRates = subRegion.getField< rateAndState::rungeKuttaStageRates >(); + + forAll< parallelDevicePolicy<> >( subRegion.size(), [=] GEOS_HOST_DEVICE ( localIndex const k ) + { + rkKernel.updateStageRates( k, stageIndex ); + rkKernel.updateStageValues( k, stageIndex+1, dt ); + } ); + } ); + } ); +} + +void QuasiDynamicEQRK32::stepRateStateODEAndComputeError( real64 const dt, DomainPartition & domain ) const +{ + forDiscretizationOnMeshTargets( domain.getMeshBodies(), [&]( string const &, + MeshLevel & mesh, + arrayView1d< string const > const & regionNames ) + + { + mesh.getElemManager().forElementSubRegions< SurfaceElementSubRegion >( regionNames, + [&]( localIndex const, + SurfaceElementSubRegion & subRegion ) + { + + string const & fricitonLawName = subRegion.template getReference< string >( viewKeyStruct::frictionLawNameString() ); + RateAndStateFriction const & frictionLaw = getConstitutiveModel< RateAndStateFriction >( subRegion, fricitonLawName ); + rateAndStateKernels::EmbeddedRungeKuttaKernel rkKernel( subRegion, frictionLaw, m_butcherTable ); + arrayView3d< real64 > const rkStageRates = subRegion.getField< rateAndState::rungeKuttaStageRates >(); + if( m_butcherTable.FSAL ) + { + forAll< parallelDevicePolicy<> >( subRegion.size(), [=] GEOS_HOST_DEVICE ( localIndex const k ) + { + // Perform last stage rate update + rkKernel.updateStageRates( k, m_butcherTable.numStages-1 ); + // Update solution to final time and compute errors + rkKernel.updateSolutionAndLocalErrorFSAL( k, dt, m_controller.absTol, m_controller.relTol ); + } ); + } + else + { + forAll< parallelDevicePolicy<> >( subRegion.size(), [=] GEOS_HOST_DEVICE ( localIndex const k ) + { + // Perform last stage rate update + rkKernel.updateStageRates( k, m_butcherTable.numStages-1 ); + // Update solution to final time and compute errors + rkKernel.updateSolutionAndLocalError( k, dt, m_controller.absTol, m_controller.relTol ); + } ); + } + } ); + } ); +} + +real64 QuasiDynamicEQRK32::updateStresses( real64 const & time_n, + real64 const & dt, + const int cycleNumber, + DomainPartition & domain ) const +{ + GEOS_LOG_LEVEL_RANK_0( 1, "Stress solver" ); + // Call member variable stress solver to update the stress state + if( m_stressSolver ) + { + // 1. Solve the momentum balance + real64 const dtStress = m_stressSolver->solverStep( time_n, dt, cycleNumber, domain ); + + return dtStress; + } + else + { + // Spring-slider shear traction computation + forDiscretizationOnMeshTargets( domain.getMeshBodies(), [&]( string const &, + MeshLevel & mesh, + arrayView1d< string const > const & regionNames ) + + { + mesh.getElemManager().forElementSubRegions< SurfaceElementSubRegion >( regionNames, + [&]( localIndex const, + SurfaceElementSubRegion & subRegion ) + { + + arrayView2d< real64 const > const deltaSlip = subRegion.getField< rateAndState::deltaSlip >(); + arrayView2d< real64 > const traction = subRegion.getField< fields::contact::traction >(); + arrayView2d< real64 const > const traction_n = subRegion.getField< fields::contact::traction_n >(); + + string const & fricitonLawName = subRegion.template getReference< string >( viewKeyStruct::frictionLawNameString() ); + RateAndStateFriction const & frictionLaw = getConstitutiveModel< RateAndStateFriction >( subRegion, fricitonLawName ); + + RateAndStateFriction::KernelWrapper frictionKernelWrapper = frictionLaw.createKernelUpdates(); + + forAll< parallelDevicePolicy<> >( subRegion.size(), [=] GEOS_HOST_DEVICE ( localIndex const k ) + { + SpringSliderParameters springSliderParameters = SpringSliderParameters( traction[k][0], + frictionKernelWrapper.getACoefficient( k ), + frictionKernelWrapper.getBCoefficient( k ), + frictionKernelWrapper.getDcCoefficient( k ) ); + + + traction[k][1] = traction_n[k][1] + springSliderParameters.tauRate * dt + - springSliderParameters.springStiffness * deltaSlip[k][0]; + traction[k][2] = traction_n[k][2] + springSliderParameters.tauRate * dt + - springSliderParameters.springStiffness * deltaSlip[k][1]; + } ); + } ); + } ); + return dt; + } +} + +void QuasiDynamicEQRK32::updateSlipVelocity( real64 const & time_n, + real64 const & dt, + DomainPartition & domain ) const +{ + GEOS_LOG_LEVEL_RANK_0( 1, "Rate and State solver" ); + integer const maxIterNewton = m_nonlinearSolverParameters.m_maxIterNewton; + real64 const newtonTol = m_nonlinearSolverParameters.m_newtonTol; + forDiscretizationOnMeshTargets( domain.getMeshBodies(), [&]( string const &, + MeshLevel & mesh, + arrayView1d< string const > const & regionNames ) + + { + mesh.getElemManager().forElementSubRegions< SurfaceElementSubRegion >( regionNames, + [&]( localIndex const, + SurfaceElementSubRegion & subRegion ) + { + // solve rate and state equations. + rateAndStateKernels::createAndLaunch< rateAndStateKernels::ExplicitRateAndStateKernel, parallelDevicePolicy<> >( subRegion, viewKeyStruct::frictionLawNameString(), m_shearImpedance, + maxIterNewton, newtonTol, time_n, dt ); + } ); + } ); +} + +void QuasiDynamicEQRK32::saveState( DomainPartition & domain ) const +{ + forDiscretizationOnMeshTargets( domain.getMeshBodies(), [&]( string const &, + MeshLevel & mesh, + arrayView1d< string const > const & regionNames ) + + { + mesh.getElemManager().forElementSubRegions< SurfaceElementSubRegion >( regionNames, + [&]( localIndex const, + SurfaceElementSubRegion & subRegion ) + { + arrayView1d< real64 const > const stateVariable = subRegion.getField< rateAndState::stateVariable >(); + arrayView2d< real64 const > const slipVelocity = subRegion.getField< rateAndState::slipVelocity >(); + arrayView2d< real64 const > const deltaSlip = subRegion.getField< rateAndState::deltaSlip >(); + arrayView2d< real64 const > const dispJump = subRegion.getField< contact::dispJump >(); + arrayView2d< real64 const > const traction = subRegion.getField< contact::traction >(); + + arrayView1d< real64 > const stateVariable_n = subRegion.getField< rateAndState::stateVariable_n >(); + arrayView2d< real64 > const slipVelocity_n = subRegion.getField< rateAndState::slipVelocity_n >(); + arrayView2d< real64 > const deltaSlip_n = subRegion.getField< rateAndState::deltaSlip >(); + arrayView2d< real64 > const dispJump_n = subRegion.getField< contact::dispJump_n >(); + arrayView2d< real64 > const traction_n = subRegion.getField< contact::traction_n >(); + + forAll< parallelDevicePolicy<> >( subRegion.size(), [=] GEOS_HOST_DEVICE ( localIndex const k ) + { + stateVariable_n[k] = stateVariable[k]; + LvArray::tensorOps::copy< 2 >( deltaSlip_n[k], deltaSlip[k] ); + LvArray::tensorOps::copy< 2 >( slipVelocity_n[k], slipVelocity[k] ); + LvArray::tensorOps::copy< 3 >( dispJump_n[k], dispJump[k] ); + LvArray::tensorOps::copy< 3 >( traction_n[k], traction[k] ); + } ); + } ); + } ); +} + +real64 QuasiDynamicEQRK32::setNextDt( real64 const & currentDt, DomainPartition & domain ) +{ + + // Spring-slider shear traction computation + forDiscretizationOnMeshTargets( domain.getMeshBodies(), [&]( string const &, + MeshLevel const & mesh, + arrayView1d< string const > const & regionNames ) + + { + mesh.getElemManager().forElementSubRegions< SurfaceElementSubRegion >( regionNames, + [&]( localIndex const, + SurfaceElementSubRegion const & subRegion ) + { + arrayView2d< real64 const > const error = subRegion.getField< rateAndState::error >(); + + RAJA::ReduceSum< parallelDeviceReduce, real64 > scaledl2ErrorSquared( 0.0 ); + integer const N = subRegion.size(); + forAll< parallelDevicePolicy<> >( subRegion.size(), [=] GEOS_HOST_DEVICE ( localIndex const k ) + { + scaledl2ErrorSquared += LvArray::tensorOps::l2NormSquared< 3 >( error[k] ); + } ); + m_controller.errors[0] = LvArray::math::sqrt( MpiWrapper::sum( scaledl2ErrorSquared.get() / (3.0*N) )); + } ); + } ); + + // Compute update factor to currentDt using PID error controller + limiter + real64 const dtFactor = m_controller.computeUpdateFactor( m_butcherTable.algHighOrder, m_butcherTable.algLowOrder ); + real64 const nextDt = dtFactor*currentDt; + // Check if step was acceptable + m_successfulStep = (dtFactor >= m_controller.acceptSafety) ? true : false; + if( m_successfulStep ) + { + m_controller.errors[2] = m_controller.errors[1]; + m_controller.errors[1] = m_controller.errors[0]; + GEOS_LOG_LEVEL_RANK_0( 1, GEOS_FMT( "Adaptive time step successful. The next dt will be {:.2e} s", nextDt )); + } + else + { + GEOS_LOG_LEVEL_RANK_0( 1, GEOS_FMT( "Adaptive time step failed. The next dt will be {:.2e} s", nextDt )); + } + + return nextDt; +} + +REGISTER_CATALOG_ENTRY( PhysicsSolverBase, QuasiDynamicEQRK32, string const &, dataRepository::Group * const ) + +} // namespace geos diff --git a/src/coreComponents/physicsSolvers/inducedSeismicity/QuasiDynamicEQRK32.hpp b/src/coreComponents/physicsSolvers/inducedSeismicity/QuasiDynamicEQRK32.hpp new file mode 100644 index 00000000000..0979fd22dd7 --- /dev/null +++ b/src/coreComponents/physicsSolvers/inducedSeismicity/QuasiDynamicEQRK32.hpp @@ -0,0 +1,235 @@ +/* + * ------------------------------------------------------------------------------------------------------------ + * SPDX-License-Identifier: LGPL-2.1-only + * + * Copyright (c) 2016-2024 Lawrence Livermore National Security LLC + * Copyright (c) 2018-2024 Total, S.A + * Copyright (c) 2018-2024 The Board of Trustees of the Leland Stanford Junior University + * Copyright (c) 2023-2024 Chevron + * Copyright (c) 2019- GEOS/GEOSX Contributors + * All rights reserved + * + * See top level LICENSE, COPYRIGHT, CONTRIBUTORS, NOTICE, and ACKNOWLEDGEMENTS files for details. + * ------------------------------------------------------------------------------------------------------------ + */ + +#ifndef GEOS_PHYSICSSOLVERS_INDUCED_QUASIDYNAMICEQRK32_HPP +#define GEOS_PHYSICSSOLVERS_INDUCED_QUASIDYNAMICEQRK32_HPP + +#include "physicsSolvers/PhysicsSolverBase.hpp" +#include "kernels/RateAndStateKernels.hpp" + +namespace geos +{ + +class QuasiDynamicEQRK32 : public PhysicsSolverBase +{ +public: + /// The default nullary constructor is disabled to avoid compiler auto-generation: + QuasiDynamicEQRK32() = delete; + + /// The constructor needs a user-defined "name" and a parent Group (to place this instance in the tree structure of classes) + QuasiDynamicEQRK32( const string & name, + Group * const parent ); + + /// Destructor + virtual ~QuasiDynamicEQRK32() override; + + static string catalogName() { return "QuasiDynamicEQRK32"; } + + /** + * @return Get the final class Catalog name + */ + virtual string getCatalogName() const override { return catalogName(); } + + /// This method ties properties with their supporting mesh + virtual void registerDataOnMesh( Group & meshBodies ) override; + + struct viewKeyStruct : public PhysicsSolverBase::viewKeyStruct + { + /// stress solver name + constexpr static char const * stressSolverNameString() { return "stressSolverName"; } + /// Friction law name string + constexpr static char const * frictionLawNameString() { return "frictionLawName"; } + /// Friction law name string + constexpr static char const * shearImpedanceString() { return "shearImpedance"; } + /// target slip increment + constexpr static char const * timeStepTol() { return "timeStepTol"; } + }; + + virtual real64 solverStep( real64 const & time_n, + real64 const & dt, + integer const cycleNumber, + DomainPartition & domain ) override final; + + + virtual real64 setNextDt( real64 const & currentDt, + DomainPartition & domain ) override final; + + /** + * @brief Computes stage rates for the initial Runge-Kutta substage and updates slip and state + * @param dt + * @param domain + */ + void stepRateStateODEInitialSubstage( real64 const dt, DomainPartition & domain ) const; + + /** + * @brief Computes stage rates at the Runge-Kutta substage specified by stageIndex and updates slip and state + * @param stageIndex + * @param dt + * @param domain + */ + void stepRateStateODESubstage( integer const stageIndex, + real64 const dt, + DomainPartition & domain ) const; + + /** + * @brief Updates slip and state to t + dt and approximates the error + * @param dt + * @param domain + */ + void stepRateStateODEAndComputeError( real64 const dt, DomainPartition & domain ) const; + + real64 updateStresses( real64 const & time_n, + real64 const & dt, + const int cycleNumber, + DomainPartition & domain ) const; + + /** + * @brief Updates rate-and-state slip velocity + * @param domain + */ + void updateSlipVelocity( real64 const & time_n, + real64 const & dt, + DomainPartition & domain ) const; + + /** + * @brief save the current state + * @param domain + */ + void saveState( DomainPartition & domain ) const; + +private: + + virtual void postInputInitialization() override; + + + /// pointer to stress solver + PhysicsSolverBase * m_stressSolver; + + /// stress solver name + string m_stressSolverName; + + /// shear impedance + real64 m_shearImpedance; + + /// Runge-Kutta Butcher table (specifies the embedded RK method) + // TODO: The specific type should not be hardcoded! + // Should be possible to change RK-method based on the table. + rateAndStateKernels::BogackiShampine32Table m_butcherTable; + + bool m_successfulStep; // Flag indicating if the adative time step was accepted + + /** + * @brief Proportional-integral-derivative controller used for updating time step + * based error estimate in the current and previous time steps. + */ + class PIDController + { +public: + + GEOS_HOST_DEVICE + PIDController( std::array< const real64, 3 > const & cparams, + const real64 atol, + const real64 rtol, + const real64 safety ): + controlParameters{ cparams }, + absTol( atol ), + relTol( rtol ), + acceptSafety( safety ), + errors{ {0.0, 0.0, 0.0} } + {} + + /// Default copy constructor + PIDController( PIDController const & ) = default; + + /// Default move constructor + PIDController( PIDController && ) = default; + + /// Deleted default constructor + PIDController() = delete; + + /// Deleted copy assignment operator + PIDController & operator=( PIDController const & ) = delete; + + /// Deleted move assignment operator + PIDController & operator=( PIDController && ) = delete; + + /// Parameters for the PID error controller + const std::array< const real64, 3 > controlParameters; // Controller parameters + + real64 const absTol; // absolut tolerence + + real64 const relTol; // relative tolerence + + real64 const acceptSafety; // Acceptance safety + + std::array< real64, 3 > errors; // Errors for current and two previous updates + // stored as [n+1, n, n-1] + + real64 computeUpdateFactor( integer const algHighOrder, integer const algLowOrder ) + { + // PID error controller + limiter + real64 const k = LvArray::math::min( algHighOrder, algLowOrder ) + 1.0; + real64 const eps0 = 1.0/(errors[0] + std::numeric_limits< real64 >::epsilon()); // n + 1 + real64 const eps1 = 1.0/(errors[1] + std::numeric_limits< real64 >::epsilon()); // n + real64 const eps2 = 1.0/(errors[2] + std::numeric_limits< real64 >::epsilon()); // n-1 + // Compute update factor eps0^(beta0/k)*eps1^(beta1/k)*eps2^(beta2/k) where + // beta0 - beta2 are the control parameters. Also apply limiter to smoothen changes. + // Limiter is 1.0 + atan(x - 1.0). Here use atan(x) = atan2(x, 1.0). + return 1.0 + LvArray::math::atan2( pow( eps0, controlParameters[0] / k ) * + pow( eps1, controlParameters[1] / k ) * + pow( eps2, controlParameters[2] / k ) - 1.0, 1.0 ); + } + }; + + PIDController m_controller; + + + class SpringSliderParameters + { +public: + + GEOS_HOST_DEVICE + SpringSliderParameters( real64 const normalTraction, real64 const a, real64 const b, real64 const Dc ): + tauRate( 1e-4 ), + springStiffness( 0.0 ) + { + real64 const criticalStiffness = normalTraction * (b - a) / Dc; + springStiffness = 0.9 * criticalStiffness; + } + + /// Default copy constructor + SpringSliderParameters( SpringSliderParameters const & ) = default; + + /// Default move constructor + SpringSliderParameters( SpringSliderParameters && ) = default; + + /// Deleted default constructor + SpringSliderParameters() = delete; + + /// Deleted copy assignment operator + SpringSliderParameters & operator=( SpringSliderParameters const & ) = delete; + + /// Deleted move assignment operator + SpringSliderParameters & operator=( SpringSliderParameters && ) = delete; + + real64 tauRate; + + real64 springStiffness; + }; +}; + +} /* namespace geos */ + +#endif /* GEOS_PHYSICSSOLVERS_INDUCED_QUASIDYNAMICEQRK32_HPP */ diff --git a/src/coreComponents/physicsSolvers/inducedSeismicity/kernels/RateAndStateKernels.hpp b/src/coreComponents/physicsSolvers/inducedSeismicity/kernels/RateAndStateKernels.hpp index cf2fc306b6b..87d05f3b443 100644 --- a/src/coreComponents/physicsSolvers/inducedSeismicity/kernels/RateAndStateKernels.hpp +++ b/src/coreComponents/physicsSolvers/inducedSeismicity/kernels/RateAndStateKernels.hpp @@ -27,20 +27,38 @@ namespace geos namespace rateAndStateKernels { + +// TBD: Pass the kernel and add getters for relevant fields to make this function general purpose and avoid +// wrappers? +GEOS_HOST_DEVICE +static void projectSlipRateBase( localIndex const k, + real64 const frictionCoefficient, + real64 const shearImpedance, + arrayView2d< real64 const > const traction, + arrayView1d< real64 const > const slipRate, + arrayView2d< real64 > const slipVelocity ) +{ + // Project slip rate onto shear traction to get slip velocity components + real64 const frictionForce = traction[k][0] * frictionCoefficient; + real64 const projectionScaling = 1.0 / ( shearImpedance + frictionForce / slipRate[k] ); + slipVelocity[k][0] = projectionScaling * traction[k][1]; + slipVelocity[k][1] = projectionScaling * traction[k][2]; +} + /** - * @class RateAndStateKernel + * @class ImplicitFixedStressRateAndStateKernel * * @brief * * @details */ -class RateAndStateKernel +class ImplicitFixedStressRateAndStateKernel { public: - RateAndStateKernel( SurfaceElementSubRegion & subRegion, - constitutive::RateAndStateFriction const & frictionLaw, - real64 const shearImpedance ): + ImplicitFixedStressRateAndStateKernel( SurfaceElementSubRegion & subRegion, + constitutive::RateAndStateFriction const & frictionLaw, + real64 const shearImpedance ): m_slipRate( subRegion.getField< fields::rateAndState::slipRate >() ), m_stateVariable( subRegion.getField< fields::rateAndState::stateVariable >() ), m_stateVariable_n( subRegion.getField< fields::rateAndState::stateVariable_n >() ), @@ -90,7 +108,7 @@ class RateAndStateKernel stack.jacobian[0][0] = dFriction[0]; // derivative of Eq 1 w.r.t. stateVariable stack.jacobian[0][1] = dFriction[1]; // derivative of Eq 1 w.r.t. slipRate stack.jacobian[1][0] = dStateEvolutionLaw[0]; // derivative of Eq 2 w.r.t. stateVariable - stack.jacobian[1][1] = dStateEvolutionLaw[1]; // derivative of Eq 2 w.r.t. m_slipRate + stack.jacobian[1][1] = dStateEvolutionLaw[1]; // derivative of Eq 2 w.r.t. slipRate } GEOS_HOST_DEVICE @@ -109,11 +127,8 @@ class RateAndStateKernel GEOS_HOST_DEVICE void projectSlipRate( localIndex const k ) const { - // Project slip rate onto shear traction to get slip velocity components - real64 const frictionForce = m_traction[k][0] * m_frictionLaw.frictionCoefficient( k, m_slipRate[k], m_stateVariable[k] ); - real64 const projectionScaling = 1.0 / ( m_shearImpedance + frictionForce / m_slipRate[k] ); - m_slipVelocity[k][0] = projectionScaling * m_traction[k][1]; - m_slipVelocity[k][1] = projectionScaling * m_traction[k][2]; + real64 const frictionCoefficient = m_frictionLaw.frictionCoefficient( k, m_slipRate[k], m_stateVariable[k] ); + projectSlipRateBase( k, frictionCoefficient, m_shearImpedance, m_traction, m_slipRate, m_slipVelocity ); } GEOS_HOST_DEVICE @@ -144,18 +159,124 @@ class RateAndStateKernel }; +/** + * @class ExplicitRateAndStateKernel + * + * @brief + * + * @details + */ +class ExplicitRateAndStateKernel +{ +public: + + ExplicitRateAndStateKernel( SurfaceElementSubRegion & subRegion, + constitutive::RateAndStateFriction const & frictionLaw, + real64 const shearImpedance ): + m_slipRate( subRegion.getField< fields::rateAndState::slipRate >() ), + m_stateVariable( subRegion.getField< fields::rateAndState::stateVariable >() ), + m_traction( subRegion.getField< fields::contact::traction >() ), + m_slipVelocity( subRegion.getField< fields::rateAndState::slipVelocity >() ), + m_shearImpedance( shearImpedance ), + m_frictionLaw( frictionLaw.createKernelUpdates() ) + {} + + /** + * @struct StackVariables + * @brief Kernel variables located on the stack + */ + struct StackVariables + { +public: + + GEOS_HOST_DEVICE + StackVariables( ) + {} + + real64 jacobian; + real64 rhs; + + }; + + GEOS_HOST_DEVICE + void setup( localIndex const k, + real64 const dt, + StackVariables & stack ) const + { + GEOS_UNUSED_VAR( dt ); + real64 const normalTraction = m_traction[k][0]; + real64 const shearTractionMagnitude = LvArray::math::sqrt( m_traction[k][1] * m_traction[k][1] + m_traction[k][2] * m_traction[k][2] ); + + // Slip rate is bracketed between [0, shear traction magnitude / shear impedance] + // If slip rate is outside the bracket, re-initialize to the middle value + real64 const upperBound = shearTractionMagnitude/m_shearImpedance; + real64 const bracketedSlipRate = m_slipRate[k] > upperBound ? 0.5*upperBound : m_slipRate[k]; + + stack.rhs = shearTractionMagnitude - m_shearImpedance *bracketedSlipRate - normalTraction * m_frictionLaw.frictionCoefficient( k, bracketedSlipRate, m_stateVariable[k] ); + stack.jacobian = -m_shearImpedance - normalTraction * m_frictionLaw.dFrictionCoefficient_dSlipRate( k, bracketedSlipRate, m_stateVariable[k] ); + } + + GEOS_HOST_DEVICE + void solve( localIndex const k, + StackVariables & stack ) const + { + m_slipRate[k] -= stack.rhs/stack.jacobian; + + // Slip rate is bracketed between [0, shear traction magnitude / shear impedance] + // Check that the update did not end outside of the bracket. + real64 const shearTractionMagnitude = LvArray::math::sqrt( m_traction[k][1] * m_traction[k][1] + m_traction[k][2] * m_traction[k][2] ); + real64 const upperBound = shearTractionMagnitude/m_shearImpedance; + if( m_slipRate[k] > upperBound ) m_slipRate[k] = 0.5*upperBound; + + } + + + GEOS_HOST_DEVICE + camp::tuple< int, real64 > checkConvergence( StackVariables const & stack, + real64 const tol ) const + { + real64 const residualNorm = LvArray::math::abs( stack.rhs ); + int const converged = residualNorm < tol ? 1 : 0; + camp::tuple< int, real64 > result { converged, residualNorm }; + return result; + } + + GEOS_HOST_DEVICE + void projectSlipRate( localIndex const k ) const + { + real64 const frictionCoefficient = m_frictionLaw.frictionCoefficient( k, m_slipRate[k], m_stateVariable[k] ); + projectSlipRateBase( k, frictionCoefficient, m_shearImpedance, m_traction, m_slipRate, m_slipVelocity ); + } + +private: + + arrayView1d< real64 > const m_slipRate; + + arrayView1d< real64 > const m_stateVariable; + + arrayView2d< real64 const > const m_traction; + + arrayView2d< real64 > const m_slipVelocity; + + real64 const m_shearImpedance; + + constitutive::RateAndStateFriction::KernelWrapper m_frictionLaw; + +}; /** * @brief Performs the kernel launch + * @tparam KernelType The Rate-and-state kernel to launch * @tparam POLICY the policy used in the RAJA kernels */ -template< typename POLICY > +template< typename KernelType, typename POLICY > static void createAndLaunch( SurfaceElementSubRegion & subRegion, string const & frictionLawNameKey, real64 const shearImpedance, - integer const maxNewtonIter, + integer const maxIterNewton, + real64 const newtonTol, real64 const time_n, real64 const dt ) { @@ -165,20 +286,20 @@ createAndLaunch( SurfaceElementSubRegion & subRegion, string const & frictionaLawName = subRegion.getReference< string >( frictionLawNameKey ); constitutive::RateAndStateFriction const & frictionLaw = subRegion.getConstitutiveModel< constitutive::RateAndStateFriction >( frictionaLawName ); - RateAndStateKernel kernel( subRegion, frictionLaw, shearImpedance ); + KernelType kernel( subRegion, frictionLaw, shearImpedance ); // Newton loop (outside of the kernel launch) bool allConverged = false; - for( integer iter = 0; iter < maxNewtonIter; iter++ ) + for( integer iter = 0; iter < maxIterNewton; iter++ ) { RAJA::ReduceMin< parallelDeviceReduce, int > converged( 1 ); RAJA::ReduceMax< parallelDeviceReduce, real64 > residualNorm( 0.0 ); forAll< POLICY >( subRegion.size(), [=] GEOS_HOST_DEVICE ( localIndex const k ) { - RateAndStateKernel::StackVariables stack; + typename KernelType::StackVariables stack; kernel.setup( k, dt, stack ); kernel.solve( k, stack ); - auto const [elementConverged, elementResidualNorm] = kernel.checkConvergence( stack, 1.0e-6 ); + auto const [elementConverged, elementResidualNorm] = kernel.checkConvergence( stack, newtonTol ); converged.min( elementConverged ); residualNorm.max( elementResidualNorm ); } ); @@ -202,8 +323,260 @@ createAndLaunch( SurfaceElementSubRegion & subRegion, } ); } +/** + * @brief Butcher table for embedded RK3(2) method using Kuttas third order + * method for the high-order update, and an explicit trapezoidal rule + * based on the first and third stage rates for the low-order update. + */ +struct Kutta32Table +{ + integer constexpr static algHighOrder = 3; // High-order update order + integer constexpr static algLowOrder = 2; // Low-order update order + integer constexpr static numStages = 3; // Number of stages + real64 const a[2][2] = { { 1.0/2.0, 0.0 }, // Coefficients for stage value updates + { -1.0, 2.0 } }; // (lower-triangular part of table). + real64 const c[3] = { 0.0, 1.0/2.0, 1.0 }; // Coefficients for time increments of substages + real64 const b[3] = { 1.0/6.0, 4.0/6.0, 1.0/6.0 }; // Quadrature weights used to step the solution to next time + real64 const bStar[3] = { 1.0/2.0, 0.0, 1.0/2.0 }; // Quadrature weights used for low-order comparision solution + real64 constexpr static FSAL = false; // Not first same as last +}; + +/** + * @brief Butcher table for the BogackiShampine 3(2) method. + */ +struct BogackiShampine32Table +{ + integer constexpr static algHighOrder = 3; // High-order update order + integer constexpr static algLowOrder = 2; // Low-order update order + integer constexpr static numStages = 4; // Number of stages + real64 const a[3][3] = { { 1.0/2.0, 0.0, 0.0 }, // Coefficients for stage value updates + { 0.0, 3.0/4.0, 0.0 }, // (lower-triangular part of table). + { 2.0/9.0, 1.0/3.0, 4.0/9.0 } }; + real64 const c[4] = { 0.0, 1.0/2.0, 3.0/4.0, 1.0 }; // Coefficients for time increments of substages + real64 const b[4] = { 2.0/9.0, 1.0/3.0, 4.0/9.0, 0.0 }; // Quadrature weights used to step the solution to next time + real64 const bStar[4] = { 7.0/24.0, 1.0/4.0, 1.0/3.0, 1.0/8.0}; // Quadrature weights used for low-order comparision solution + bool constexpr static FSAL = true; // First same as last (can reuse the last stage rate in next + // update) +}; + +/** + * @brief Runge-Kutta method used to time integrate slip and state. Uses of a high order + * update used to integrate the solutions, and a lower order update to estimate the error + * in the time step. + * + * @tparam Butcher table defining the Runge-Kutta method. + */ +template< typename TABLE_TYPE > class EmbeddedRungeKuttaKernel +{ + +public: + EmbeddedRungeKuttaKernel( SurfaceElementSubRegion & subRegion, + constitutive::RateAndStateFriction const & frictionLaw, + TABLE_TYPE butcherTable ): + m_stateVariable( subRegion.getField< fields::rateAndState::stateVariable >() ), + m_stateVariable_n( subRegion.getField< fields::rateAndState::stateVariable_n >() ), + m_slipRate( subRegion.getField< fields::rateAndState::slipRate >() ), + m_slipVelocity( subRegion.getField< fields::rateAndState::slipVelocity >() ), + m_slipVelocity_n( subRegion.getField< fields::rateAndState::slipVelocity_n >() ), + m_deltaSlip( subRegion.getField< fields::rateAndState::deltaSlip >() ), + m_deltaSlip_n( subRegion.getField< fields::rateAndState::deltaSlip_n >() ), + m_dispJump( subRegion.getField< fields::contact::dispJump >() ), + m_dispJump_n( subRegion.getField< fields::contact::dispJump_n >() ), + m_error( subRegion.getField< fields::rateAndState::error >() ), + m_stageRates( subRegion.getField< fields::rateAndState::rungeKuttaStageRates >() ), + m_frictionLaw( frictionLaw.createKernelUpdates() ), + m_butcherTable( butcherTable ) + {} + + /** + * @brief Initialize slip and state buffers + */ + GEOS_HOST_DEVICE + void initialize( localIndex const k ) const + { + LvArray::tensorOps::copy< 2 >( m_slipVelocity[k], m_slipVelocity_n[k] ); + m_slipRate[k] = LvArray::tensorOps::l2Norm< 2 >( m_slipVelocity_n[k] ); + m_stateVariable[k] = m_stateVariable_n[k]; + } + + /** + * @brief Re-uses the last stage rate from the previous time step as the first + * in the next update. Only valid for FSAL (first-same-as-last) Runge-Kutta methods. + */ + GEOS_HOST_DEVICE + void updateStageRatesFSAL( localIndex const k ) const + { + LvArray::tensorOps::copy< 3 >( m_stageRates[k][0], m_stageRates[k][m_butcherTable.numStages-1] ); + } + + /** + * @brief Updates the stage rates rates (the right-hand-side of the ODEs for slip and state) + */ + GEOS_HOST_DEVICE + void updateStageRates( localIndex const k, integer const stageIndex ) const + { + m_stageRates[k][stageIndex][0] = m_slipVelocity[k][0]; + m_stageRates[k][stageIndex][1] = m_slipVelocity[k][1]; + m_stageRates[k][stageIndex][2] = m_frictionLaw.stateEvolution( k, m_slipRate[k], m_stateVariable[k] ); + } + + /** + * @brief Update stage values (slip, state and displacement jump) to a Runge-Kutta substage. + */ + GEOS_HOST_DEVICE + void updateStageValues( localIndex const k, integer const stageIndex, real64 const dt ) const + { + + real64 stateVariableIncrement = 0.0; + real64 deltaSlipIncrement[2] = {0.0, 0.0}; + + for( integer i = 0; i < stageIndex; i++ ) + { + deltaSlipIncrement[0] += m_butcherTable.a[stageIndex-1][i] * m_stageRates[k][i][0]; + deltaSlipIncrement[1] += m_butcherTable.a[stageIndex-1][i] * m_stageRates[k][i][1]; + stateVariableIncrement += m_butcherTable.a[stageIndex-1][i] * m_stageRates[k][i][2]; + } + m_deltaSlip[k][0] = m_deltaSlip_n[k][0] + dt*deltaSlipIncrement[0]; + m_deltaSlip[k][1] = m_deltaSlip_n[k][1] + dt*deltaSlipIncrement[1]; + m_stateVariable[k] = m_stateVariable_n[k] + dt*stateVariableIncrement; + + m_dispJump[k][1] = m_dispJump_n[k][1] + m_deltaSlip[k][0]; + m_dispJump[k][2] = m_dispJump_n[k][2] + m_deltaSlip[k][1]; + } + + /** + * @brief Updates slip, state and displacement jump to the next time computes error the local error + * in the time step + */ + GEOS_HOST_DEVICE + void updateSolutionAndLocalError( localIndex const k, real64 const dt, real64 const absTol, real64 const relTol ) const + { + + real64 deltaSlipIncrement[2] = {0.0, 0.0}; + real64 deltaSlipIncrementLowOrder[2] = {0.0, 0.0}; + + real64 stateVariableIncrement = 0.0; + real64 stateVariableIncrementLowOrder = 0.0; + + for( localIndex i = 0; i < m_butcherTable.numStages; i++ ) + { + + // High order update of solution + deltaSlipIncrement[0] += m_butcherTable.b[i] * m_stageRates[k][i][0]; + deltaSlipIncrement[1] += m_butcherTable.b[i] * m_stageRates[k][i][1]; + stateVariableIncrement += m_butcherTable.b[i] * m_stageRates[k][i][2]; + + // Low order update for error + deltaSlipIncrementLowOrder[0] += m_butcherTable.bStar[i] * m_stageRates[k][i][0]; + deltaSlipIncrementLowOrder[1] += m_butcherTable.bStar[i] * m_stageRates[k][i][1]; + stateVariableIncrementLowOrder += m_butcherTable.bStar[i] * m_stageRates[k][i][2]; + } + + m_deltaSlip[k][0] = m_deltaSlip_n[k][0] + dt * deltaSlipIncrement[0]; + m_deltaSlip[k][1] = m_deltaSlip_n[k][1] + dt * deltaSlipIncrement[1]; + m_stateVariable[k] = m_stateVariable_n[k] + dt * stateVariableIncrement; + + real64 const deltaSlipLowOrder[2] = {m_deltaSlip_n[k][0] + dt * deltaSlipIncrementLowOrder[0], + m_deltaSlip_n[k][1] + dt * deltaSlipIncrementLowOrder[1]}; + real64 const stateVariableLowOrder = m_stateVariable_n[k] + dt * stateVariableIncrementLowOrder; + + m_dispJump[k][1] = m_dispJump_n[k][1] + m_deltaSlip[k][0]; + m_dispJump[k][2] = m_dispJump_n[k][2] + m_deltaSlip[k][1]; + + // Compute error + m_error[k][0] = computeError( m_deltaSlip[k][0], deltaSlipLowOrder[0], absTol, relTol ); + m_error[k][1] = computeError( m_deltaSlip[k][1], deltaSlipLowOrder[1], absTol, relTol ); + m_error[k][2] = computeError( m_stateVariable[k], stateVariableLowOrder, absTol, relTol ); + } + + /** + * @brief Updates slip, state and displacement jump to the next time computes error the local error + * in the time step. Uses the FSAL (first-same-as-last) property. + */ + GEOS_HOST_DEVICE + void updateSolutionAndLocalErrorFSAL( localIndex const k, real64 const dt, real64 const absTol, real64 const relTol ) const + { + + real64 deltaSlipIncrementLowOrder[2] = {0.0, 0.0}; + real64 stateVariableIncrementLowOrder = 0.0; + + for( localIndex i = 0; i < m_butcherTable.numStages; i++ ) + { + // In FSAL algorithms the last RK substage update coincides with the + // high-order update. Only need to compute increments for the the + // low-order updates for error computation. + deltaSlipIncrementLowOrder[0] += m_butcherTable.bStar[i] * m_stageRates[k][i][0]; + deltaSlipIncrementLowOrder[1] += m_butcherTable.bStar[i] * m_stageRates[k][i][1]; + stateVariableIncrementLowOrder += m_butcherTable.bStar[i] * m_stageRates[k][i][2]; + } + + real64 const deltaSlipLowOrder[2] = {m_deltaSlip_n[k][0] + dt * deltaSlipIncrementLowOrder[0], + m_deltaSlip_n[k][1] + dt * deltaSlipIncrementLowOrder[1]}; + real64 const stateVariableLowOrder = m_stateVariable_n[k] + dt * stateVariableIncrementLowOrder; + + m_dispJump[k][1] = m_dispJump_n[k][1] + m_deltaSlip[k][0]; + m_dispJump[k][2] = m_dispJump_n[k][2] + m_deltaSlip[k][1]; + + // Compute error + m_error[k][0] = computeError( m_deltaSlip[k][0], deltaSlipLowOrder[0], absTol, relTol ); + m_error[k][1] = computeError( m_deltaSlip[k][1], deltaSlipLowOrder[1], absTol, relTol ); + m_error[k][2] = computeError( m_stateVariable[k], stateVariableLowOrder, absTol, relTol ); + } + + /** + * @brief Computes the relative error scaled by error tolerances + */ + GEOS_HOST_DEVICE + real64 computeError( real64 const highOrderApprox, real64 const lowOrderApprox, real64 const absTol, real64 const relTol ) const + { + return (highOrderApprox - lowOrderApprox) / + ( absTol + relTol * LvArray::math::max( LvArray::math::abs( highOrderApprox ), LvArray::math::abs( lowOrderApprox ) )); + } + +private: + + /// Current state variable + arrayView1d< real64 > const m_stateVariable; + + /// State variable at t = t_n + arrayView1d< real64 > const m_stateVariable_n; + + /// Current slip rate (magnitude of slip velocity) + arrayView1d< real64 > const m_slipRate; + + /// Current slip velocity + arrayView2d< real64 > const m_slipVelocity; + + /// Slip velocity at time t_n + arrayView2d< real64 > const m_slipVelocity_n; + + /// Current slip change + arrayView2d< real64 > const m_deltaSlip; + + /// Slip change at time t_n + arrayView2d< real64 > const m_deltaSlip_n; + + /// Current displacment jump + arrayView2d< real64 > const m_dispJump; + + /// Displacment jump at time t_n + arrayView2d< real64 > const m_dispJump_n; + + /// Local error for each solution component stored as slip1, slip2, state + arrayView2d< real64 > const m_error; + + /// Stage rates for each solution component stored as slip1, slip2, state + arrayView3d< real64 > const m_stageRates; + + /// Friction law used for rate-and-state updates + constitutive::RateAndStateFriction::KernelWrapper m_frictionLaw; + + /// Butcher table used for explicit time stepping of slip and state + TABLE_TYPE m_butcherTable; +}; + } /* namespace rateAndStateKernels */ -}/* namespace geos */ +} /* namespace geos */ #endif /* GEOS_PHYSICSSOLVERS_RATEANDSTATEKERNELS_HPP_ */ diff --git a/src/coreComponents/physicsSolvers/inducedSeismicity/rateAndStateFields.hpp b/src/coreComponents/physicsSolvers/inducedSeismicity/rateAndStateFields.hpp index 8031ab1c344..df5e63a2d19 100644 --- a/src/coreComponents/physicsSolvers/inducedSeismicity/rateAndStateFields.hpp +++ b/src/coreComponents/physicsSolvers/inducedSeismicity/rateAndStateFields.hpp @@ -36,10 +36,26 @@ DECLARE_FIELD( slipRate, "slipRate", array1d< real64 >, 1.0e-6, - LEVEL_0, + NOPLOT, WRITE_AND_READ, "Slip rate" ); +DECLARE_FIELD( slipVelocity, + "slipVelocity", + array2d< real64 >, + 0.70710678118e-6, + LEVEL_0, + WRITE_AND_READ, + "Slip velocity" ); + +DECLARE_FIELD( slipVelocity_n, + "slipVelocity_n", + array2d< real64 >, + 0.70710678118e-6, + NOPLOT, + WRITE_AND_READ, + "Slip velocity at previous time step" ); + DECLARE_FIELD( stateVariable, "stateVariable", array1d< real64 >, @@ -48,21 +64,14 @@ DECLARE_FIELD( stateVariable, WRITE_AND_READ, "Rate- and state-dependent friction state variable" ); -DECLARE_FIELD( slipVelocity, - "slipVelocity", - array2d< real64 >, - 1.0e-6, - LEVEL_0, - WRITE_AND_READ, - "Slip velocity" ); - DECLARE_FIELD( stateVariable_n, "stateVariable_n", array1d< real64 >, 0.6, NOPLOT, WRITE_AND_READ, - "Rate- and state-dependent friction state variable at previous time step" ); + "Initial rate- and state-dependent friction state variable at this time step" ); + DECLARE_FIELD( deltaSlip, "deltaSlip", @@ -72,6 +81,31 @@ DECLARE_FIELD( deltaSlip, WRITE_AND_READ, "Slip increment" ); +DECLARE_FIELD( deltaSlip_n, + "deltaSlip_n", + array2d< real64 >, + 0.0, + NOPLOT, + WRITE_AND_READ, + "Initial slip increment at this time step" ); + + +DECLARE_FIELD( rungeKuttaStageRates, + "rungeKuttaStageRates", + array3d< real64 >, + 0.0, + NOPLOT, + WRITE_AND_READ, + "Runge-Kutta stage rates for rate-and-state variables" ); + + +DECLARE_FIELD( error, + "error", + array2d< real64 >, + 0.0, + LEVEL_0, + WRITE_AND_READ, + "Error for rate-and-state fields" ); } diff --git a/src/coreComponents/schema/schema.xsd b/src/coreComponents/schema/schema.xsd index 7a6e7967c03..e130b1e56d4 100644 --- a/src/coreComponents/schema/schema.xsd +++ b/src/coreComponents/schema/schema.xsd @@ -399,6 +399,10 @@ + + + + @@ -2109,6 +2113,11 @@ the relative residual norm satisfies: + + @@ -2125,6 +2134,11 @@ the relative residual norm satisfies: + + @@ -2139,6 +2153,11 @@ the relative residual norm satisfies: + + @@ -2147,6 +2166,11 @@ the relative residual norm satisfies: + + @@ -2157,6 +2181,11 @@ the relative residual norm satisfies: + + @@ -2195,7 +2224,7 @@ the relative residual norm satisfies: - + @@ -2262,6 +2291,7 @@ the relative residual norm satisfies: + @@ -2647,6 +2677,8 @@ Level 0 outputs no specific information for this solver. Higher levels require m + + @@ -2734,6 +2766,8 @@ Level 0 outputs no specific information for this solver. Higher levels require m + + @@ -3211,7 +3245,7 @@ Level 0 outputs no specific information for this solver. Higher levels require m - + @@ -3234,7 +3268,7 @@ Local- Add jump stabilization on interior of macro elements--> - + @@ -3538,6 +3572,41 @@ Level 0 outputs no specific information for this solver. Higher levels require m + + + + + + + + + + + + + + + + + + + + + + + + diff --git a/src/coreComponents/schema/schema.xsd.other b/src/coreComponents/schema/schema.xsd.other index 33386f19956..3abf5b863f0 100644 --- a/src/coreComponents/schema/schema.xsd.other +++ b/src/coreComponents/schema/schema.xsd.other @@ -486,7 +486,7 @@ - + @@ -538,6 +538,7 @@ + @@ -571,7 +572,7 @@ - + @@ -608,7 +609,7 @@ - + @@ -659,7 +660,7 @@ - + @@ -700,7 +701,7 @@ - + @@ -733,7 +734,7 @@ - + @@ -744,7 +745,7 @@ - + @@ -757,7 +758,7 @@ - + @@ -770,7 +771,7 @@ - + @@ -786,7 +787,7 @@ - + @@ -820,7 +821,7 @@ - + @@ -883,7 +884,7 @@ - + @@ -914,7 +915,7 @@ - + @@ -927,7 +928,7 @@ - + @@ -940,7 +941,7 @@ - + @@ -953,7 +954,7 @@ - + @@ -966,7 +967,7 @@ - + @@ -981,7 +982,7 @@ - + @@ -992,7 +993,7 @@ - + @@ -1005,7 +1006,7 @@ - + @@ -1016,7 +1017,7 @@ - + @@ -1027,7 +1028,18 @@ - + + + + + + + + + + + + @@ -1038,7 +1050,7 @@ - + @@ -1051,7 +1063,7 @@ - + @@ -1062,7 +1074,7 @@ - + @@ -1073,7 +1085,7 @@ - + @@ -1086,7 +1098,7 @@ - + @@ -1101,7 +1113,7 @@ - + @@ -1116,7 +1128,7 @@ - + @@ -1129,7 +1141,7 @@ - + @@ -1144,7 +1156,7 @@ - + @@ -1155,7 +1167,7 @@ - + @@ -1168,7 +1180,7 @@ - + @@ -1181,7 +1193,7 @@ - + @@ -1196,7 +1208,7 @@ - + @@ -1212,7 +1224,7 @@ - + @@ -1227,7 +1239,7 @@ - + @@ -1244,7 +1256,7 @@ - + @@ -1261,7 +1273,7 @@ - + @@ -1278,7 +1290,7 @@ - + @@ -1293,7 +1305,7 @@ - + @@ -1306,7 +1318,7 @@ - + @@ -1345,7 +1357,7 @@ - + @@ -1374,7 +1386,7 @@ - + @@ -1467,7 +1479,7 @@ - + @@ -3091,7 +3103,7 @@ - + @@ -3119,7 +3131,7 @@ - + @@ -3138,11 +3150,11 @@ - + - + @@ -3152,7 +3164,7 @@ - + @@ -3162,11 +3174,11 @@ - + - + @@ -3176,7 +3188,7 @@ - + @@ -3186,7 +3198,7 @@ - + @@ -3196,7 +3208,7 @@ - + @@ -3220,7 +3232,7 @@ - + @@ -3238,7 +3250,7 @@ - + @@ -3250,7 +3262,7 @@ - + @@ -3262,7 +3274,7 @@ - + @@ -3270,11 +3282,11 @@ - + - + @@ -3297,7 +3309,7 @@ - + @@ -3323,7 +3335,7 @@ - + @@ -3344,7 +3356,7 @@ - + @@ -3374,7 +3386,7 @@ - + @@ -3388,7 +3400,7 @@ - + @@ -3415,7 +3427,7 @@ - + @@ -3454,7 +3466,7 @@ - +