fix: make new gravity treatment from #3337 an option

.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: ''

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".

src/coreComponents/physicsSolvers/fluidFlow/CompositionalMultiphaseBase.cpp

@@ -78,6 +78,7 @@
   m_allowCompDensChopping( 1 ),
   m_useTotalMassEquation( 1 ),
   m_useSimpleAccumulation( 1 ),
+  m_useNewGravity( 0 ),
   m_minCompDens( isothermalCompositionalMultiphaseBaseKernels::minDensForDivision )
 {
 //START_SPHINX_INCLUDE_00
@@ -164,6 +165,12 @@
     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 @@
       isothermalCompositionalMultiphaseBaseKernels::KernelLaunchSelector1
       < isothermalCompositionalMultiphaseFVMKernels::CFLFluxKernel >( numComps,
                                                                       numPhases,
+                                                                      m_useNewGravity,
                                                                       dt,
                                                                       stencilWrapper,
                                                                       compFlowAccessors.get( fields::flow::pressure{} ),

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;
 

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(),

src/coreComponents/physicsSolvers/fluidFlow/kernels/compositional/C1PPUPhaseFlux.hpp

@@ -79,6 +79,7 @@
   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 @@
     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 );

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 @@
 
 /******************************** 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 @@
     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 @@
   }
 }
 
-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 @@
   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,30 +188,32 @@
                                    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 );
   } );
 }
 
 #define INST_CFLFluxKernel( NC, STENCILWRAPPER_TYPE ) \
   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, \

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,

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],

src/coreComponents/physicsSolvers/fluidFlow/kernels/compositional/FluxComputeKernel.hpp

@@ -277,6 +277,7 @@
               ( m_numPhases,
               ip,
               m_kernelFlags.isSet( KernelFlags::CapPressure ),
+              m_kernelFlags.isSet( KernelFlags::NewGravity ),
               seri, sesri, sei,
               trans,
               dTrans_dPres,
@@ -303,6 +304,7 @@
               ( m_numPhases,
               ip,
               m_kernelFlags.isSet( KernelFlags::CapPressure ),
+              m_kernelFlags.isSet( KernelFlags::NewGravity ),
               seri, sesri, sei,
               trans,
               dTrans_dPres,
@@ -329,6 +331,7 @@
               ( m_numPhases,
               ip,
               m_kernelFlags.isSet( KernelFlags::CapPressure ),
+              m_kernelFlags.isSet( KernelFlags::NewGravity ),
               seri, sesri, sei,
               trans,
               dTrans_dPres,
@@ -352,7 +355,8 @@
 
           // 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 @@
                    string const & dofKey,
                    integer const hasCapPressure,
                    integer const useTotalMassEquation,
+                   integer const useNewGravity,
                    UpwindingParameters upwindingParams,
                    string const & solverName,
                    ElementRegionManager const & elemManager,
@@ -547,6 +552,8 @@
         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 );

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