diff --git a/src/coreComponents/physicsSolvers/PhysicsSolverBase.cpp b/src/coreComponents/physicsSolvers/PhysicsSolverBase.cpp index e0eefbde8c..e2bc95c717 100644 --- a/src/coreComponents/physicsSolvers/PhysicsSolverBase.cpp +++ b/src/coreComponents/physicsSolvers/PhysicsSolverBase.cpp @@ -802,7 +802,7 @@ real64 PhysicsSolverBase::nonlinearImplicitStep( real64 const & time_n, if( isNewtonConverged ) { - isConfigurationLoopConverged = updateConfiguration( domain ); + isConfigurationLoopConverged = updateConfiguration( domain, configurationLoopIter ); if( isConfigurationLoopConverged ) { @@ -1334,7 +1334,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; } diff --git a/src/coreComponents/physicsSolvers/PhysicsSolverBase.hpp b/src/coreComponents/physicsSolvers/PhysicsSolverBase.hpp index 95a87f1a87..1dc86196a4 100644 --- a/src/coreComponents/physicsSolvers/PhysicsSolverBase.hpp +++ b/src/coreComponents/physicsSolvers/PhysicsSolverBase.hpp @@ -568,9 +568,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 diff --git a/src/coreComponents/physicsSolvers/multiphysics/CoupledSolver.hpp b/src/coreComponents/physicsSolvers/multiphysics/CoupledSolver.hpp index 8517edb3df..97f64113ea 100644 --- a/src/coreComponents/physicsSolvers/multiphysics/CoupledSolver.hpp +++ b/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() diff --git a/src/coreComponents/physicsSolvers/solidMechanics/contact/SolidMechanicsAugmentedLagrangianContact.cpp b/src/coreComponents/physicsSolvers/solidMechanics/contact/SolidMechanicsAugmentedLagrangianContact.cpp index 58f96fbee2..c5d79759cc 100644 --- a/src/coreComponents/physicsSolvers/solidMechanics/contact/SolidMechanicsAugmentedLagrangianContact.cpp +++ b/src/coreComponents/physicsSolvers/solidMechanics/contact/SolidMechanicsAugmentedLagrangianContact.cpp @@ -810,7 +810,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; diff --git a/src/coreComponents/physicsSolvers/solidMechanics/contact/SolidMechanicsAugmentedLagrangianContact.hpp b/src/coreComponents/physicsSolvers/solidMechanics/contact/SolidMechanicsAugmentedLagrangianContact.hpp index 04ad72ab3e..65e934c114 100644 --- a/src/coreComponents/physicsSolvers/solidMechanics/contact/SolidMechanicsAugmentedLagrangianContact.hpp +++ b/src/coreComponents/physicsSolvers/solidMechanics/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; /** diff --git a/src/coreComponents/physicsSolvers/solidMechanics/contact/SolidMechanicsEmbeddedFractures.cpp b/src/coreComponents/physicsSolvers/solidMechanics/contact/SolidMechanicsEmbeddedFractures.cpp index b5f1ac23d6..b7a3798f87 100644 --- a/src/coreComponents/physicsSolvers/solidMechanics/contact/SolidMechanicsEmbeddedFractures.cpp +++ b/src/coreComponents/physicsSolvers/solidMechanics/contact/SolidMechanicsEmbeddedFractures.cpp @@ -799,7 +799,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; diff --git a/src/coreComponents/physicsSolvers/solidMechanics/contact/SolidMechanicsEmbeddedFractures.hpp b/src/coreComponents/physicsSolvers/solidMechanics/contact/SolidMechanicsEmbeddedFractures.hpp index 7bc61b0e52..07be4773ab 100644 --- a/src/coreComponents/physicsSolvers/solidMechanics/contact/SolidMechanicsEmbeddedFractures.hpp +++ b/src/coreComponents/physicsSolvers/solidMechanics/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; } diff --git a/src/coreComponents/physicsSolvers/solidMechanics/contact/SolidMechanicsLagrangeContact.cpp b/src/coreComponents/physicsSolvers/solidMechanics/contact/SolidMechanicsLagrangeContact.cpp index 9d5a7b1f02..9c250ed350 100644 --- a/src/coreComponents/physicsSolvers/solidMechanics/contact/SolidMechanicsLagrangeContact.cpp +++ b/src/coreComponents/physicsSolvers/solidMechanics/contact/SolidMechanicsLagrangeContact.cpp @@ -71,6 +71,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 >(); } @@ -187,6 +192,11 @@ void SolidMechanicsLagrangeContact::initializePreSubGroups() } ); } + if( m_useLocalYieldAcceleration ) + { + GEOS_LOG_LEVEL_RANK_0( 1, GEOS_FMT( "{}: local yield acceleration enabled", getName() ) ); + } + } void SolidMechanicsLagrangeContact::setupSystem( DomainPartition & domain, @@ -208,6 +218,11 @@ void SolidMechanicsLagrangeContact::setupSystem( DomainPartition & domain, { createPreconditioner( domain ); } + + if( m_useLocalYieldAcceleration ) + { + initializeAccelerationVariables( domain ); + } } void SolidMechanicsLagrangeContact::implicitStepSetup( real64 const & time_n, @@ -2218,7 +2233,8 @@ bool SolidMechanicsLagrangeContact::resetConfigurationToDefault( DomainPartition return false; } -bool SolidMechanicsLagrangeContact::updateConfiguration( DomainPartition & domain ) +bool SolidMechanicsLagrangeContact::updateConfiguration( DomainPartition & domain, + integer const configurationLoopIter ) { GEOS_MARK_FUNCTION; @@ -2277,6 +2293,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], @@ -2291,6 +2308,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); @@ -2313,6 +2333,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 = {}", @@ -2345,6 +2370,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; diff --git a/src/coreComponents/physicsSolvers/solidMechanics/contact/SolidMechanicsLagrangeContact.hpp b/src/coreComponents/physicsSolvers/solidMechanics/contact/SolidMechanicsLagrangeContact.hpp index 15cefe4e14..783d8ee3c4 100644 --- a/src/coreComponents/physicsSolvers/solidMechanics/contact/SolidMechanicsLagrangeContact.hpp +++ b/src/coreComponents/physicsSolvers/solidMechanics/contact/SolidMechanicsLagrangeContact.hpp @@ -133,7 +133,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; @@ -206,8 +207,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 */