diff --git a/.integrated_tests.yaml b/.integrated_tests.yaml index 97f13dd906..0ce10cbe33 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 e5e337912c..8f34487782 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 9c7d0108c7..9eb10dbab9 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 278d8d6b45..1e65761384 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 2e99cfe6f4..8eec4acbef 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 a254aa9f97..07e5d54df9 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 b416573af0..e2f32348bc 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 8449d57e43..a2242e202f 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 1b4183c8c1..62df1c7408 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 d28ecf41b4..acd5c371d7 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 f36124d6ed..340d5c7f5f 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 ec9b3c5578..bce2f8cd26 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 c8ac5256d9..58ba51a9f9 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 14f702792d..138175d5c8 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 de9e92c8b8..71801e9ee5 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 17395f6232..8d7190cc70 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; }