Freeze porosity when doing poromechanics initialization (#3056)

GEOS-DEV · Apr 19, 2024 · 5408b7d · 5408b7d
1 parent 65064b0
commit 5408b7d
Show file tree

Hide file tree

Showing 14 changed files with 45 additions and 4 deletions.
diff --git a/integratedTests b/integratedTests
diff --git a/src/coreComponents/constitutive/solid/PorousSolid.hpp b/src/coreComponents/constitutive/solid/PorousSolid.hpp
@@ -82,6 +82,7 @@ class PorousSolidUpdates : public CoupledSolidUpdates< SOLID_TYPE, BiotPorosity,
                                        real64 ( & dTotalStress_dPressure )[6],
                                        real64 ( & dTotalStress_dTemperature )[6],
                                        DiscretizationOps & stiffness,
+                                       integer const performStressInitialization,
                                        real64 & porosity,
                                        real64 & porosity_n,
                                        real64 & dPorosity_dVolStrain,
@@ -116,6 +117,15 @@ class PorousSolidUpdates : public CoupledSolidUpdates< SOLID_TYPE, BiotPorosity,
                      dPorosity_dPressure,
                      dPorosity_dTemperature );
 
+    // skip porosity update when doing poromechanics initialization
+    if( performStressInitialization )
+    {
+      porosity = porosityInit;
+      dPorosity_dVolStrain = 0.0;
+      dPorosity_dPressure = 0.0;
+      dPorosity_dTemperature = 0.0;
+    }
+
     // Save the derivative of solid density wrt pressure for the computation of the body force
     dSolidDensity_dPressure = m_porosityUpdate.dGrainDensity_dPressure();
   }

diff --git a/src/coreComponents/physicsSolvers/multiphysics/MultiphasePoromechanics.cpp b/src/coreComponents/physicsSolvers/multiphysics/MultiphasePoromechanics.cpp
@@ -224,6 +224,7 @@ void MultiphasePoromechanics< FLOW_SOLVER >::assembleElementBasedTerms( real64 c
                                                                                                     this->flowSolver()->numFluidComponents(),
                                                                                                     this->flowSolver()->numFluidPhases(),
                                                                                                     this->flowSolver()->useTotalMassEquation(),
+                                                                                                    this->m_performStressInitialization,
                                                                                                     FlowSolverBase::viewKeyStruct::fluidNamesString() );
     }
     else
@@ -242,6 +243,7 @@ void MultiphasePoromechanics< FLOW_SOLVER >::assembleElementBasedTerms( real64 c
                                                                                       this->flowSolver()->numFluidPhases(),
                                                                                       this->flowSolver()->useSimpleAccumulation(),
                                                                                       this->flowSolver()->useTotalMassEquation(),
+                                                                                      this->m_performStressInitialization,
                                                                                       FlowSolverBase::viewKeyStruct::fluidNamesString() );
     }
   } );

diff --git a/src/coreComponents/physicsSolvers/multiphysics/PoromechanicsSolver.hpp b/src/coreComponents/physicsSolvers/multiphysics/PoromechanicsSolver.hpp
@@ -389,7 +389,8 @@ class PoromechanicsSolver : public CoupledSolver< FLOW_SOLVER, MECHANICS_SOLVER
     }
 
     /// After the solid mechanics solver
-    if( solverType == static_cast< integer >( SolverType::SolidMechanics ) )
+    if( solverType == static_cast< integer >( SolverType::SolidMechanics )
+        && !m_performStressInitialization ) // do not update during poromechanics initialization
     {
       // compute the average of the mean total stress increment over quadrature points
       averageMeanTotalStressIncrement( domain );

diff --git a/src/coreComponents/physicsSolvers/multiphysics/SinglePhasePoromechanics.cpp b/src/coreComponents/physicsSolvers/multiphysics/SinglePhasePoromechanics.cpp
@@ -223,6 +223,7 @@ void SinglePhasePoromechanics< FLOW_SOLVER, MECHANICS_SOLVER >::assembleElementB
                                                                                                      localRhs,
                                                                                                      dt,
                                                                                                      flowDofKey,
+                                                                                                     this->m_performStressInitialization,
                                                                                                      FlowSolverBase::viewKeyStruct::fluidNamesString() );
     }
     else
@@ -237,6 +238,7 @@ void SinglePhasePoromechanics< FLOW_SOLVER, MECHANICS_SOLVER >::assembleElementB
                                                                                        localRhs,
                                                                                        dt,
                                                                                        flowDofKey,
+                                                                                       this->m_performStressInitialization,
                                                                                        FlowSolverBase::viewKeyStruct::fluidNamesString() );
     }
   } );

diff --git a/src/coreComponents/physicsSolvers/multiphysics/SinglePhasePoromechanicsEmbeddedFractures.hpp b/src/coreComponents/physicsSolvers/multiphysics/SinglePhasePoromechanicsEmbeddedFractures.hpp
@@ -185,6 +185,7 @@ real64 SinglePhasePoromechanicsEmbeddedFractures::assemblyLaunch( MeshLevel & me
                                 dt,
                                 gravityVectorData,
                                 flowDofKey,
+                                m_performStressInitialization,
                                 FlowSolverBase::viewKeyStruct::fluidNamesString() );
 
   real64 const maxForce =

diff --git a/...reComponents/physicsSolvers/multiphysics/poromechanicsKernels/MultiphasePoromechanics.hpp b/...reComponents/physicsSolvers/multiphysics/poromechanicsKernels/MultiphasePoromechanics.hpp
@@ -102,6 +102,7 @@ class MultiphasePoromechanics :
                            localIndex const numPhases,
                            integer const useSimpleAccumulation,
                            integer const useTotalMassEquation,
+                           integer const performStressInitialization,
                            string const fluidModelKey );
 
   //*****************************************************************************
@@ -354,6 +355,8 @@ class MultiphasePoromechanics :
 
   /// Use total mass equation flag
   integer const m_useTotalMassEquation;
+
+  integer const m_performStressInitialization;
 };
 
 using MultiphasePoromechanicsKernelFactory =
@@ -369,6 +372,7 @@ using MultiphasePoromechanicsKernelFactory =
                                 localIndex const,
                                 integer const,
                                 integer const,
+                                integer const,
                                 string const >;
 
 /**

diff --git a/...ponents/physicsSolvers/multiphysics/poromechanicsKernels/MultiphasePoromechanics_impl.hpp b/...ponents/physicsSolvers/multiphysics/poromechanicsKernels/MultiphasePoromechanics_impl.hpp
@@ -54,6 +54,7 @@ MultiphasePoromechanics( NodeManager const & nodeManager,
                          localIndex const numPhases,
                          integer const useSimpleAccumulation,
                          integer const useTotalMassEquation,
+                         integer const performStressInitialization,
                          string const fluidModelKey ):
   Base( nodeManager,
         edgeManager,
@@ -79,7 +80,8 @@ MultiphasePoromechanics( NodeManager const & nodeManager,
   m_numComponents( numComponents ),
   m_numPhases( numPhases ),
   m_useSimpleAccumulation( useSimpleAccumulation ),
-  m_useTotalMassEquation( useTotalMassEquation )
+  m_useTotalMassEquation( useTotalMassEquation ),
+  m_performStressInitialization( performStressInitialization )
 {
   GEOS_ERROR_IF_GT_MSG( m_numComponents, maxNumComponents,
                         "MultiphasePoromechanics solver allows at most " <<
@@ -162,6 +164,7 @@ smallStrainUpdate( localIndex const k,
                                                        stack.dTotalStress_dPressure,
                                                        stack.dTotalStress_dTemperature,
                                                        stack.stiffness,
+                                                       m_performStressInitialization,
                                                        porosity,
                                                        porosity_n,
                                                        dPorosity_dVolStrain,

diff --git a/...eComponents/physicsSolvers/multiphysics/poromechanicsKernels/SinglePhasePoromechanics.hpp b/...eComponents/physicsSolvers/multiphysics/poromechanicsKernels/SinglePhasePoromechanics.hpp
@@ -94,6 +94,7 @@ class SinglePhasePoromechanics :
                             real64 const inputDt,
                             real64 const (&gravityVector)[3],
                             string const inputFlowDofKey,
+                            integer const performStressInitialization,
                             string const fluidModelKey );
 
   //*****************************************************************************
@@ -257,6 +258,7 @@ class SinglePhasePoromechanics :
   /// Derivative of fluid density wrt pressure
   arrayView2d< real64 const > const m_dFluidDensity_dPressure;
 
+  integer const m_performStressInitialization;
 };
 
 using SinglePhasePoromechanicsKernelFactory =
@@ -268,6 +270,7 @@ using SinglePhasePoromechanicsKernelFactory =
                                 real64 const,
                                 real64 const (&)[3],
                                 string const,
+                                integer const,
                                 string const >;
 
 /**

diff --git a/...onents/physicsSolvers/multiphysics/poromechanicsKernels/SinglePhasePoromechanics_impl.hpp b/...onents/physicsSolvers/multiphysics/poromechanicsKernels/SinglePhasePoromechanics_impl.hpp
@@ -50,6 +50,7 @@ SinglePhasePoromechanics( NodeManager const & nodeManager,
                           real64 const inputDt,
                           real64 const (&gravityVector)[3],
                           string const inputFlowDofKey,
+                          integer const performStressInitialization,
                           string const fluidModelKey ):
   Base( nodeManager,
         edgeManager,
@@ -68,7 +69,8 @@ SinglePhasePoromechanics( NodeManager const & nodeManager,
         fluidModelKey ),
   m_fluidDensity( elementSubRegion.template getConstitutiveModel< constitutive::SingleFluidBase >( elementSubRegion.template getReference< string >( fluidModelKey ) ).density() ),
   m_fluidDensity_n( elementSubRegion.template getConstitutiveModel< constitutive::SingleFluidBase >( elementSubRegion.template getReference< string >( fluidModelKey ) ).density_n() ),
-  m_dFluidDensity_dPressure( elementSubRegion.template getConstitutiveModel< constitutive::SingleFluidBase >( elementSubRegion.template getReference< string >( fluidModelKey ) ).dDensity_dPressure() )
+  m_dFluidDensity_dPressure( elementSubRegion.template getConstitutiveModel< constitutive::SingleFluidBase >( elementSubRegion.template getReference< string >( fluidModelKey ) ).dDensity_dPressure() ),
+  m_performStressInitialization( performStressInitialization )
 {}
 
 template< typename SUBREGION_TYPE,
@@ -100,6 +102,7 @@ smallStrainUpdate( localIndex const k,
                                                        stack.dTotalStress_dPressure,
                                                        stack.dTotalStress_dTemperature,
                                                        stack.stiffness,
+                                                       m_performStressInitialization,
                                                        porosity,
                                                        porosity_n,
                                                        dPorosity_dVolStrain,

diff --git a/...nents/physicsSolvers/multiphysics/poromechanicsKernels/ThermalMultiphasePoromechanics.hpp b/...nents/physicsSolvers/multiphysics/poromechanicsKernels/ThermalMultiphasePoromechanics.hpp
@@ -76,6 +76,7 @@ class ThermalMultiphasePoromechanics :
   using Base::m_numComponents;
   using Base::m_numPhases;
   using Base::m_useTotalMassEquation;
+  using Base::m_performStressInitialization;
   using Base::m_solidDensity;
   using Base::m_fluidPhaseMassDensity;
   using Base::m_dFluidPhaseMassDensity;
@@ -110,6 +111,7 @@ class ThermalMultiphasePoromechanics :
                                   localIndex const numComponents,
                                   localIndex const numPhases,
                                   integer const useTotalMassEquation,
+                                  integer const performStressInitialization,
                                   string const fluidModelKey );
 
   /**
@@ -341,6 +343,7 @@ using ThermalMultiphasePoromechanicsKernelFactory =
                                 localIndex const,
                                 localIndex const,
                                 integer const,
+                                integer const,
                                 string const >;
 
 } // namespace thermalporomechanicsKernels

diff --git a/.../physicsSolvers/multiphysics/poromechanicsKernels/ThermalMultiphasePoromechanics_impl.hpp b/.../physicsSolvers/multiphysics/poromechanicsKernels/ThermalMultiphasePoromechanics_impl.hpp
@@ -52,6 +52,7 @@ ThermalMultiphasePoromechanics( NodeManager const & nodeManager,
                                 localIndex const numComponents,
                                 localIndex const numPhases,
                                 integer const useTotalMassEquation,
+                                integer const performStressInitialization,
                                 string const fluidModelKey ):
   Base( nodeManager,
         edgeManager,
@@ -71,6 +72,7 @@ ThermalMultiphasePoromechanics( NodeManager const & nodeManager,
         numPhases,
         false, // do not use simple accumulation form
         useTotalMassEquation,
+        performStressInitialization,
         fluidModelKey ),
   m_rockInternalEnergy_n( inputConstitutiveType.getInternalEnergy_n() ),
   m_rockInternalEnergy( inputConstitutiveType.getInternalEnergy() ),
@@ -137,6 +139,7 @@ smallStrainUpdate( localIndex const k,
                                                        stack.dTotalStress_dPressure,
                                                        stack.dTotalStress_dTemperature,
                                                        stack.stiffness,
+                                                       m_performStressInitialization,
                                                        porosity,
                                                        porosity_n,
                                                        dPorosity_dVolStrain,

diff --git a/...ents/physicsSolvers/multiphysics/poromechanicsKernels/ThermalSinglePhasePoromechanics.hpp b/...ents/physicsSolvers/multiphysics/poromechanicsKernels/ThermalSinglePhasePoromechanics.hpp
@@ -70,6 +70,7 @@ class ThermalSinglePhasePoromechanics :
   using Base::m_solidDensity;
   using Base::m_flowDofNumber;
   using Base::m_dt;
+  using Base::m_performStressInitialization;
 
   /**
    * @brief Constructor
@@ -90,6 +91,7 @@ class ThermalSinglePhasePoromechanics :
                                    real64 const inputDt,
                                    real64 const (&gravityVector)[3],
                                    string const inputFlowDofKey,
+                                   integer const performStressInitialization,
                                    string const fluidModelKey );
 
   //*****************************************************************************
@@ -305,6 +307,7 @@ using ThermalSinglePhasePoromechanicsKernelFactory =
                                 real64 const,
                                 real64 const (&)[3],
                                 string const,
+                                integer const,
                                 string const >;
 
 } // namespace thermalPoromechanicsKernels

diff --git a/...physicsSolvers/multiphysics/poromechanicsKernels/ThermalSinglePhasePoromechanics_impl.hpp b/...physicsSolvers/multiphysics/poromechanicsKernels/ThermalSinglePhasePoromechanics_impl.hpp
@@ -46,6 +46,7 @@ ThermalSinglePhasePoromechanics( NodeManager const & nodeManager,
                                  real64 const inputDt,
                                  real64 const (&gravityVector)[3],
                                  string const inputFlowDofKey,
+                                 integer const performStressInitialization,
                                  string const fluidModelKey ):
   Base( nodeManager,
         edgeManager,
@@ -61,6 +62,7 @@ ThermalSinglePhasePoromechanics( NodeManager const & nodeManager,
         inputDt,
         gravityVector,
         inputFlowDofKey,
+        performStressInitialization,
         fluidModelKey ),
   m_dFluidDensity_dTemperature( elementSubRegion.template getConstitutiveModel< constitutive::SingleFluidBase >( elementSubRegion.template getReference< string >(
                                                                                                                    fluidModelKey ) ).dDensity_dTemperature() ),
@@ -121,6 +123,7 @@ smallStrainUpdate( localIndex const k,
                                                        stack.dTotalStress_dPressure,
                                                        stack.dTotalStress_dTemperature,
                                                        stack.stiffness,
+                                                       m_performStressInitialization,
                                                        porosity,
                                                        porosity_n,
                                                        dPorosity_dVolStrain,
+ −		poromechanics/PoroElastic_staircase_co2_3d_fim_01/0to33_restart_000000033.root
+2 −2		poromechanics/PoroElastic_staircase_co2_3d_fim_01/0to33_restart_000000033/rank_0000000.hdf5
+ −		poromechanics/PoroElastic_staircase_co2_3d_fim_04/0to33_restart_000000033.root
+2 −2		poromechanics/PoroElastic_staircase_co2_3d_fim_04/0to33_restart_000000033/rank_0000000.hdf5
+2 −2		poromechanics/PoroElastic_staircase_co2_3d_fim_04/0to33_restart_000000033/rank_0000001.hdf5
+2 −2		poromechanics/PoroElastic_staircase_co2_3d_fim_04/0to33_restart_000000033/rank_0000002.hdf5
+2 −2		poromechanics/PoroElastic_staircase_co2_3d_fim_04/0to33_restart_000000033/rank_0000003.hdf5
+ −		poromechanics/PoroElastic_staircase_co2_3d_sequential_01/0to33_restart_000000033.root
+2 −2		poromechanics/PoroElastic_staircase_co2_3d_sequential_01/0to33_restart_000000033/rank_0000000.hdf5
+ −		poromechanics/PoroElastic_staircase_co2_3d_sequential_04/0to33_restart_000000033.root
+2 −2		poromechanics/PoroElastic_staircase_co2_3d_sequential_04/0to33_restart_000000033/rank_0000000.hdf5
+2 −2		poromechanics/PoroElastic_staircase_co2_3d_sequential_04/0to33_restart_000000033/rank_0000001.hdf5
+2 −2		poromechanics/PoroElastic_staircase_co2_3d_sequential_04/0to33_restart_000000033/rank_0000002.hdf5
+2 −2		poromechanics/PoroElastic_staircase_co2_3d_sequential_04/0to33_restart_000000033/rank_0000003.hdf5
+ −		poromechanics/PoroElastic_staircase_singlephase_3d_fim_01/0to11_restart_000000011.root
+2 −2		poromechanics/PoroElastic_staircase_singlephase_3d_fim_01/0to11_restart_000000011/rank_0000000.hdf5
+ −		poromechanics/PoroElastic_staircase_singlephase_3d_fim_04/0to11_restart_000000011.root
+2 −2		poromechanics/PoroElastic_staircase_singlephase_3d_fim_04/0to11_restart_000000011/rank_0000000.hdf5
+2 −2		poromechanics/PoroElastic_staircase_singlephase_3d_fim_04/0to11_restart_000000011/rank_0000001.hdf5
+2 −2		poromechanics/PoroElastic_staircase_singlephase_3d_fim_04/0to11_restart_000000011/rank_0000002.hdf5
+2 −2		poromechanics/PoroElastic_staircase_singlephase_3d_fim_04/0to11_restart_000000011/rank_0000003.hdf5
+ −		poromechanics/PoroElastic_staircase_singlephase_3d_sequential_01/0to11_restart_000000011.root
+2 −2		poromechanics/PoroElastic_staircase_singlephase_3d_sequential_01/0to11_restart_000000011/rank_0000000.hdf5
+ −		poromechanics/PoroElastic_staircase_singlephase_3d_sequential_04/0to11_restart_000000011.root
+2 −2		poromechanics/PoroElastic_staircase_singlephase_3d_sequential_04/0to11_restart_000000011/rank_0000000.hdf5
+2 −2		poromechanics/PoroElastic_staircase_singlephase_3d_sequential_04/0to11_restart_000000011/rank_0000001.hdf5
+2 −2		poromechanics/PoroElastic_staircase_singlephase_3d_sequential_04/0to11_restart_000000011/rank_0000002.hdf5
+2 −2		poromechanics/PoroElastic_staircase_singlephase_3d_sequential_04/0to11_restart_000000011/rank_0000003.hdf5
+ −		thermoPoromechanics/ThermoPoroElastic_staircase_co2_smoke_01/0to33_restart_000000033.root
+2 −2		thermoPoromechanics/ThermoPoroElastic_staircase_co2_smoke_01/0to33_restart_000000033/rank_0000000.hdf5
+ −		thermoPoromechanics/ThermoPoroElastic_staircase_co2_smoke_04/0to33_restart_000000033.root
+2 −2		thermoPoromechanics/ThermoPoroElastic_staircase_co2_smoke_04/0to33_restart_000000033/rank_0000000.hdf5
+2 −2		thermoPoromechanics/ThermoPoroElastic_staircase_co2_smoke_04/0to33_restart_000000033/rank_0000001.hdf5
+2 −2		thermoPoromechanics/ThermoPoroElastic_staircase_co2_smoke_04/0to33_restart_000000033/rank_0000002.hdf5
+2 −2		thermoPoromechanics/ThermoPoroElastic_staircase_co2_smoke_04/0to33_restart_000000033/rank_0000003.hdf5