fix: Add HU simplified 2-phase version, refactor IHU implementation and fix bugs where possible, change phase/component flux logic

src/coreComponents/finiteVolume/FluxApproximationBase.hpp

@@ -39,7 +39,8 @@ enum class UpwindingScheme : integer
 {
   PPU,    ///< PPU upwinding
   C1PPU,  ///< C1-PPU upwinding from https://doi.org/10.1016/j.advwatres.2017.07.028
-  IHU ///< IHU as in https://link.springer.com/content/pdf/10.1007/s10596-019-09835-6.pdf
+  IHU,    ///< IHU as in https://link.springer.com/content/pdf/10.1007/s10596-019-09835-6.pdf
+  HU2PH   ///< HU simplified 2-phase version
 };
 
 /**
@@ -48,7 +49,8 @@ enum class UpwindingScheme : integer
 ENUM_STRINGS( UpwindingScheme,
               "PPU",
               "C1PPU",
-              "IHU" );
+              "IHU",
+              "HU2PH" );
 
 /**
  * @struct UpwindingParameters
@@ -57,7 +59,7 @@ ENUM_STRINGS( UpwindingScheme,
  */
 struct UpwindingParameters
 {
-  /// PPU or C1-PPU or IHU
+  /// PPU or C1-PPU or IHU or HU2PH
   UpwindingScheme upwindingScheme;
 
   /// C1-PPU smoothing tolerance

src/coreComponents/physicsSolvers/fluidFlow/CMakeLists.txt

@@ -62,6 +62,7 @@ set( physicsSolvers_headers
      fluidFlow/kernels/compositional/GlobalComponentFractionKernel.hpp
      fluidFlow/kernels/compositional/HydrostaticPressureKernel.hpp
      fluidFlow/kernels/compositional/IHUPhaseFlux.hpp
+     fluidFlow/kernels/compositional/HU2PhaseFlux.hpp
      fluidFlow/kernels/compositional/KernelLaunchSelectors.hpp
      fluidFlow/kernels/compositional/PhaseComponentFlux.hpp
      fluidFlow/kernels/compositional/PhaseMobilityKernel.hpp

src/coreComponents/physicsSolvers/fluidFlow/CompositionalMultiphaseFVM.cpp

@@ -166,6 +166,14 @@ void CompositionalMultiphaseFVM::initializePreSubGroups()
                                                 : LinearSolverParameters::MGR::StrategyType::compositionalMultiphaseFVM;
 
   checkDiscretizationName();
+
+  DomainPartition & domain = this->getGroupByPath< DomainPartition >( "/Problem/domain" );
+  NumericalMethodsManager const & numericalMethodManager = domain.getNumericalMethodManager();
+  FiniteVolumeManager const & fvManager = numericalMethodManager.getFiniteVolumeManager();
+  FluxApproximationBase const & fluxApprox = fvManager.getFluxApproximation( m_discretizationName );
+  GEOS_ERROR_IF( fluxApprox.upwindingParams().upwindingScheme == UpwindingScheme::HU2PH && m_numPhases != 2,
+                 GEOS_FMT( "{}: upwinding scheme {} only supports 2-phase flow",
+                           getName(), EnumStrings< UpwindingScheme >::toString( UpwindingScheme::HU2PH )));
 }
 
 void CompositionalMultiphaseFVM::setupDofs( DomainPartition const & domain,

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

@@ -26,7 +26,6 @@
 #include "constitutive/capillaryPressure/layouts.hpp"
 #include "mesh/ElementRegionManager.hpp"
 #include "physicsSolvers/fluidFlow/kernels/compositional/PotGrad.hpp"
-#include "physicsSolvers/fluidFlow/kernels/compositional/PhaseComponentFlux.hpp"
 
 
 namespace geos
@@ -67,7 +66,6 @@ struct C1PPUPhaseFlux
    * @param dPhaseMassDens derivative of phase mass density wrt pressure, temperature, comp fraction
    * @param phaseCapPressure phase capillary pressure
    * @param dPhaseCapPressure_dPhaseVolFrac derivative of phase capillary pressure wrt phase volume fraction
-   * @param k_up uptream index for this phase
    * @param potGrad potential gradient for this phase
    * @param phaseFlux phase flux
    * @param dPhaseFlux_dP derivative of phase flux wrt pressure
@@ -91,21 +89,15 @@ struct C1PPUPhaseFlux
            ElementViewConst< arrayView3d< real64 const, compflow::USD_PHASE_DC > > const & dPhaseMob,
            ElementViewConst< arrayView2d< real64 const, compflow::USD_PHASE > > const & phaseVolFrac,
            ElementViewConst< arrayView3d< real64 const, compflow::USD_PHASE_DC > > const & dPhaseVolFrac,
-           ElementViewConst< arrayView4d< real64 const, constitutive::multifluid::USD_PHASE_COMP > > const & phaseCompFrac,
-           ElementViewConst< arrayView5d< real64 const, constitutive::multifluid::USD_PHASE_COMP_DC > > const & dPhaseCompFrac,
            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,
            ElementViewConst< arrayView3d< real64 const, constitutive::cappres::USD_CAPPRES > > const & phaseCapPressure,
            ElementViewConst< arrayView4d< real64 const, constitutive::cappres::USD_CAPPRES_DS > > const & dPhaseCapPressure_dPhaseVolFrac,
-           localIndex & k_up,
            real64 & potGrad,
            real64 ( &phaseFlux ),
            real64 ( & dPhaseFlux_dP )[numFluxSupportPoints],
-           real64 ( & dPhaseFlux_dC )[numFluxSupportPoints][numComp],
-           real64 ( & compFlux )[numComp],
-           real64 ( & dCompFlux_dP )[numFluxSupportPoints][numComp],
-           real64 ( & dCompFlux_dC )[numFluxSupportPoints][numComp][numComp] )
+           real64 ( & dPhaseFlux_dC )[numFluxSupportPoints][numComp] )
   {
     real64 dPotGrad_dTrans = 0.0;
     real64 dPresGrad_dP[numFluxSupportPoints]{};
@@ -214,13 +206,6 @@ struct C1PPUPhaseFlux
     }
 
     potGrad = potGrad * Ttrans;
-
-    // choose upstream cell for composition upwinding
-    k_up = (phaseFlux >= 0) ? 0 : 1;
-
-    //distribute on phaseComponentFlux here
-    PhaseComponentFlux::compute( ip, k_up, seri, sesri, sei, phaseCompFrac, dPhaseCompFrac, dCompFrac_dCompDens, phaseFlux
-                                 , dPhaseFlux_dP, dPhaseFlux_dC, compFlux, dCompFlux_dP, dCompFlux_dC );
   }
 };
 

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

@@ -35,6 +35,8 @@
 #include "physicsSolvers/fluidFlow/kernels/compositional/PPUPhaseFlux.hpp"
 #include "physicsSolvers/fluidFlow/kernels/compositional/C1PPUPhaseFlux.hpp"
 #include "physicsSolvers/fluidFlow/kernels/compositional/IHUPhaseFlux.hpp"
+#include "physicsSolvers/fluidFlow/kernels/compositional/HU2PhaseFlux.hpp"
+#include "physicsSolvers/fluidFlow/kernels/compositional/PhaseComponentFlux.hpp"
 
 namespace geos
 {
@@ -228,6 +230,7 @@ class FluxComputeKernel : public FluxComputeKernelBase
                     StackVariables & stack,
                     FUNC && compFluxKernelOp = NoOpFunc{} ) const
   {
+    using namespace isothermalCompositionalMultiphaseFVMKernelUtilities;
 
     // first, compute the transmissibilities at this face
     m_stencilWrapper.computeWeights( iconn,
@@ -269,12 +272,11 @@ class FluxComputeKernel : public FluxComputeKernelBase
           real64 phaseFlux = 0.0;
           real64 dPhaseFlux_dP[numFluxSupportPoints]{};
           real64 dPhaseFlux_dC[numFluxSupportPoints][numComp]{};
-
-          localIndex k_up = -1;
+          real64 dPhaseFlux_dTrans = 0.0; // not really used
 
           if( m_kernelFlags.isSet( KernelFlags::C1PPU ) )
           {
-            isothermalCompositionalMultiphaseFVMKernelUtilities::C1PPUPhaseFlux::compute< numComp, numFluxSupportPoints >
+            C1PPUPhaseFlux::compute< numComp, numFluxSupportPoints >
               ( m_numPhases,
               ip,
               m_kernelFlags.isSet( KernelFlags::CapPressure ),
@@ -286,22 +288,17 @@ class FluxComputeKernel : public FluxComputeKernelBase
               m_gravCoef,
               m_phaseMob, m_dPhaseMob,
               m_phaseVolFrac, m_dPhaseVolFrac,
-              m_phaseCompFrac, m_dPhaseCompFrac,
               m_dCompFrac_dCompDens,
               m_phaseMassDens, m_dPhaseMassDens,
               m_phaseCapPressure, m_dPhaseCapPressure_dPhaseVolFrac,
-              k_up,
               potGrad,
               phaseFlux,
               dPhaseFlux_dP,
-              dPhaseFlux_dC,
-              compFlux,
-              dCompFlux_dP,
-              dCompFlux_dC );
+              dPhaseFlux_dC );
           }
           else if( m_kernelFlags.isSet( KernelFlags::IHU ) )
           {
-            isothermalCompositionalMultiphaseFVMKernelUtilities::IHUPhaseFlux::compute< numComp, numFluxSupportPoints >
+            IHUPhaseFlux::compute< numComp, numFluxSupportPoints >
               ( m_numPhases,
               ip,
               m_kernelFlags.isSet( KernelFlags::CapPressure ),
@@ -313,22 +310,39 @@ class FluxComputeKernel : public FluxComputeKernelBase
               m_gravCoef,
               m_phaseMob, m_dPhaseMob,
               m_phaseVolFrac, m_dPhaseVolFrac,
-              m_phaseCompFrac, m_dPhaseCompFrac,
               m_dCompFrac_dCompDens,
               m_phaseMassDens, m_dPhaseMassDens,
               m_phaseCapPressure, m_dPhaseCapPressure_dPhaseVolFrac,
-              k_up,
               potGrad,
               phaseFlux,
               dPhaseFlux_dP,
-              dPhaseFlux_dC,
-              compFlux,
-              dCompFlux_dP,
-              dCompFlux_dC );
+              dPhaseFlux_dC );
+          }
+          else if( m_kernelFlags.isSet( KernelFlags::HU2PH ) )
+          {
+            HU2PhaseFlux::compute< numComp, numFluxSupportPoints >
+              ( m_numPhases,
+              ip,
+              m_kernelFlags.isSet( KernelFlags::CapPressure ),
+              m_kernelFlags.isSet( KernelFlags::CheckPhasePresenceInGravity ),
+              seri, sesri, sei,
+              trans,
+              dTrans_dPres,
+              m_pres,
+              m_gravCoef,
+              m_phaseMob, m_dPhaseMob,
+              m_phaseVolFrac, m_dPhaseVolFrac,
+              m_dCompFrac_dCompDens,
+              m_phaseMassDens, m_dPhaseMassDens,
+              m_phaseCapPressure, m_dPhaseCapPressure_dPhaseVolFrac,
+              potGrad,
+              phaseFlux,
+              dPhaseFlux_dP,
+              dPhaseFlux_dC );
           }
           else
           {
-            isothermalCompositionalMultiphaseFVMKernelUtilities::PPUPhaseFlux::compute< numComp, numFluxSupportPoints >
+            PPUPhaseFlux::compute< numComp, numFluxSupportPoints >
               ( m_numPhases,
               ip,
               m_kernelFlags.isSet( KernelFlags::CapPressure ),
@@ -340,21 +354,25 @@ class FluxComputeKernel : public FluxComputeKernelBase
               m_gravCoef,
               m_phaseMob, m_dPhaseMob,
               m_phaseVolFrac, m_dPhaseVolFrac,
-              m_phaseCompFrac, m_dPhaseCompFrac,
               m_dCompFrac_dCompDens,
               m_phaseMassDens, m_dPhaseMassDens,
               m_phaseCapPressure, m_dPhaseCapPressure_dPhaseVolFrac,
-              k_up,
               potGrad,
               phaseFlux,
               dPhaseFlux_dP,
               dPhaseFlux_dC,
-              compFlux,
-              dCompFlux_dP,
-              dCompFlux_dC,
-              dCompFlux_dTrans );
+              dPhaseFlux_dTrans );
           }
 
+          // choose upstream cell for composition upwinding
+          localIndex const k_up = (phaseFlux >= 0) ? 0 : 1;
+
+          // distribute on phaseComponentFlux here
+          PhaseComponentFlux::compute( ip, k_up, seri, sesri, sei,
+                                       m_phaseCompFrac, m_dPhaseCompFrac, m_dCompFrac_dCompDens,
+                                       phaseFlux, dPhaseFlux_dP, dPhaseFlux_dC, dPhaseFlux_dTrans,
+                                       compFlux, dCompFlux_dP, dCompFlux_dC, dCompFlux_dTrans );
+
           // 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, m_kernelFlags.isSet( KernelFlags::CheckPhasePresenceInGravity ),

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

@@ -56,9 +56,9 @@ enum class KernelFlags
   /// Flag indicating whether C1-PPU is used or not
   C1PPU = 1 << 5, // 32
   /// Flag indicating whether IHU is used or not
-  IHU = 1 << 6 // 64
-        /// Add more flags like that if needed:
-        // Flag8 = 1 << 7  //128
+  IHU = 1 << 6, // 64
+  /// Flag indicating whether HU 2-phase simplified version is used or not
+  HU2PH = 1 << 7 // 128
 };
 
 /******************************** FluxComputeKernelBase ********************************/
@@ -176,6 +176,81 @@ class FluxComputeKernelBase
   BitFlags< KernelFlags > const m_kernelFlags;
 };
 
+namespace helpers
+{
+template< typename VIEWTYPE >
+using ElementViewConst = ElementRegionManager::ElementViewConst< VIEWTYPE >;
+
+template< localIndex numComp, localIndex numFluxSupportPoints >
+GEOS_HOST_DEVICE
+static void calculateMeanDensity( localIndex const ip,
+                                  localIndex const (&seri)[numFluxSupportPoints],
+                                  localIndex const (&sesri)[numFluxSupportPoints],
+                                  localIndex const (&sei)[numFluxSupportPoints],
+                                  integer const checkPhasePresenceInGravity,
+                                  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_dPres)[numFluxSupportPoints], real64 (& dDensMean_dComp)[numFluxSupportPoints][numComp] )
+{
+  using Deriv = constitutive::multifluid::DerivativeOffset;
+
+  densMean = 0;
+  integer denom = 0;
+  real64 dDens_dC[numComp]{};
+  for( localIndex 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( checkPhasePresenceInGravity && !phaseExists )
+    {
+      dDensMean_dPres[i] = 0.0;
+      for( localIndex jc = 0; jc < numComp; ++jc )
+      {
+        dDensMean_dComp[i][jc] = 0.0;
+      }
+      continue;
+    }
+
+    // density
+    real64 const density = phaseMassDens[er][esr][ei][0][ip];
+    real64 const dDens_dPres = 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_dPres[i] = dDens_dPres;
+    for( localIndex jc = 0; jc < numComp; ++jc )
+    {
+      dDensMean_dComp[i][jc] = dDens_dC[jc];
+    }
+    denom++;
+  }
+  if( denom > 1 )
+  {
+    densMean /= denom;
+    for( localIndex i = 0; i < numFluxSupportPoints; ++i )
+    {
+      dDensMean_dPres[i] /= denom;
+      for( integer jc = 0; jc < numComp; ++jc )
+      {
+        dDensMean_dComp[i][jc] /= denom;
+      }
+    }
+  }
+}
+
+}
+
 } // namespace isothermalCompositionalMultiphaseFVMKernels
 
 } // namespace geos