Make sequential poromechanics to conserve mass (#2952)

Co-authored-by: Pavel Tomin <“[email protected]”>
Co-authored-by: Matteo Cusini <[email protected]>
Co-authored-by: Randolph R. Settgast <[email protected]>
Co-authored-by: Nicola Castelletto <[email protected]>

src/coreComponents/physicsSolvers/fluidFlow/CompositionalMultiphaseBase.cpp

@@ -359,10 +359,15 @@ void CompositionalMultiphaseBase::registerDataOnMesh( Group & meshBodies )
       subRegion.registerField< dPhaseMobility >( getName() ).
         reference().resizeDimension< 1, 2 >( m_numPhases, m_numComponents + 2 ); // dP, dT, dC
 
+      // needed for time step selector
       subRegion.registerField< phaseVolumeFraction_n >( getName() ).
         reference().resizeDimension< 1 >( m_numPhases );
-      subRegion.registerField< phaseMobility_n >( getName() ).
-        reference().resizeDimension< 1 >( m_numPhases );
+
+      subRegion.registerField< compAmount >( getName() ).
+        reference().resizeDimension< 1 >( m_numComponents );
+      subRegion.registerField< compAmount_n >( getName() ).
+        reference().resizeDimension< 1 >( m_numComponents );
+
     } );
 
     FaceManager & faceManager = mesh.getFaceManager();
@@ -697,6 +702,57 @@ void CompositionalMultiphaseBase::updateCapPressureModel( ObjectManagerBase & da
   }
 }
 
+void CompositionalMultiphaseBase::updateCompAmount( ElementSubRegionBase & subRegion ) const
+{
+  GEOS_MARK_FUNCTION;
+
+  string const & solidName = subRegion.template getReference< string >( viewKeyStruct::solidNamesString() );
+  CoupledSolidBase const & porousMaterial = getConstitutiveModel< CoupledSolidBase >( subRegion, solidName );
+  arrayView2d< real64 const > const porosity = porousMaterial.getPorosity();
+  arrayView1d< real64 const > const volume = subRegion.getElementVolume();
+  arrayView2d< real64 const, compflow::USD_COMP > const compDens = subRegion.getField< fields::flow::globalCompDensity >();
+  arrayView2d< real64, compflow::USD_COMP > const compAmount = subRegion.getField< fields::flow::compAmount >();
+
+  integer const numComp = m_numComponents;
+
+  forAll< parallelDevicePolicy<> >( subRegion.size(), [=] GEOS_HOST_DEVICE ( localIndex const ei )
+  {
+    for( integer ic = 0; ic < numComp; ++ic )
+    {
+      compAmount[ei][ic] = porosity[ei][0] * volume[ei] * compDens[ei][ic];
+    }
+  } );
+}
+
+void CompositionalMultiphaseBase::updateEnergy( ElementSubRegionBase & subRegion ) const
+{
+  GEOS_MARK_FUNCTION;
+
+  string const & solidName = subRegion.template getReference< string >( viewKeyStruct::solidNamesString() );
+  CoupledSolidBase const & porousMaterial = getConstitutiveModel< CoupledSolidBase >( subRegion, solidName );
+  arrayView2d< real64 const > const porosity = porousMaterial.getPorosity();
+  arrayView2d< real64 const > rockInternalEnergy = porousMaterial.getInternalEnergy();
+  arrayView1d< real64 const > const volume = subRegion.getElementVolume();
+  arrayView2d< real64 const, compflow::USD_PHASE > const phaseVolFrac = subRegion.getField< fields::flow::phaseVolumeFraction >();
+  string const & fluidName = getConstitutiveName< MultiFluidBase >( subRegion );
+  MultiFluidBase & fluid = subRegion.getConstitutiveModel< MultiFluidBase >( fluidName );
+  arrayView3d< real64 const, multifluid::USD_PHASE > const phaseDens = fluid.phaseDensity();
+  arrayView3d< real64 const, multifluid::USD_PHASE > const phaseInternalEnergy = fluid.phaseInternalEnergy();
+
+  arrayView1d< real64 > const energy = subRegion.getField< fields::flow::energy >();
+
+  integer const numPhases = m_numPhases;
+
+  forAll< parallelDevicePolicy<> >( subRegion.size(), [=] GEOS_HOST_DEVICE ( localIndex const ei )
+  {
+    energy[ei] = volume[ei] * (1.0 - porosity[ei][0]) * rockInternalEnergy[ei][0];
+    for( integer ip = 0; ip < numPhases; ++ip )
+    {
+      energy[ei] += volume[ei] * porosity[ei][0] * phaseVolFrac[ei][ip] * phaseDens[ei][0][ip] * phaseInternalEnergy[ei][0][ip];
+    }
+  } );
+}
+
 void CompositionalMultiphaseBase::updateSolidInternalEnergyModel( ObjectManagerBase & dataGroup ) const
 {
   arrayView1d< real64 const > const temp = dataGroup.getField< fields::flow::temperature >();
@@ -715,12 +771,13 @@ void CompositionalMultiphaseBase::updateSolidInternalEnergyModel( ObjectManagerB
                                       temp );
 }
 
-real64 CompositionalMultiphaseBase::updateFluidState( ObjectManagerBase & subRegion ) const
+real64 CompositionalMultiphaseBase::updateFluidState( ElementSubRegionBase & subRegion ) const
 {
   GEOS_MARK_FUNCTION;
 
   updateGlobalComponentFraction( subRegion );
   updateFluidModel( subRegion );
+  updateCompAmount( subRegion );
   real64 const maxDeltaPhaseVolFrac = updatePhaseVolumeFraction( subRegion );
   updateRelPermModel( subRegion );
   updatePhaseMobility( subRegion );
@@ -768,7 +825,6 @@ void CompositionalMultiphaseBase::initializeFluidState( MeshLevel & mesh,
         compDens[ei][ic] = totalDens[ei][0] * compFrac[ei][ic];
       }
     } );
-
   } );
 
   // with initial component densities defined - check if they need to be corrected to avoid zero diags etc
@@ -796,6 +852,7 @@ void CompositionalMultiphaseBase::initializeFluidState( MeshLevel & mesh,
     porousSolid.scaleReferencePorosity( netToGross );
     saveConvergedState( subRegion ); // necessary for a meaningful porosity update in sequential schemes
     updatePorosityAndPermeability( subRegion );
+    updateCompAmount( subRegion );
     updatePhaseVolumeFraction( subRegion );
 
     // Now, we initialize and update each constitutive model one by one
@@ -883,6 +940,8 @@ void CompositionalMultiphaseBase::initializeFluidState( MeshLevel & mesh,
       SolidInternalEnergy const & solidInternalEnergyMaterial =
         getConstitutiveModel< SolidInternalEnergy >( subRegion, solidInternalEnergyName );
       solidInternalEnergyMaterial.saveConvergedState();
+
+      updateEnergy( subRegion );
     }
 
     // Step 4.7: if the diffusion and/or dispersion is/are supported, initialize the two models
@@ -1251,19 +1310,13 @@ CompositionalMultiphaseBase::implicitStepSetup( real64 const & GEOS_UNUSED_PARAM
         updateSolidInternalEnergyModel( subRegion );
       }
 
-      // after the update, save the new saturation and phase mobilities
+      // after the update, save the new saturation
       arrayView2d< real64 const, compflow::USD_PHASE > const phaseVolFrac =
         subRegion.template getField< fields::flow::phaseVolumeFraction >();
       arrayView2d< real64, compflow::USD_PHASE > const phaseVolFrac_n =
         subRegion.template getField< fields::flow::phaseVolumeFraction_n >();
       phaseVolFrac_n.setValues< parallelDevicePolicy<> >( phaseVolFrac );
 
-      arrayView2d< real64 const, compflow::USD_PHASE > const phaseMob =
-        subRegion.template getField< fields::flow::phaseMobility >();
-      arrayView2d< real64, compflow::USD_PHASE > const phaseMob_n =
-        subRegion.template getField< fields::flow::phaseMobility_n >();
-      phaseMob_n.setValues< parallelDevicePolicy<> >( phaseMob );
-
     } );
   } );
 }
@@ -2381,6 +2434,13 @@ void CompositionalMultiphaseBase::saveConvergedState( ElementSubRegionBase & sub
   arrayView2d< real64, compflow::USD_COMP > const & compDens_n =
     subRegion.template getField< fields::flow::globalCompDensity_n >();
   compDens_n.setValues< parallelDevicePolicy<> >( compDens );
+
+  arrayView2d< real64 const, compflow::USD_COMP > const & compAmount =
+    subRegion.template getField< fields::flow::compAmount >();
+  arrayView2d< real64, compflow::USD_COMP > const & compAmount_n =
+    subRegion.template getField< fields::flow::compAmount_n >();
+  compAmount_n.setValues< parallelDevicePolicy<> >( compAmount );
+
   if( m_isFixedStressPoromechanicsUpdate )
   {
     arrayView2d< real64, compflow::USD_COMP > const & compDens_k =
@@ -2438,6 +2498,9 @@ void CompositionalMultiphaseBase::updateState( DomainPartition & domain )
 {
   GEOS_MARK_FUNCTION;
 
+  if( m_keepFlowVariablesConstantDuringInitStep )
+    return;
+
   real64 maxDeltaPhaseVolFrac = 0.0;
   forDiscretizationOnMeshTargets( domain.getMeshBodies(), [&]( string const &,
                                                                MeshLevel & mesh,
@@ -2456,6 +2519,7 @@ void CompositionalMultiphaseBase::updateState( DomainPartition & domain )
       if( m_isThermal )
       {
         updateSolidInternalEnergyModel( subRegion );
+        updateEnergy( subRegion );
       }
     } );
   } );

src/coreComponents/physicsSolvers/fluidFlow/CompositionalMultiphaseBase.hpp

@@ -133,6 +133,18 @@ class CompositionalMultiphaseBase : public FlowSolverBase
    */
   void updateCapPressureModel( ObjectManagerBase & dataGroup ) const;
 
+  /**
+   * @brief Update components mass/moles
+   * @param subRegion the subregion storing the required fields
+   */
+  void updateCompAmount( ElementSubRegionBase & subRegion ) const;
+
+  /**
+   * @brief Update energy
+   * @param subRegion the subregion storing the required fields
+   */
+  void updateEnergy( ElementSubRegionBase & subRegion ) const;
+
   /**
    * @brief Update all relevant solid internal energy models using current values of temperature
    * @param dataGroup the group storing the required fields
@@ -145,7 +157,7 @@ class CompositionalMultiphaseBase : public FlowSolverBase
    */
   virtual void updatePhaseMobility( ObjectManagerBase & dataGroup ) const = 0;
 
-  real64 updateFluidState( ObjectManagerBase & dataGroup ) const;
+  real64 updateFluidState( ElementSubRegionBase & subRegion ) const;
 
   virtual void saveConvergedState( ElementSubRegionBase & subRegion ) const override final;
 

src/coreComponents/physicsSolvers/fluidFlow/CompositionalMultiphaseBaseFields.hpp

@@ -119,6 +119,7 @@ DECLARE_FIELD( dPhaseMobility,
                NO_WRITE,
                "Derivative of phase volume fraction with respect to pressure, temperature, global component density" );
 
+// this is needed for time step selector
 DECLARE_FIELD( phaseVolumeFraction_n,
                "phaseVolumeFraction_n",
                array2dLayoutPhase,
@@ -127,14 +128,6 @@ DECLARE_FIELD( phaseVolumeFraction_n,
                WRITE_AND_READ,
                "Phase volume fraction at the previous converged time step" );
 
-DECLARE_FIELD( phaseMobility_n,
-               "phaseMobility_n",
-               array2dLayoutPhase,
-               0,
-               NOPLOT,
-               WRITE_AND_READ,
-               "Phase mobility at the previous converged time step" );
-
 DECLARE_FIELD( phaseOutflux,
                "phaseOutflux",
                array2dLayoutPhase,
@@ -175,6 +168,22 @@ DECLARE_FIELD( globalCompDensityScalingFactor,
                NO_WRITE,
                "Scaling factors for global component densities" );
 
+DECLARE_FIELD( compAmount,
+               "compAmount",
+               array2dLayoutComp,
+               0,
+               LEVEL_0,
+               WRITE_AND_READ,
+               "Component amount" );
+
+DECLARE_FIELD( compAmount_n,
+               "compAmount_n",
+               array2dLayoutComp,
+               0,
+               LEVEL_0,
+               WRITE_AND_READ,
+               "Component amount at the previous converged time step" );
+
 }
 
 }

src/coreComponents/physicsSolvers/fluidFlow/FlowSolverBase.cpp

@@ -160,6 +160,11 @@ void FlowSolverBase::registerDataOnMesh( Group & meshBodies )
       {
         subRegion.registerField< fields::flow::temperature_k >( getName() ); // needed for the fixed-stress porosity update
       }
+      if( m_isThermal )
+      {
+        subRegion.registerField< fields::flow::energy >( getName() );
+        subRegion.registerField< fields::flow::energy_n >( getName() );
+      }
     } );
 
     elemManager.forElementSubRegionsComplete< SurfaceElementSubRegion >( [&]( localIndex const,
@@ -218,6 +223,13 @@ void FlowSolverBase::saveConvergedState( ElementSubRegionBase & subRegion ) cons
   arrayView1d< real64 > const temp_n = subRegion.template getField< fields::flow::temperature_n >();
   temp_n.setValues< parallelDevicePolicy<> >( temp );
 
+  if( m_isThermal )
+  {
+    arrayView1d< real64 const > const energy = subRegion.template getField< fields::flow::energy >();
+    arrayView1d< real64 > const energy_n = subRegion.template getField< fields::flow::energy_n >();
+    energy_n.setValues< parallelDevicePolicy<> >( energy );
+  }
+
   if( m_isFixedStressPoromechanicsUpdate )
   {
     arrayView1d< real64 > const pres_k = subRegion.template getField< fields::flow::pressure_k >();

src/coreComponents/physicsSolvers/fluidFlow/FlowSolverBaseFields.hpp

@@ -96,14 +96,6 @@ DECLARE_FIELD( temperature,
                WRITE_AND_READ,
                "Temperature" );
 
-DECLARE_FIELD( faceTemperature,
-               "faceTemperature",
-               array1d< real64 >,
-               0,
-               LEVEL_0,
-               WRITE_AND_READ,
-               "Face temperature" );
-
 DECLARE_FIELD( temperature_n,
                "temperature_n",
                array1d< real64 >,
@@ -116,7 +108,7 @@ DECLARE_FIELD( temperature_k,
                "temperature_k",
                array1d< real64 >,
                0,
-               LEVEL_0,
+               NOPLOT,
                NO_WRITE,
                "Temperature at the previous sequential iteration" );
 
@@ -128,6 +120,14 @@ DECLARE_FIELD( initialTemperature,
                WRITE_AND_READ,
                "Initial temperature" );
 
+DECLARE_FIELD( faceTemperature,
+               "faceTemperature",
+               array1d< real64 >,
+               0,
+               LEVEL_0,
+               WRITE_AND_READ,
+               "Face temperature" );
+
 DECLARE_FIELD( netToGross,
                "netToGross",
                array1d< real64 >,
@@ -240,6 +240,38 @@ DECLARE_FIELD( temperatureScalingFactor,
                NO_WRITE,
                "Scaling factors for temperature" );
 
+DECLARE_FIELD( mass,
+               "mass",
+               array1d< real64 >,
+               0,
+               LEVEL_0,
+               WRITE_AND_READ,
+               "Mass" );
+
+DECLARE_FIELD( mass_n,
+               "mass_n",
+               array1d< real64 >,
+               0,
+               NOPLOT,
+               WRITE_AND_READ,
+               "Mass at the previous converged time step" );
+
+DECLARE_FIELD( energy,
+               "energy",
+               array1d< real64 >,
+               0,
+               LEVEL_0,
+               WRITE_AND_READ,
+               "Energy" );
+
+DECLARE_FIELD( energy_n,
+               "energy_n",
+               array1d< real64 >,
+               0,
+               NOPLOT,
+               WRITE_AND_READ,
+               "Energy at the previous converged time step" );
+
 }
 
 }