feat: Configuration loop acceleration for SolidMechanicsLagrangeContact

src/coreComponents/physicsSolvers/PhysicsSolverBase.cpp

@@ -810,7 +810,7 @@ real64 PhysicsSolverBase::nonlinearImplicitStep( real64 const & time_n,
 
       if( isNewtonConverged )
       {
-        isConfigurationLoopConverged = updateConfiguration( domain );
+        isConfigurationLoopConverged = updateConfiguration( domain, configurationLoopIter );
 
         if( isConfigurationLoopConverged )
         {
@@ -1342,7 +1342,8 @@ void PhysicsSolverBase::updateState( DomainPartition & GEOS_UNUSED_PARAM( domain
   GEOS_ERROR( "PhysicsSolverBase::updateState called!. Should be overridden." );
 }
 
-bool PhysicsSolverBase::updateConfiguration( DomainPartition & GEOS_UNUSED_PARAM( domain ) )
+bool PhysicsSolverBase::updateConfiguration( DomainPartition & GEOS_UNUSED_PARAM( domain ),
+                                             integer const GEOS_UNUSED_PARAM( configurationLoopIter ) )
 {
   return true;
 }

src/coreComponents/physicsSolvers/PhysicsSolverBase.hpp

@@ -579,9 +579,11 @@ class PhysicsSolverBase : public ExecutableGroup
   /**
    * @brief updates the configuration (if needed) based on the state after a converged Newton loop.
    * @param domain the domain containing the mesh and fields
+   * @param configurationLoopIter current configuration iteration number
    * @return a bool that states whether the configuration used to solve the nonlinear loop is still valid or not.
    */
-  virtual bool updateConfiguration( DomainPartition & domain );
+  virtual bool updateConfiguration( DomainPartition & domain,
+                                    integer configurationLoopIter );
 
   /**
    * @brief

src/coreComponents/physicsSolvers/contact/SolidMechanicsAugmentedLagrangianContact.cpp

@@ -746,7 +746,8 @@ void SolidMechanicsAugmentedLagrangianContact::updateState( DomainPartition & do
   GEOS_UNUSED_VAR( domain );
 }
 
-bool SolidMechanicsAugmentedLagrangianContact::updateConfiguration( DomainPartition & domain )
+bool SolidMechanicsAugmentedLagrangianContact::updateConfiguration( DomainPartition & domain,
+                                                                    integer const GEOS_UNUSED_PARAM( configurationLoopIter ) )
 {
   GEOS_MARK_FUNCTION;
 

src/coreComponents/physicsSolvers/contact/SolidMechanicsAugmentedLagrangianContact.hpp

@@ -88,7 +88,8 @@ class SolidMechanicsAugmentedLagrangianContact : public ContactSolverBase
 
   void updateState( DomainPartition & domain ) override final;
 
-  virtual bool updateConfiguration( DomainPartition & domain ) override final;
+  virtual bool updateConfiguration( DomainPartition & domain,
+                                    integer configurationLoopIter ) override final;
 
 
   /**

src/coreComponents/physicsSolvers/contact/SolidMechanicsEmbeddedFractures.cpp

@@ -786,7 +786,8 @@ void SolidMechanicsEmbeddedFractures::updateState( DomainPartition & domain )
   } );
 }
 
-bool SolidMechanicsEmbeddedFractures::updateConfiguration( DomainPartition & domain )
+bool SolidMechanicsEmbeddedFractures::updateConfiguration( DomainPartition & domain,
+                                                           integer const GEOS_UNUSED_PARAM( configurationLoopIter ) )
 {
   int hasConfigurationConverged = true;
 

src/coreComponents/physicsSolvers/contact/SolidMechanicsEmbeddedFractures.hpp

@@ -118,7 +118,8 @@ class SolidMechanicsEmbeddedFractures : public ContactSolverBase
                         real64 const dt,
                         DomainPartition & domain );
 
-  virtual bool updateConfiguration( DomainPartition & domain ) override final;
+  virtual bool updateConfiguration( DomainPartition & domain,
+                                    integer configurationLoopIter ) override final;
 
   bool useStaticCondensation() const { return m_useStaticCondensation; }
 

src/coreComponents/physicsSolvers/contact/SolidMechanicsLagrangeContact.cpp

@@ -72,6 +72,11 @@ SolidMechanicsLagrangeContact::SolidMechanicsLagrangeContact( const string & nam
     setApplyDefaultValue( 1.0 ).
     setDescription( "It be used to increase the scale of the stabilization entries. A value < 1.0 results in larger entries in the stabilization matrix." );
 
+  registerWrapper( viewKeyStruct::useLocalYieldAccelerationString(), &m_useLocalYieldAcceleration ).
+    setInputFlag( InputFlags::OPTIONAL ).
+    setApplyDefaultValue( 0 ).
+    setDescription( "Flag to enable local acceleration for yield (to accelerate configuration loop convergence)." );
+
   addLogLevel< logInfo::Configuration >();
 
   LinearSolverParameters & linSolParams = m_linearSolverParameters.get();
@@ -171,6 +176,12 @@ void SolidMechanicsLagrangeContact::initializePreSubGroups()
     } );
   }
 
+  if( m_useLocalYieldAcceleration )
+  {
+    GEOS_LOG_LEVEL_RANK_0( 1, GEOS_FMT( "{}: local yield acceleration enabled", getName() ) );
+    initializeAccelerationVariables( domain );
+  }
+
 }
 
 void SolidMechanicsLagrangeContact::setupSystem( DomainPartition & domain,
@@ -2202,7 +2213,8 @@ bool SolidMechanicsLagrangeContact::resetConfigurationToDefault( DomainPartition
   return false;
 }
 
-bool SolidMechanicsLagrangeContact::updateConfiguration( DomainPartition & domain )
+bool SolidMechanicsLagrangeContact::updateConfiguration( DomainPartition & domain,
+                                                         integer const configurationLoopIter )
 {
   GEOS_MARK_FUNCTION;
 
@@ -2261,6 +2273,7 @@ bool SolidMechanicsLagrangeContact::updateConfiguration( DomainPartition & domai
             else if( traction[kfe][0] > normalTractionTolerance[kfe] )
             {
               fractureState[kfe] = FractureState::Open;
+
               if( getLogLevel() >= 10 )
                 GEOS_LOG( GEOS_FMT( "{}: {} -> {}: traction = {}, normalTractionTolerance = {}",
                                     kfe, originalFractureState, fractureState[kfe],
@@ -2275,6 +2288,9 @@ bool SolidMechanicsLagrangeContact::updateConfiguration( DomainPartition & domai
                 frictionWrapper.computeLimitTangentialTractionNorm( traction[kfe][0],
                                                                     dLimitTangentialTractionNorm_dTraction );
 
+              // store to use in acceleration when enabled
+              real64 const currentTau_unscaled = currentTau;
+
               if( originalFractureState == FractureState::Stick && currentTau >= limitTau )
               {
                 currentTau *= (1.0 - m_slidingCheckTolerance);
@@ -2297,6 +2313,11 @@ bool SolidMechanicsLagrangeContact::updateConfiguration( DomainPartition & domai
               else
               {
                 fractureState[kfe] = FractureState::Stick;
+
+                if( m_useLocalYieldAcceleration )
+                {
+                  tryLocalYieldAcceleration( configurationLoopIter, kfe, currentTau_unscaled, limitTau, currentTau, fractureState[kfe] );
+                }
               }
               if( getLogLevel() >= 10 && fractureState[kfe] != originalFractureState )
                 GEOS_LOG( GEOS_FMT( "{}: {} -> {}: currentTau = {}, limitTau = {}",
@@ -2329,6 +2350,90 @@ bool SolidMechanicsLagrangeContact::updateConfiguration( DomainPartition & domai
   return changedArea <= m_nonlinearSolverParameters.m_configurationTolerance * totalArea;
 }
 
+void SolidMechanicsLagrangeContact::tryLocalYieldAcceleration( integer const configurationLoopIter,
+                                                               localIndex const kfe,
+                                                               real64 const currentTau_unscaled,
+                                                               real64 const limitTau,
+                                                               real64 const currentTau,
+                                                               integer & fractureState )
+{
+  using namespace fields::contact;
+
+  if( configurationLoopIter == 0 )
+  {
+    m_x0[kfe] = currentTau_unscaled - limitTau;
+  }
+  else if( configurationLoopIter == 1 )
+  {
+    m_x1_tilde[kfe] = currentTau_unscaled - limitTau;
+    m_x1[kfe] = m_x1_tilde[kfe];
+    m_omega0[kfe] = 1.0;
+  }
+  else
+  {
+    // only apply acceleration if within a fraction of the limitTau
+    real64 acceleration_buffer = (limitTau - currentTau) / limitTau;
+
+    if( acceleration_buffer > 0.1 / (configurationLoopIter - 1) )
+    {
+      // do not apply acceleration, just update previous values
+      m_x0[kfe] = m_x1_tilde[kfe];
+      m_x1_tilde[kfe] = currentTau_unscaled - limitTau;
+      m_x1[kfe] = m_x1_tilde[kfe];
+    }
+    else
+    {
+      // apply acceleration
+      m_x2_tilde[kfe] = currentTau_unscaled - limitTau;
+
+      m_omega1[kfe] = -1.0 * m_omega0[kfe] * (m_x1_tilde[kfe] - m_x0[kfe]) / ((m_x1_tilde[kfe] - m_x0[kfe]) - (m_x2_tilde[kfe] - m_x1[kfe]));
+
+      m_x2[kfe] = (1.0 - m_omega1[kfe]) * m_x1[kfe] + m_omega1[kfe] * m_x2_tilde[kfe];
+
+      real64 acc_currentTau_unscaled = m_x2[kfe] + limitTau;
+      real64 acc_currentTau_scaled = acc_currentTau_unscaled * (1.0 - m_slidingCheckTolerance);
+
+      if( acc_currentTau_scaled > limitTau )
+      {
+        fractureState = FractureState::NewSlip;
+      }
+      else
+      {
+        m_x0[kfe] = m_x1[kfe];
+        m_x1[kfe] = m_x2[kfe];
+        m_x1_tilde[kfe] = m_x2_tilde[kfe];
+        m_omega0[kfe] = m_omega1[kfe];
+      }
+    }
+  }
+}
+
+void SolidMechanicsLagrangeContact::initializeAccelerationVariables( DomainPartition & domain )
+{
+  // calculate total number of face elements
+  localIndex total_size = 0;
+  forDiscretizationOnMeshTargets( domain.getMeshBodies(), [&] ( string const &,
+                                                                MeshLevel & mesh,
+                                                                arrayView1d< string const > const & regionNames )
+  {
+    ElementRegionManager & elemManager = mesh.getElemManager();
+    elemManager.forElementSubRegions< FaceElementSubRegion >( regionNames, [&]( localIndex const,
+                                                                                FaceElementSubRegion & subRegion )
+    {
+      total_size += subRegion.size();
+    } );
+  } );
+
+  // allocate arrays
+  m_x0.resize( total_size );
+  m_x1.resize( total_size );
+  m_x1_tilde.resize( total_size );
+  m_x2.resize( total_size );
+  m_x2_tilde.resize( total_size );
+  m_omega0.resize( total_size );
+  m_omega1.resize( total_size );
+}
+
 bool SolidMechanicsLagrangeContact::isFractureAllInStickCondition( DomainPartition const & domain ) const
 {
   globalIndex numStick, numNewSlip, numSlip, numOpen;

src/coreComponents/physicsSolvers/contact/SolidMechanicsLagrangeContact.hpp

@@ -138,7 +138,8 @@ class SolidMechanicsLagrangeContact : public ContactSolverBase
 
   bool resetConfigurationToDefault( DomainPartition & domain ) const override final;
 
-  bool updateConfiguration( DomainPartition & domain ) override final;
+  bool updateConfiguration( DomainPartition & domain,
+                            integer configurationLoopIter ) override final;
 
   bool isFractureAllInStickCondition( DomainPartition const & domain ) const;
 
@@ -207,8 +208,29 @@ class SolidMechanicsLagrangeContact : public ContactSolverBase
     constexpr static char const * transMultiplierString() { return "penaltyStiffnessTransMultiplier"; }
 
     constexpr static char const * stabilizationScalingCoefficientString() { return "stabilizationScalingCoefficient"; }
+
+    constexpr static char const * useLocalYieldAccelerationString() { return "useLocalYieldAcceleration"; }
   };
 
+  /// Member variables and functions needed for yield acceleration. Naming convention follows ( Jiang & Tchelepi, 2019 )
+  array1d< real64 > m_x0; // Accelerated variable @ outer iteration v ( two iterations ago )
+  array1d< real64 > m_x1; // Accelerated variable @ outer iteration v + 1 ( previous iteration )
+  array1d< real64 > m_x1_tilde; // Unaccelerated variable @ outer iteration v + 1 ( previous iteration )
+  array1d< real64 > m_x2; // Accelerated variable @ outer iteration v + 2 ( current iteration )
+  array1d< real64 > m_x2_tilde; // Unaccelerated variable @ outer iteration v + 1 ( current iteration )
+  array1d< real64 > m_omega0; // Old relaxation factor
+  array1d< real64 > m_omega1; // New relaxation factor
+  integer m_useLocalYieldAcceleration; // flag for applying acceleration to yield
+
+  void initializeAccelerationVariables( DomainPartition & domain );
+
+  void tryLocalYieldAcceleration( integer configurationLoopIter,
+                                  localIndex kfe,
+                                  real64 currentTau_unscaled,
+                                  real64 limitTau,
+                                  real64 currentTau,
+                                  integer & fractureState );
+
 };
 
 } /* namespace geos */

src/coreComponents/physicsSolvers/multiphysics/CoupledSolver.hpp

@@ -359,12 +359,13 @@ class CoupledSolver : public PhysicsSolverBase
     return isConverged;
   }
 
-  virtual bool updateConfiguration( DomainPartition & domain ) override
+  virtual bool updateConfiguration( DomainPartition & domain,
+                                    integer const configurationLoopIter ) override
   {
     bool result = true;
     forEachArgInTuple( m_solvers, [&]( auto & solver, auto )
     {
-      result &= solver->updateConfiguration( domain );
+      result &= solver->updateConfiguration( domain, configurationLoopIter );
     } );
     return result;
   }
@@ -689,7 +690,9 @@ class CoupledSolver : public PhysicsSolverBase
     }
 
     if( m_nonlinearSolverParameters.m_nonlinearAccelerationType != NonlinearSolverParameters::NonlinearAccelerationType::None )
+    {
       validateNonlinearAcceleration();
+    }
   }
 
   virtual void validateNonlinearAcceleration()