feat: make possible to use solid mechanics solver to perform poromechanics initialization

.integrated_tests.yaml

@@ -1,6 +1,6 @@
 baselines:
   bucket: geosx
-  baseline: integratedTests/baseline_integratedTests-pr2125-10859-1db8623
+  baseline: integratedTests/baseline_integratedTests-pr3396-10875-13b594f
 
 allow_fail:
   all: ''

BASELINE_NOTES.md

@@ -6,6 +6,10 @@ This file is designed to track changes to the integrated test baselines.
 Any developer who updates the baseline ID in the .integrated_tests.yaml file is expected to create an entry in this file with the pull request number, date, and their justification for rebaselining.
 These notes should be in reverse-chronological order, and use the following time format: (YYYY-MM-DD).
 
+PR #3396 (2024-03-21)
+=====================
+Use solid mechanics solver directly to perform poromechanics initialization.
+
 PR #2125 (2024-03-20)
 =====================
 Phase-field nucleation model.

src/coreComponents/physicsSolvers/fluidFlow/FlowSolverBase.hpp

@@ -168,6 +168,8 @@ class FlowSolverBase : public PhysicsSolverBase
     GEOS_ERROR( "Poroelastic fluxes with conforming fractures not yet implemented." );
   }
 
+  void initializeState( DomainPartition & domain );
+
   virtual void initializeFluidState( MeshLevel & mesh, string_array const & regionNames ) { GEOS_UNUSED_VAR( mesh, regionNames ); }
 
   virtual void initializeThermalState( MeshLevel & mesh, string_array const & regionNames ) { GEOS_UNUSED_VAR( mesh, regionNames ); }
@@ -236,8 +238,6 @@ class FlowSolverBase : public PhysicsSolverBase
 
   virtual void initializePostInitialConditionsPreSubGroups() override;
 
-  void initializeState( DomainPartition & domain );
-
   virtual void computeHydrostaticEquilibrium( DomainPartition & domain ) { GEOS_UNUSED_VAR( domain ); }
 
   void initializePorosityAndPermeability( MeshLevel & mesh, string_array const & regionNames );

src/coreComponents/physicsSolvers/multiphysics/CoupledReservoirAndWellsBase.hpp

@@ -231,6 +231,8 @@ class CoupledReservoirAndWellsBase : public CoupledSolver< RESERVOIR_SOLVER, WEL
     }
   }
 
+  void initializeState( DomainPartition & domain ) const { return reservoirSolver()->initializeState( domain ); }
+
   void
   assembleFluxTerms( real64 const dt,
                      DomainPartition const & domain,

src/coreComponents/physicsSolvers/multiphysics/MultiphasePoromechanics.hpp

@@ -41,7 +41,7 @@ class MultiphasePoromechanics : public PoromechanicsSolver< FLOW_SOLVER, MECHANI
   using Base::m_stabilizationType;
   using Base::m_stabilizationRegionNames;
   using Base::m_stabilizationMultiplier;
-  using Base::getLogLevel;
+  using Base::updateBulkDensity;
 
   /**
    * @brief main constructor for MultiphasePoromechanics Objects

src/coreComponents/physicsSolvers/multiphysics/PoromechanicsInitialization.cpp

@@ -21,7 +21,6 @@
 
 #include "events/tasks/TasksManager.hpp"
 #include "physicsSolvers/PhysicsSolverManager.hpp"
-#include "physicsSolvers/fluidFlow/SinglePhaseBase.hpp"
 #include "physicsSolvers/solidMechanics/SolidMechanicsStatistics.hpp"
 #include "physicsSolvers/multiphysics/MultiphasePoromechanics.hpp"
 #include "physicsSolvers/multiphysics/MultiphasePoromechanicsConformingFractures.hpp"
@@ -87,8 +86,9 @@ postInputInitialization()
     TasksManager & tasksManager = problemManager.getGroup< TasksManager >( "Tasks" );
 
     GEOS_THROW_IF( !tasksManager.hasGroup( m_solidMechanicsStatisticsName ),
-                   GEOS_FMT( "{}: statistics task named {} not found",
+                   GEOS_FMT( "{}: {} task named {} not found",
                              getWrapperDataContext( viewKeyStruct::solidMechanicsStatisticsNameString() ),
+                             SolidMechanicsStatistics::catalogName(),
                              m_solidMechanicsStatisticsName ),
                    InputError );
 
@@ -116,7 +116,9 @@ execute( real64 const time_n,
 
   m_solidMechanicsStateResetTask.execute( time_n, dt, cycleNumber, eventCounter, eventProgress, domain );
 
-  m_poromechanicsSolver->execute( time_n, dt, cycleNumber, eventCounter, eventProgress, domain );
+  m_poromechanicsSolver->flowSolver()->initializeState( domain );
+  m_poromechanicsSolver->updateBulkDensity( domain );
+  m_poromechanicsSolver->solidMechanicsSolver()->execute( time_n, dt, cycleNumber, eventCounter, eventProgress, domain );
 
   GEOS_LOG_LEVEL_RANK_0( logInfo::SolverInitialization, GEOS_FMT( "Task `{}`: at time {}s, physics solver `{}` has completed stress initialization",//1
                                                                   getName(), time_n + dt, m_poromechanicsSolverName ) );

src/coreComponents/physicsSolvers/multiphysics/PoromechanicsInitialization.hpp

@@ -91,8 +91,6 @@ class PoromechanicsInitialization : public TaskBase
 
   void postInputInitialization() override;
 
-//  void registerDataOnMesh( Group & meshBodies ) override;
-
   /// Name of the poromechanics solver
   string m_poromechanicsSolverName;
 

src/coreComponents/physicsSolvers/multiphysics/PoromechanicsSolver.hpp

@@ -81,18 +81,14 @@ class PoromechanicsSolver : public CoupledSolver< FLOW_SOLVER, MECHANICS_SOLVER
   PoromechanicsSolver( const string & name,
                        dataRepository::Group * const parent )
     : Base( name, parent ),
-    m_isThermal( 0 )
+    m_isThermal( 0 ),
+    m_performStressInitialization( false )
   {
     this->registerWrapper( viewKeyStruct::isThermalString(), &m_isThermal ).
       setApplyDefaultValue( 0 ).
       setInputFlag( dataRepository::InputFlags::OPTIONAL ).
       setDescription( "Flag indicating whether the problem is thermal or not. Set isThermal=\"1\" to enable the thermal coupling" );
 
-    this->registerWrapper( viewKeyStruct::performStressInitializationString(), &m_performStressInitialization ).
-      setApplyDefaultValue( false ).
-      setInputFlag( dataRepository::InputFlags::FALSE ).
-      setDescription( "Flag to indicate that the solver is going to perform stress initialization" );
-
     this->registerWrapper( viewKeyStruct::stabilizationTypeString(), &m_stabilizationType ).
       setInputFlag( dataRepository::InputFlags::OPTIONAL ).
       setDescription( "StabilizationType. Options are:\n" +
@@ -220,12 +216,9 @@ class PoromechanicsSolver : public CoupledSolver< FLOW_SOLVER, MECHANICS_SOLVER
           setRestartFlags( dataRepository::RestartFlags::NO_WRITE ).
           setSizedFromParent( 0 );
 
-        if( this->getNonlinearSolverParameters().m_couplingType == NonlinearSolverParameters::CouplingType::Sequential )
-        {
-          // register the bulk density for use in the solid mechanics solver
-          // ideally we would resize it here as well, but the solid model name is not available yet (see below)
-          subRegion.registerField< fields::poromechanics::bulkDensity >( this->getName() );
-        }
+        // register the bulk density for use in the solid mechanics solver
+        // ideally we would resize it here as well, but the solid model name is not available yet (see below)
+        subRegion.registerField< fields::poromechanics::bulkDensity >( this->getName() );
 
         if( m_stabilizationType == stabilization::StabilizationType::Global || m_stabilizationType == stabilization::StabilizationType::Local )
         {
@@ -302,9 +295,10 @@ class PoromechanicsSolver : public CoupledSolver< FLOW_SOLVER, MECHANICS_SOLVER
    * @brief Utility function to set the stress initialization flag
    * @param[in] performStressInitialization true if the solver has to initialize stress, false otherwise
    */
-  void setStressInitialization( integer const performStressInitialization )
+  void setStressInitialization( bool const performStressInitialization )
   {
     m_performStressInitialization = performStressInitialization;
+    solidMechanicsSolver()->setStressInitialization( performStressInitialization );
   }
 
   struct viewKeyStruct : Base::viewKeyStruct
@@ -315,9 +309,6 @@ class PoromechanicsSolver : public CoupledSolver< FLOW_SOLVER, MECHANICS_SOLVER
     /// Flag to indicate that the simulation is thermal
     constexpr static char const * isThermalString() { return "isThermal"; }
 
-    /// Flag to indicate that the solver is going to perform stress initialization
-    constexpr static char const * performStressInitializationString() { return "performStressInitialization"; }
-
     /// Type of pressure stabilization
     constexpr static char const * stabilizationTypeString() {return "stabilizationType"; }
 
@@ -394,6 +385,23 @@ class PoromechanicsSolver : public CoupledSolver< FLOW_SOLVER, MECHANICS_SOLVER
 
   }
 
+  void updateBulkDensity( DomainPartition & domain )
+  {
+    this->template forDiscretizationOnMeshTargets( domain.getMeshBodies(), [&]( string const &,
+                                                                                MeshLevel & mesh,
+                                                                                string_array const & regionNames )
+    {
+      mesh.getElemManager().forElementSubRegions< CellElementSubRegion >( regionNames, [&]( localIndex const,
+                                                                                            auto & subRegion )
+      {
+        // update the bulk density
+        // TODO: ideally, we would not recompute the bulk density, but a more general "rhs" containing the body force and the
+        // pressure/temperature terms
+        updateBulkDensity( subRegion );
+      } );
+    } );
+  }
+
 protected:
 
   template< typename CONSTITUTIVE_BASE,
@@ -561,20 +569,7 @@ class PoromechanicsSolver : public CoupledSolver< FLOW_SOLVER, MECHANICS_SOLVER
     /// After the flow solver
     if( solverType == static_cast< integer >( SolverType::Flow ) )
     {
-      this->template forDiscretizationOnMeshTargets( domain.getMeshBodies(), [&]( string const &,
-                                                                                  MeshLevel & mesh,
-                                                                                  string_array const & regionNames )
-      {
-
-        mesh.getElemManager().forElementSubRegions< CellElementSubRegion >( regionNames, [&]( localIndex const,
-                                                                                              auto & subRegion )
-        {
-          // update the bulk density
-          // TODO: ideally, we would not recompute the bulk density, but a more general "rhs" containing the body force and the
-          // pressure/temperature terms
-          updateBulkDensity( subRegion );
-        } );
-      } );
+      updateBulkDensity( domain );
     }
 
     /// After the solid mechanics solver
@@ -668,7 +663,7 @@ class PoromechanicsSolver : public CoupledSolver< FLOW_SOLVER, MECHANICS_SOLVER
   integer m_isThermal;
 
   /// Flag to indicate that the solver is going to perform stress initialization
-  integer m_performStressInitialization;
+  bool m_performStressInitialization;
 
   /// Type of stabilization used
   stabilization::StabilizationType m_stabilizationType;

src/coreComponents/physicsSolvers/multiphysics/SinglePhasePoromechanics.hpp

@@ -41,6 +41,7 @@ class SinglePhasePoromechanics : public PoromechanicsSolver< FLOW_SOLVER, MECHAN
   using Base::m_stabilizationType;
   using Base::m_stabilizationRegionNames;
   using Base::m_stabilizationMultiplier;
+  using Base::updateBulkDensity;
 
   /**
    * @brief main constructor for SinglePhasePoromechanics objects

src/coreComponents/physicsSolvers/solidMechanics/SolidMechanicsLagrangianFEM.cpp

@@ -60,7 +60,8 @@ SolidMechanicsLagrangianFEM::SolidMechanicsLagrangianFEM( const string & name,
   m_maxForce( 0.0 ),
   m_maxNumResolves( 10 ),
   m_strainTheory( 0 ),
-  m_isFixedStressPoromechanicsUpdate( false )
+  m_isFixedStressPoromechanicsUpdate( false ),
+  m_performStressInitialization( false )
 {
 
   registerWrapper( viewKeyStruct::newmarkGammaString(), &m_newmarkGamma ).
@@ -1006,6 +1007,8 @@ void SolidMechanicsLagrangianFEM::setupSystem( DomainPartition & domain,
                                                ParallelVector & solution,
                                                bool const setSparsity )
 {
+  GEOS_LOG( "SolidMechanicsLagrangianFEM::setupSystem" );
+
   GEOS_MARK_FUNCTION;
   PhysicsSolverBase::setupSystem( domain, dofManager, localMatrix, rhs, solution, setSparsity );
 
@@ -1021,25 +1024,6 @@ void SolidMechanicsLagrangianFEM::setupSystem( DomainPartition & domain,
     arrayView1d< globalIndex const > const
     dofNumber = nodeManager.getReference< globalIndex_array >( dofManager.getKey( solidMechanics::totalDisplacement::key() ) );
 
-    if( m_contactRelationName != viewKeyStruct::noContactRelationNameString() )
-    {
-      ElementRegionManager const & elemManager = mesh.getElemManager();
-      string_array allFaceElementRegions;
-      elemManager.forElementRegions< SurfaceElementRegion >( [&]( SurfaceElementRegion const & elemRegion )
-      {
-        allFaceElementRegions.emplace_back( elemRegion.getName() );
-      } );
-
-      finiteElement::
-        fillSparsity< FaceElementSubRegion,
-                      solidMechanicsLagrangianFEMKernels::ImplicitSmallStrainQuasiStatic >( mesh,
-                                                                                            allFaceElementRegions,
-                                                                                            this->getDiscretizationName(),
-                                                                                            dofNumber,
-                                                                                            dofManager.rankOffset(),
-                                                                                            sparsityPattern );
-
-    }
     finiteElement::
       fillSparsity< CellElementSubRegion,
                     solidMechanicsLagrangianFEMKernels::ImplicitSmallStrainQuasiStatic >( mesh,
@@ -1072,7 +1056,7 @@ void SolidMechanicsLagrangianFEM::assembleSystem( real64 const GEOS_UNUSED_PARAM
                                                                 MeshLevel & mesh,
                                                                 string_array const & regionNames )
   {
-    if( m_isFixedStressPoromechanicsUpdate )
+    if( m_isFixedStressPoromechanicsUpdate || m_performStressInitialization )
     {
       set< string > poromechanicsRegions;
       set< string > mechanicsRegions;

src/coreComponents/physicsSolvers/solidMechanics/SolidMechanicsLagrangianFEM.hpp

@@ -278,6 +278,15 @@ class SolidMechanicsLagrangianFEM : public PhysicsSolverBase
     return m_rigidBodyModes;
   }
 
+  /*
+   * @brief Utility function to set the stress initialization flag
+   * @param[in] performStressInitialization true if the solver has to initialize stress, false otherwise
+   */
+  void setStressInitialization( bool const performStressInitialization )
+  {
+    m_performStressInitialization = performStressInitialization;
+  }
+
 protected:
   virtual void postInputInitialization() override;
 
@@ -293,8 +302,11 @@ class SolidMechanicsLagrangianFEM : public PhysicsSolverBase
   real64 m_maxForce = 0.0;
   integer m_maxNumResolves;
   integer m_strainTheory;
-//  MPI_iCommData m_iComm;
+
+  /// Flag to indicate that the solver is running in fixed stress (sequential) mode
   bool m_isFixedStressPoromechanicsUpdate;
+  /// Flag to indicate that the solver is going to perform stress initialization
+  bool m_performStressInitialization;
 
   /// Rigid body modes
   array1d< ParallelVector > m_rigidBodyModes;

src/coreComponents/physicsSolvers/solidMechanics/contact/SolidMechanicsLagrangeContact.cpp

@@ -640,7 +640,7 @@ void SolidMechanicsLagrangeContact::assembleSystem( real64 const time,
   assembleContact( domain, dofManager, localMatrix, localRhs );
 
   // for sequential: add (fixed) pressure force contribution into residual (no derivatives)
-  if( m_isFixedStressPoromechanicsUpdate )
+  if( m_isFixedStressPoromechanicsUpdate || m_performStressInitialization )
   {
     forDiscretizationOnMeshTargets( domain.getMeshBodies(), [&]( string const &,
                                                                  MeshLevel const & mesh,