refactor: mixture density computation in multiphase poromechanics sol…

…vers (#3342)

* Simplifying body force computation

src/coreComponents/physicsSolvers/multiphysics/MultiphasePoromechanics.cpp

@@ -296,7 +296,7 @@ void MultiphasePoromechanics< FLOW_SOLVER, MECHANICS_SOLVER >::updateBulkDensity
   // update the bulk density
   poromechanicsKernels::
     MultiphaseBulkDensityKernelFactory::
-    createAndLaunch< parallelDevicePolicy<> >( this->flowSolver()->numFluidPhases(),
+    createAndLaunch< parallelDevicePolicy<> >( this->flowSolver()->numFluidComponents(),
                                                fluid,
                                                solid,
                                                subRegion );

src/coreComponents/physicsSolvers/multiphysics/poromechanicsKernels/MultiphasePoromechanics.hpp

@@ -342,8 +342,10 @@ class MultiphasePoromechanics :
   arrayView3d< real64 const, compflow::USD_COMP_DC > const m_dGlobalCompFraction_dGlobalCompDensity;
 
   // Views on component densities
+  bool m_useMass;
   arrayView2d< real64 const, compflow::USD_COMP > m_compDens;
   arrayView2d< real64 const, compflow::USD_COMP > m_compDens_n;
+  arrayView1d< real64 const > m_componentMolarWeight;
 
   /// Number of components
   localIndex const m_numComponents;
@@ -386,20 +388,21 @@ class MultiphaseBulkDensityKernel
 
   /**
    * @brief Constructor
-   * @param[in] numPhases the number of fluid phases
+   * @param[in] numComponents the number of fluid components (species)
    * @param[in] fluid the fluid model
    * @param[in] solid the porous solid model
    * @param[in] subRegion the element subregion
    */
-  MultiphaseBulkDensityKernel( integer const numPhases,
+  MultiphaseBulkDensityKernel( integer const numComponents,
                                constitutive::MultiFluidBase const & fluid,
                                constitutive::CoupledSolidBase const & solid,
                                ElementSubRegionBase & subRegion )
-    : m_numPhases( numPhases ),
+    : m_numComponents( numComponents ),
     m_bulkDensity( subRegion.getField< fields::poromechanics::bulkDensity >() ),
-    m_fluidPhaseVolFrac( subRegion.getField< fields::flow::phaseVolumeFraction >() ),
+    m_compDens( subRegion.getField< fields::flow::globalCompDensity >() ),
+    m_useMass( fluid.getMassFlag()),
     m_rockDensity( solid.getDensity() ),
-    m_fluidPhaseDensity( fluid.phaseMassDensity() ),
+    m_componentMolarWeight( fluid.componentMolarWeights() ),
     m_porosity( solid.getPorosity() )
   {}
 
@@ -413,12 +416,12 @@ class MultiphaseBulkDensityKernel
                 localIndex const q ) const
   {
     m_bulkDensity[ei][q] = 0.0;
-    for( integer ip = 0; ip < m_numPhases; ++ip )
+    for( integer ic = 0; ic < m_numComponents; ++ic )
     {
-      m_bulkDensity[ei][q] += m_fluidPhaseVolFrac[ei][ip] * m_fluidPhaseDensity[ei][q][ip];
+      m_bulkDensity[ei][q] = m_bulkDensity[ei][q] + ( m_useMass ? m_compDens[ei][ic] : m_compDens[ei][ic] * m_componentMolarWeight[ic] );
     }
     m_bulkDensity[ei][q] *= m_porosity[ei][q];
-    m_bulkDensity[ei][q] += ( 1 - m_porosity[ei][q] ) * m_rockDensity[ei][q];
+    m_bulkDensity[ei][q] = m_bulkDensity[ei][q] + ( 1 - m_porosity[ei][q] ) * m_rockDensity[ei][q];
   }
 
   /**
@@ -446,20 +449,23 @@ class MultiphaseBulkDensityKernel
 
 protected:
 
-  // number of fluid phases
-  integer const m_numPhases;
+  // number of fluid components (species)
+  integer const m_numComponents;
 
   // the bulk density
   arrayView2d< real64 > const m_bulkDensity;
 
-  // the fluid phase saturation
-  arrayView2d< real64 const, compflow::USD_PHASE > const m_fluidPhaseVolFrac;
+  // the fluid component densities
+  arrayView2d< real64 const, compflow::USD_PHASE > const m_compDens;
+
+  // the flag for selecting mass formulation instead of molar formulation
+  bool m_useMass;
 
   // the rock density
   arrayView2d< real64 const > const m_rockDensity;
 
-  // the fluid density
-  arrayView3d< real64 const, constitutive::multifluid::USD_PHASE > const m_fluidPhaseDensity;
+  // the molar component weights
+  arrayView1d< real64 const > m_componentMolarWeight;
 
   // the porosity
   arrayView2d< real64 const > const m_porosity;
@@ -476,19 +482,19 @@ class MultiphaseBulkDensityKernelFactory
   /**
    * @brief Create a new kernel and launch
    * @tparam POLICY the policy used in the RAJA kernel
-   * @param[in] numPhases number of phases
+   * @param[in] numComponents number of phases components (species)
    * @param[in] fluid the fluid model
    * @param[in] solid the porous solid model
    * @param[in] subRegion the element subregion
    */
   template< typename POLICY >
   static void
-  createAndLaunch( integer const numPhases,
+  createAndLaunch( integer const numComponents,
                    constitutive::MultiFluidBase const & fluid,
                    constitutive::CoupledSolidBase const & solid,
                    ElementSubRegionBase & subRegion )
   {
-    MultiphaseBulkDensityKernel kernel( numPhases, fluid, solid, subRegion );
+    MultiphaseBulkDensityKernel kernel( numComponents, fluid, solid, subRegion );
     MultiphaseBulkDensityKernel::launch< POLICY >( subRegion.size(),
                                                    subRegion.getField< fields::poromechanics::bulkDensity >().size( 1 ),
                                                    kernel );

src/coreComponents/physicsSolvers/multiphysics/poromechanicsKernels/MultiphasePoromechanics_impl.hpp

@@ -104,6 +104,9 @@ MultiphasePoromechanics( NodeManager const & nodeManager,
 
     m_fluidPhaseMassDensity = fluid.phaseMassDensity();
     m_dFluidPhaseMassDensity = fluid.dPhaseMassDensity();
+
+    m_useMass = fluid.getMassFlag();
+    m_componentMolarWeight = fluid.componentMolarWeights();
   }
 }
 
@@ -214,48 +217,26 @@ computeBodyForce( localIndex const k,
                   StackVariables & stack,
                   FUNC && bodyForceKernelOp ) const
 {
-  using Deriv = constitutive::multifluid::DerivativeOffset;
-
   GEOS_UNUSED_VAR( dPorosity_dTemperature );
 
-  arraySlice1d< real64 const, constitutive::multifluid::USD_PHASE - 2 > const phaseMassDensity = m_fluidPhaseMassDensity[k][q];
-  arraySlice2d< real64 const, constitutive::multifluid::USD_PHASE_DC - 2 > const dPhaseMassDensity = m_dFluidPhaseMassDensity[k][q];
-  arraySlice1d< real64 const, compflow::USD_PHASE - 1 > const phaseVolFrac = m_fluidPhaseVolFrac[k];
-  arraySlice2d< real64 const, compflow::USD_PHASE_DC - 1 > const dPhaseVolFrac = m_dFluidPhaseVolFrac[k];
-  arraySlice2d< real64 const, compflow::USD_COMP_DC - 1 > const dGlobalCompFrac_dGlobalCompDensity = m_dGlobalCompFraction_dGlobalCompDensity[k];
-
   // Step 1: compute fluid total mass density and its derivatives
 
   real64 totalMassDensity = 0.0;
-  real64 dTotalMassDensity_dPressure = 0.0;
   real64 dTotalMassDensity_dComponents[maxNumComponents]{};
-  real64 dPhaseMassDensity_dComponents[maxNumComponents]{};
 
-  for( integer ip = 0; ip < m_numPhases; ++ip )
+  arraySlice1d< real64 const, compflow::USD_COMP - 1 > const compDens = m_compDens[k];
+  for( integer ic = 0; ic < m_numComponents; ++ic )
   {
-    totalMassDensity += phaseVolFrac( ip ) * phaseMassDensity( ip );
-    dTotalMassDensity_dPressure += dPhaseVolFrac( ip, Deriv::dP ) * phaseMassDensity( ip )
-                                   + phaseVolFrac( ip ) * dPhaseMassDensity( ip, Deriv::dP );
-
-    applyChainRule( m_numComponents,
-                    dGlobalCompFrac_dGlobalCompDensity,
-                    dPhaseMassDensity[ip],
-                    dPhaseMassDensity_dComponents,
-                    Deriv::dC );
-    for( integer jc = 0; jc < m_numComponents; ++jc )
-    {
-      dTotalMassDensity_dComponents[jc] += dPhaseVolFrac( ip, Deriv::dC+jc ) * phaseMassDensity( ip )
-                                           + phaseVolFrac( ip ) * dPhaseMassDensity_dComponents[jc];
-    }
+    totalMassDensity = totalMassDensity + ( m_useMass ? compDens[ic] : compDens[ic] * m_componentMolarWeight[ic] );
+    dTotalMassDensity_dComponents[ic] = m_useMass ? 1.0 : m_componentMolarWeight[ic];
   }
 
   // Step 2: compute mixture density as an average between total mass density and solid density
 
   real64 const mixtureDensity = ( 1.0 - porosity ) * m_solidDensity( k, q ) + porosity * totalMassDensity;
   real64 const dMixtureDens_dVolStrainIncrement = dPorosity_dVolStrain * ( -m_solidDensity( k, q ) + totalMassDensity );
   real64 const dMixtureDens_dPressure = dPorosity_dPressure * ( -m_solidDensity( k, q ) + totalMassDensity )
-                                        + ( 1.0 - porosity ) * dSolidDensity_dPressure
-                                        + porosity * dTotalMassDensity_dPressure;
+                                        + ( 1.0 - porosity ) * dSolidDensity_dPressure;
   LvArray::tensorOps::scale< maxNumComponents >( dTotalMassDensity_dComponents, porosity );
 
   // Step 3: finally, get the body force

src/coreComponents/physicsSolvers/multiphysics/poromechanicsKernels/ThermalMultiphasePoromechanics_impl.hpp

@@ -188,8 +188,6 @@ computeBodyForce( localIndex const k,
                   real64 const & dSolidDensity_dPressure,
                   StackVariables & stack ) const
 {
-  using Deriv = constitutive::multifluid::DerivativeOffset;
-
   Base::computeBodyForce( k, q,
                           porosity,
                           dPorosity_dVolStrain,
@@ -201,27 +199,11 @@ computeBodyForce( localIndex const k,
   {
     GEOS_UNUSED_VAR( mixtureDensity );
 
-    arraySlice1d< real64 const, constitutive::multifluid::USD_PHASE - 2 > const phaseMassDensity = m_fluidPhaseMassDensity[k][q];
-    arraySlice2d< real64 const, constitutive::multifluid::USD_PHASE_DC - 2 > const dPhaseMassDensity = m_dFluidPhaseMassDensity[k][q];
-    arraySlice1d< real64 const, compflow::USD_PHASE - 1 > const phaseVolFrac = m_fluidPhaseVolFrac[k];
-    arraySlice2d< real64 const, compflow::USD_PHASE_DC - 1 > const dPhaseVolFrac = m_dFluidPhaseVolFrac[k];
-
-    // Step 1: compute fluid total mass density and its derivatives
-
-    real64 dTotalMassDensity_dTemperature = 0.0;
-
-    for( integer ip = 0; ip < m_numPhases; ++ip )
-    {
-      dTotalMassDensity_dTemperature += dPhaseVolFrac( ip, Deriv::dT ) * phaseMassDensity( ip )
-                                        + phaseVolFrac( ip ) * dPhaseMassDensity( ip, Deriv::dT );
-    }
-
-    // Step 2: compute the derivative of the bulk density (an average between total mass density and solid density) wrt temperature
-
-    real64 const dMixtureDens_dTemperature = dPorosity_dTemperature * ( -m_solidDensity( k, q ) + totalMassDensity )
-                                             + porosity * dTotalMassDensity_dTemperature;
+    // Step 1: compute the derivative of the mixture density (an average between total mass density and solid density) wrt temperature
+    //         TODO include solid density derivative with respect to temperature
+    real64 const dMixtureDens_dTemperature = dPorosity_dTemperature * ( -m_solidDensity( k, q ) + totalMassDensity );
 
-    // Step 3: finally, get the body force
+    // Step 2: finally, get the body force
 
     LvArray::tensorOps::scaledCopy< 3 >( stack.dBodyForce_dTemperature, m_gravityVector, dMixtureDens_dTemperature );