Made density averaging optional (need gpu validation).

GEOS-DEV · Dec 10, 2024 · 8d0cd41 · 8d0cd41
1 parent 3ddd347
commit 8d0cd41
Show file tree

Hide file tree

Showing 3 changed files with 47 additions and 6 deletions.
diff --git a/src/coreComponents/physicsSolvers/fluidFlow/ImmiscibleMultiphaseFlow.cpp b/src/coreComponents/physicsSolvers/fluidFlow/ImmiscibleMultiphaseFlow.cpp
@@ -68,6 +68,12 @@ ImmiscibleMultiphaseFlow::ImmiscibleMultiphaseFlow( const string & name,
     setApplyDefaultValue( 1 ).
     setDescription( "Flag indicating whether total mass equation is used" );
 
+  this->registerWrapper( viewKeyStruct::gravityDensitySchemeString(), &m_gravityDensityScheme ).
+    setSizedFromParent( 0 ).
+    setInputFlag( InputFlags::OPTIONAL ).
+    setApplyDefaultValue( GravityDensityScheme::ArithmeticAverage ).
+    setDescription( "Scheme for density treatment in gravity" );
+
   this->registerWrapper( viewKeyStruct::solutionChangeScalingFactorString(), &m_solutionChangeScalingFactor ).
     setSizedFromParent( 0 ).
     setInputFlag( InputFlags::OPTIONAL ).
@@ -632,6 +638,7 @@ void ImmiscibleMultiphaseFlow::assembleFluxTerms( real64 const dt,
                                                                                    dofKey,
                                                                                    m_hasCapPressure,
                                                                                    m_useTotalMassEquation,
+                                                                                   m_gravityDensityScheme == GravityDensityScheme::PhasePresence,
                                                                                    getName(),
                                                                                    mesh.getElemManager(),
                                                                                    stencilWrapper,

diff --git a/src/coreComponents/physicsSolvers/fluidFlow/ImmiscibleMultiphaseFlow.hpp b/src/coreComponents/physicsSolvers/fluidFlow/ImmiscibleMultiphaseFlow.hpp
@@ -27,6 +27,27 @@
 namespace geos
 {
 
+/**
+ * @brief Options for density treatment in gravity
+ */
+enum class GravityDensityScheme : integer
+{
+  ArithmeticAverage, ///< average phase density is computed using simple arithmetic average:
+                     ///  rho_ave = 0.5 * (rho_i + rho_j)
+  PhasePresence,     ///< average phase density is computed using checking for phase presence:
+                     ///  rho_ave = 0.5 * (rho_i + rho_j) if phase is present in both cells i and j
+                     ///          = rho_i if phase is present in only cell i
+                     ///          = rho_j if phase is present in only cell j
+};
+
+/**
+ * @brief Strings for options for density treatment in gravity
+ */
+ENUM_STRINGS( GravityDensityScheme,
+              "ArithmeticAverage",
+              "PhasePresence" );
+
+
 //START_SPHINX_INCLUDE_00
 /**
  * @class ImmiscibleMultiphaseFlow
@@ -222,6 +243,9 @@ class ImmiscibleMultiphaseFlow : public FlowSolverBase
     static constexpr char const * relPermNamesString() { return "relPermNames"; }
     static constexpr char const * elemDofFieldString() { return "elemDofField"; }
 
+    // density averaging scheme
+    static constexpr char const * gravityDensitySchemeString()    { return "gravityDensityScheme"; }
+
     // time stepping controls
     static constexpr char const * solutionChangeScalingFactorString() { return "solutionChangeScalingFactor"; }
     static constexpr char const * targetRelativePresChangeString() { return "targetRelativePressureChangeInTimeStep"; }
@@ -287,6 +311,9 @@ class ImmiscibleMultiphaseFlow : public FlowSolverBase
   /// flag to determine whether or not to use total velocity formulation
   integer m_useTotalMassEquation;
 
+  /// scheme for density treatment in gravity
+  GravityDensityScheme m_gravityDensityScheme;
+
   /// target (relative) change in pressure in a time step
   real64 m_targetRelativePresChange;
 

diff --git a/...nts/physicsSolvers/fluidFlow/kernels/immiscibleMultiphase/ImmiscibleMultiphaseKernels.hpp b/...nts/physicsSolvers/fluidFlow/kernels/immiscibleMultiphase/ImmiscibleMultiphaseKernels.hpp
@@ -125,7 +125,8 @@ class FaceBasedAssemblyKernelBase
                                CRSMatrixView< real64, globalIndex const > const & localMatrix,
                                arrayView1d< real64 > const & localRhs,
                                integer const hasCapPressure,
-                               integer const useTotalMassEquation )
+                               integer const useTotalMassEquation,
+                               integer const checkPhasePresenceInGravity )
     : m_numPhases ( numPhases ),
     m_rankOffset( rankOffset ),
     m_dt( dt ),
@@ -145,7 +146,8 @@ class FaceBasedAssemblyKernelBase
     m_localMatrix( localMatrix ),
     m_localRhs( localRhs ),
     m_hasCapPressure ( hasCapPressure ),
-    m_useTotalMassEquation ( useTotalMassEquation )
+    m_useTotalMassEquation ( useTotalMassEquation ),
+    m_checkPhasePresenceInGravity ( checkPhasePresenceInGravity )
   {}
 
 protected:
@@ -197,6 +199,7 @@ class FaceBasedAssemblyKernelBase
   // Flags
   integer const m_hasCapPressure;
   integer const m_useTotalMassEquation;
+  integer const m_checkPhasePresenceInGravity;
 };
 
 /***************************************** */
@@ -255,7 +258,8 @@ class FaceBasedAssemblyKernel : public FaceBasedAssemblyKernelBase
                            CRSMatrixView< real64, globalIndex const > const & localMatrix,
                            arrayView1d< real64 > const & localRhs,
                            integer const hasCapPressure,
-                           integer const useTotalMassEquation )
+                           integer const useTotalMassEquation,
+                           integer const checkPhasePresenceInGravity )
     : FaceBasedAssemblyKernelBase( numPhases,
                                    rankOffset,
                                    dofNumberAccessor,
@@ -267,7 +271,8 @@ class FaceBasedAssemblyKernel : public FaceBasedAssemblyKernelBase
                                    localMatrix,
                                    localRhs,
                                    hasCapPressure,
-                                   useTotalMassEquation ),
+                                   useTotalMassEquation,
+                                   checkPhasePresenceInGravity ),
     m_stencilWrapper( stencilWrapper ),
     m_seri( stencilWrapper.getElementRegionIndices() ),
     m_sesri( stencilWrapper.getElementSubRegionIndices() ),
@@ -429,7 +434,7 @@ class FaceBasedAssemblyKernel : public FaceBasedAssemblyKernelBase
           {
             // density
             bool const phaseExists = (m_phaseVolFrac[seri[ke]][sesri[ke]][sei[ke]][ip] > 0);
-            if( !phaseExists )
+            if( m_checkPhasePresenceInGravity && !phaseExists )
             {
               continue;
             }
@@ -758,6 +763,7 @@ class FaceBasedAssemblyKernelFactory
                    string const & dofKey,
                    integer const hasCapPressure,
                    integer const useTotalMassEquation,
+                   integer const checkPhasePresenceInGravity,
                    string const & solverName,
                    ElementRegionManager const & elemManager,
                    STENCILWRAPPER const & stencilWrapper,
@@ -780,7 +786,8 @@ class FaceBasedAssemblyKernelFactory
 
     kernelType kernel( numPhases, rankOffset, stencilWrapper, dofNumberAccessor,
                        flowAccessors, fluidAccessors, capPressureAccessors, permAccessors,
-                       dt, localMatrix, localRhs, hasCapPressure, useTotalMassEquation );
+                       dt, localMatrix, localRhs, hasCapPressure, useTotalMassEquation,
+                       checkPhasePresenceInGravity );
     kernelType::template launch< POLICY >( stencilWrapper.size(), kernel );
   }
 };