Merge remote-tracking branch 'origin/develop' into cusini/refactor/se…

…parate-friction-and-aperture
GEOS-DEV · Jul 1, 2024 · 0463ec5 · 0463ec5
2 parents 6072e5b + 32b28b4
commit 0463ec5
Show file tree

Hide file tree

Showing 56 changed files with 942 additions and 440 deletions.
diff --git a/.integrated_tests.yaml b/.integrated_tests.yaml
@@ -1,6 +1,6 @@
 baselines:
   bucket: geosx
-  baseline: integratedTests/baseline_integratedTests-pr3177-5828-eabb935
+  baseline: integratedTests/baseline_integratedTests-pr3196-5835-6563552
 
 allow_fail:
   all: ''

diff --git a/BASELINE_NOTES.md b/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 #3196 (2024-06-28)
+======================
+Added isLaggingFractureStencilWeightsUpdate to hydrofracture solve. Rebaseline because of the new input.
+
 PR #3177 (2024-06-28)
 ======================
 Added logLevel to TimeHistoryOutput. Rebaseline because of the new input flag.

diff --git a/src/coreComponents/mesh/FaceElementSubRegion.cpp b/src/coreComponents/mesh/FaceElementSubRegion.cpp
@@ -494,7 +494,6 @@ buildCollocatedEdgeBuckets( std::map< globalIndex, globalIndex > const & referen
                             arrayView2d< localIndex const > const edgeToNodes )
 {
   GEOS_ASSERT_EQ( edgeToNodes.size( 1 ), 2 );
-  static constexpr std::string_view nodeNotFound = "Internal error when trying to access the reference collocated node for global node {}.";
 
   // Checks if the node `gni` is handled as a collocated node on the curren rank.
   auto hasCollocatedNode = [&]( globalIndex const gni ) -> bool
@@ -526,6 +525,8 @@ buildCollocatedEdgeBuckets( std::map< globalIndex, globalIndex > const & referen
   std::map< std::pair< globalIndex, globalIndex >, std::set< localIndex > > collocatedEdgeBuckets;
   for( auto const & p: edgesIds )
   {
+    static constexpr std::string_view nodeNotFound = "Internal error when trying to access the reference collocated node for global node {}.";
+
     std::pair< globalIndex, globalIndex > const & nodes = p.first;
     localIndex const & edge = p.second;
 

diff --git a/src/coreComponents/physicsSolvers/SolverBase.hpp b/src/coreComponents/physicsSolvers/SolverBase.hpp
@@ -774,8 +774,6 @@ class SolverBase : public ExecutableGroup
    *  names set.
    */
   virtual void setConstitutiveNamesCallSuper( ElementSubRegionBase & subRegion ) const { GEOS_UNUSED_VAR( subRegion ); }
-  virtual void setConstitutiveNamesCallSuper( ParticleSubRegionBase & subRegion ) const { GEOS_UNUSED_VAR( subRegion ); } // particle
-                                                                                                                          // overload
 
   template< typename BASETYPE = constitutive::ConstitutiveBase, typename LOOKUP_TYPE >
   static BASETYPE const & getConstitutiveModel( dataRepository::Group const & dataGroup, LOOKUP_TYPE const & key );
@@ -840,7 +838,6 @@ class SolverBase : public ExecutableGroup
    *  names set.
    */
   virtual void setConstitutiveNames( ElementSubRegionBase & subRegion ) const { GEOS_UNUSED_VAR( subRegion ); }
-  virtual void setConstitutiveNames( ParticleSubRegionBase & subRegion ) const { GEOS_UNUSED_VAR( subRegion ); } // particle overload
 
   bool solveNonlinearSystem( real64 const & time_n,
                              real64 const & dt,

diff --git a/src/coreComponents/physicsSolvers/fluidFlow/FlowSolverBase.cpp b/src/coreComponents/physicsSolvers/fluidFlow/FlowSolverBase.cpp
@@ -91,7 +91,8 @@ FlowSolverBase::FlowSolverBase( string const & name,
   m_isThermal( 0 ),
   m_keepFlowVariablesConstantDuringInitStep( 0 ),
   m_isFixedStressPoromechanicsUpdate( false ),
-  m_isJumpStabilized( false )
+  m_isJumpStabilized( false ),
+  m_isLaggingFractureStencilWeightsUpdate( 0 )
 {
   this->registerWrapper( viewKeyStruct::isThermalString(), &m_isThermal ).
     setApplyDefaultValue( 0 ).
@@ -754,6 +755,8 @@ void FlowSolverBase::prepareStencilWeights( DomainPartition & domain ) const
     FiniteVolumeManager const & fvManager = numericalMethodManager.getFiniteVolumeManager();
     FluxApproximationBase const & fluxApprox = fvManager.getFluxApproximation( getDiscretizationName() );
     ElementRegionManager::ElementViewAccessor< arrayView1d< real64 const > > hydraulicAperture =
+      m_isLaggingFractureStencilWeightsUpdate ?
+      mesh.getElemManager().constructViewAccessor< array1d< real64 >, arrayView1d< real64 const > >( fields::flow::aperture0::key() ) :
       mesh.getElemManager().constructViewAccessor< array1d< real64 >, arrayView1d< real64 const > >( fields::flow::hydraulicAperture::key() );
 
     fluxApprox.forStencils< SurfaceElementStencil, FaceElementToCellStencil, EmbeddedSurfaceToCellStencil >( mesh, [&]( auto & stencil )

diff --git a/src/coreComponents/physicsSolvers/fluidFlow/FlowSolverBase.hpp b/src/coreComponents/physicsSolvers/fluidFlow/FlowSolverBase.hpp
@@ -157,6 +157,8 @@ class FlowSolverBase : public SolverBase
 
   virtual bool checkSequentialSolutionIncrements( DomainPartition & domain ) const override;
 
+  void enableLaggingFractureStencilWeightsUpdate(){ m_isLaggingFractureStencilWeightsUpdate = 1; };
+
 protected:
 
   /**
@@ -228,6 +230,8 @@ class FlowSolverBase : public SolverBase
 private:
   virtual void setConstitutiveNames( ElementSubRegionBase & subRegion ) const override;
 
+  // flag to determine whether or not to apply lagging update for the fracture stencil weights
+  integer m_isLaggingFractureStencilWeightsUpdate;
 
 };
 

diff --git a/src/coreComponents/physicsSolvers/multiphysics/HydrofractureSolver.cpp b/src/coreComponents/physicsSolvers/multiphysics/HydrofractureSolver.cpp
@@ -75,6 +75,11 @@ HydrofractureSolver< POROMECHANICS_SOLVER >::HydrofractureSolver( const string &
     setApplyDefaultValue( 0 ).
     setInputFlag( InputFlags::OPTIONAL );
 
+  registerWrapper( viewKeyStruct::isLaggingFractureStencilWeightsUpdateString(), &m_isLaggingFractureStencilWeightsUpdate ).
+    setApplyDefaultValue( 0 ).
+    setInputFlag( InputFlags::OPTIONAL ).
+    setDescription( "Flag to determine whether or not to apply lagging update for the fracture stencil weights. " );
+
   m_numResolves[0] = 0;
 
   // This may need to be different depending on whether poroelasticity is on or not.
@@ -112,6 +117,11 @@ void HydrofractureSolver< POROMECHANICS_SOLVER >::registerDataOnMesh( dataReposi
     } );
   } );
 #endif
+
+  if( m_isLaggingFractureStencilWeightsUpdate )
+  {
+    flowSolver()->enableLaggingFractureStencilWeightsUpdate();
+  }
 }
 
 template< typename POROMECHANICS_SOLVER >
@@ -848,13 +858,19 @@ void HydrofractureSolver< POROMECHANICS_SOLVER >::updateState( DomainPartition &
 
   Base::updateState( domain );
 
-  // remove the contribution of the hydraulic aperture from the stencil weights
-  flowSolver()->prepareStencilWeights( domain );
+  if( !m_isLaggingFractureStencilWeightsUpdate )
+  {
+    // remove the contribution of the hydraulic aperture from the stencil weights
+    flowSolver()->prepareStencilWeights( domain );
+  }
 
   updateHydraulicApertureAndFracturePermeability( domain );
 
-  // update the stencil weights using the updated hydraulic aperture
-  flowSolver()->updateStencilWeights( domain );
+  if( !m_isLaggingFractureStencilWeightsUpdate )
+  {
+    // update the stencil weights using the updated hydraulic aperture
+    flowSolver()->updateStencilWeights( domain );
+  }
 
   forDiscretizationOnMeshTargets( domain.getMeshBodies(), [&] ( string const &,
                                                                 MeshLevel & mesh,
@@ -872,6 +888,23 @@ void HydrofractureSolver< POROMECHANICS_SOLVER >::updateState( DomainPartition &
   } );
 }
 
+template< typename POROMECHANICS_SOLVER >
+void HydrofractureSolver< POROMECHANICS_SOLVER >::implicitStepComplete( real64 const & time_n,
+                                                                        real64 const & dt,
+                                                                        DomainPartition & domain )
+{
+  Base::implicitStepComplete( time_n, dt, domain );
+
+  if( m_isLaggingFractureStencilWeightsUpdate )
+  {
+    // remove the contribution of the hydraulic aperture from the stencil weights
+    flowSolver()->prepareStencilWeights( domain );
+
+    // update the stencil weights using the updated hydraulic aperture
+    flowSolver()->updateStencilWeights( domain );
+  }
+}
+
 template< typename POROMECHANICS_SOLVER >
 real64 HydrofractureSolver< POROMECHANICS_SOLVER >::setNextDt( real64 const & currentDt,
                                                                DomainPartition & domain )

diff --git a/src/coreComponents/physicsSolvers/multiphysics/HydrofractureSolver.hpp b/src/coreComponents/physicsSolvers/multiphysics/HydrofractureSolver.hpp
@@ -124,6 +124,10 @@ class HydrofractureSolver : public POROMECHANICS_SOLVER
 
   virtual void updateState( DomainPartition & domain ) override final;
 
+  virtual void implicitStepComplete( real64 const & time_n,
+                                     real64 const & dt,
+                                     DomainPartition & domain ) override final;
+
   /**@}*/
 
   void updateHydraulicApertureAndFracturePermeability( DomainPartition & domain );
@@ -170,6 +174,8 @@ class HydrofractureSolver : public POROMECHANICS_SOLVER
 
     constexpr static char const * useQuasiNewtonString() { return "useQuasiNewton"; }
 
+    static constexpr char const * isLaggingFractureStencilWeightsUpdateString() { return "isLaggingFractureStencilWeightsUpdate"; }
+
 #ifdef GEOSX_USE_SEPARATION_COEFFICIENT
     constexpr static char const * separationCoeff0String() { return "separationCoeff0"; }
     constexpr static char const * apertureAtFailureString() { return "apertureAtFailure"; }
@@ -243,6 +249,9 @@ class HydrofractureSolver : public POROMECHANICS_SOLVER
 
   integer m_useQuasiNewton;   // use Quasi-Newton (see https://arxiv.org/abs/2111.00264)
 
+  // flag to determine whether or not to apply lagging update for the fracture stencil weights
+  integer m_isLaggingFractureStencilWeightsUpdate;
+
 };
 
 ENUM_STRINGS( HydrofractureSolver< SinglePhasePoromechanics< SinglePhaseBase > >::InitializationType,

diff --git a/src/coreComponents/physicsSolvers/solidMechanics/SolidMechanicsMPM.cpp b/src/coreComponents/physicsSolvers/solidMechanics/SolidMechanicsMPM.cpp
@@ -342,8 +342,7 @@ void SolidMechanicsMPM::registerDataOnMesh( Group & meshBodies )
                                                                       [&]( localIndex const,
                                                                            ParticleSubRegionBase & subRegion )
       {
-        setConstitutiveNamesCallSuper( subRegion );
-        setConstitutiveNames( subRegion );
+        setParticlesConstitutiveNames( subRegion );
       } );
     }
 
@@ -1966,10 +1965,8 @@ void SolidMechanicsMPM::solverProfilingIf( std::string label, bool condition )
   }
 }
 
-void SolidMechanicsMPM::setConstitutiveNamesCallSuper( ParticleSubRegionBase & subRegion ) const
+void SolidMechanicsMPM::setParticlesConstitutiveNames( ParticleSubRegionBase & subRegion ) const
 {
-  SolverBase::setConstitutiveNamesCallSuper( subRegion );
-
   subRegion.registerWrapper< string >( viewKeyStruct::solidMaterialNamesString() ).
     setPlotLevel( PlotLevel::NOPLOT ).
     setRestartFlags( RestartFlags::NO_WRITE ).
@@ -1980,11 +1977,6 @@ void SolidMechanicsMPM::setConstitutiveNamesCallSuper( ParticleSubRegionBase & s
   GEOS_ERROR_IF( solidMaterialName.empty(), GEOS_FMT( "SolidBase model not found on subregion {}", subRegion.getName() ) );
 }
 
-void SolidMechanicsMPM::setConstitutiveNames( ParticleSubRegionBase & subRegion ) const
-{
-  GEOS_UNUSED_VAR( subRegion );
-}
-
 real64 SolidMechanicsMPM::computeNeighborList( ParticleManager & particleManager )
 {
   // Time this function

diff --git a/src/coreComponents/physicsSolvers/solidMechanics/SolidMechanicsMPM.hpp b/src/coreComponents/physicsSolvers/solidMechanics/SolidMechanicsMPM.hpp
@@ -363,8 +363,6 @@ class SolidMechanicsMPM : public SolverBase
 protected:
   virtual void postProcessInput() override final;
 
-  virtual void setConstitutiveNamesCallSuper( ParticleSubRegionBase & subRegion ) const override;
-
   std::vector< array2d< localIndex > > m_mappedNodes; // mappedNodes[subregion index][particle index][node index]. dims = {# of subregions,
                                                       // # of particles, # of nodes a particle on the subregion maps to}
   std::vector< array2d< real64 > > m_shapeFunctionValues; // mappedNodes[subregion][particle][nodal shape function value]. dims = {# of
@@ -457,7 +455,7 @@ class SolidMechanicsMPM : public SolverBase
     }
   };
 
-  virtual void setConstitutiveNames( ParticleSubRegionBase & subRegion ) const override;
+  void setParticlesConstitutiveNames( ParticleSubRegionBase & subRegion ) const;
 };
 
 ENUM_STRINGS( SolidMechanicsMPM::TimeIntegrationOption,

diff --git a/src/coreComponents/schema/docs/AcousticElasticSEM.rst b/src/coreComponents/schema/docs/AcousticElasticSEM.rst
@@ -11,6 +11,7 @@ initialDt                 real64             1e+99    Initial time-step value re
 logLevel                  integer            0        Log level                                                                                                                                                                                                                                                                                                                
 name                      groupName          required A name is required for any non-unique nodes                                                                                                                                                                                                                                                                              
 targetRegions             groupNameRef_array required Allowable regions that the solver may be applied to. Note that this does not indicate that the solver will be applied to these regions, only that allocation will occur such that the solver may be applied to these regions. The decision about what regions this solver will beapplied to rests in the EventManager.   
+writeLinearSystem         integer            0        Write matrix, rhs, solution to screen ( = 1) or file ( = 2).                                                                                                                                                                                                                                                             
 LinearSolverParameters    node               unique   :ref:`XML_LinearSolverParameters`                                                                                                                                                                                                                                                                                        
 NonlinearSolverParameters node               unique   :ref:`XML_NonlinearSolverParameters`                                                                                                                                                                                                                                                                                     
 ========================= ================== ======== ========================================================================================================================================================================================================================================================================================================================