diff --git a/inputFiles/constitutiveDriver/friction/frictionDriver_Coulomb.xml b/inputFiles/constitutiveDriver/friction/frictionDriver_Coulomb.xml new file mode 100644 index 00000000000..70de7cc8d4c --- /dev/null +++ b/inputFiles/constitutiveDriver/friction/frictionDriver_Coulomb.xml @@ -0,0 +1,25 @@ + + + + + + + + + + + + + diff --git a/inputFiles/constitutiveDriver/friction/frictionDriver_base.xml b/inputFiles/constitutiveDriver/friction/frictionDriver_base.xml new file mode 100644 index 00000000000..68b3b0e7148 --- /dev/null +++ b/inputFiles/constitutiveDriver/friction/frictionDriver_base.xml @@ -0,0 +1,58 @@ + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + diff --git a/inputFiles/constitutiveDriver/friction/tables/djumps.geos b/inputFiles/constitutiveDriver/friction/tables/djumps.geos new file mode 100644 index 00000000000..e1174543015 --- /dev/null +++ b/inputFiles/constitutiveDriver/friction/tables/djumps.geos @@ -0,0 +1,6 @@ +0 +1e-5 +2e-5 +3e-5 +4e-5 +5e-5 \ No newline at end of file diff --git a/inputFiles/constitutiveDriver/friction/tables/jumps.geos b/inputFiles/constitutiveDriver/friction/tables/jumps.geos new file mode 100644 index 00000000000..3c922732f7a --- /dev/null +++ b/inputFiles/constitutiveDriver/friction/tables/jumps.geos @@ -0,0 +1,6 @@ +0 +1e-3 +2e-3 +3e-3 +4e-3 +5e-3 diff --git a/inputFiles/constitutiveDriver/friction/tables/time.geos b/inputFiles/constitutiveDriver/friction/tables/time.geos new file mode 100644 index 00000000000..e8371f00609 --- /dev/null +++ b/inputFiles/constitutiveDriver/friction/tables/time.geos @@ -0,0 +1,6 @@ +0 +1 +2 +3 +4 +5 diff --git a/inputFiles/constitutiveDriver/friction/tables/tractions.geos b/inputFiles/constitutiveDriver/friction/tables/tractions.geos new file mode 100644 index 00000000000..1655949d91f --- /dev/null +++ b/inputFiles/constitutiveDriver/friction/tables/tractions.geos @@ -0,0 +1,6 @@ +0 +-1e6 +-2e6 +-3e6 +-4e6 +-5e6 diff --git a/src/coreComponents/constitutive/contact/CoulombFriction.hpp b/src/coreComponents/constitutive/contact/CoulombFriction.hpp index e5a3387743d..5bf8057fef6 100644 --- a/src/coreComponents/constitutive/contact/CoulombFriction.hpp +++ b/src/coreComponents/constitutive/contact/CoulombFriction.hpp @@ -182,6 +182,14 @@ class CoulombFriction : public FrictionBase */ KernelWrapper createKernelUpdates() const; + /// getting cohesion value + real64 getCohesion() const + { return m_cohesion; } + + /// getting friction coeff + real64 getFrictionCoeff() const + { return m_frictionCoefficient; } + /** * @struct Set of "char const *" and keys for data specified in this class. */ diff --git a/src/coreComponents/constitutiveDrivers/CMakeLists.txt b/src/coreComponents/constitutiveDrivers/CMakeLists.txt index e3b7a38416c..3b6268fe2e7 100644 --- a/src/coreComponents/constitutiveDrivers/CMakeLists.txt +++ b/src/coreComponents/constitutiveDrivers/CMakeLists.txt @@ -29,12 +29,16 @@ set( constitutiveDrivers_headers relativePermeability/RelpermDriver.hpp relativePermeability/RelpermDriverRunTest.hpp solid/TriaxialDriver.hpp + contact/FrictionDriver.hpp + contact/FrictionDriverRunTest.hpp ) # # Specify all sources # set( constitutiveDrivers_sources ConstitutiveDriver.cpp + contact/FrictionDriver.cpp + contact/FrictionDriverRunTest.cpp fluid/multiFluid/PVTDriver.cpp fluid/multiFluid/constant/PVTDriverRunTestInvariantImmiscibleFluid.cpp fluid/multiFluid/blackOil/PVTDriverRunTestDeadOilFluid.cpp diff --git a/src/coreComponents/constitutiveDrivers/contact/FrictionDriver.cpp b/src/coreComponents/constitutiveDrivers/contact/FrictionDriver.cpp new file mode 100644 index 00000000000..a57b775f1d9 --- /dev/null +++ b/src/coreComponents/constitutiveDrivers/contact/FrictionDriver.cpp @@ -0,0 +1,247 @@ +/* + * ------------------------------------------------------------------------------------------------------------ + * SPDX-License-Identifier: LGPL-2.1-only + * + * Copyright (c) 2016-2024 Lawrence Livermore National Security LLC + * Copyright (c) 2018-2024 TotalEnergies + * Copyright (c) 2018-2024 The Board of Trustees of the Leland Stanford Junior University + * Copyright (c) 2023-2024 Chevron + * Copyright (c) 2019- GEOS/GEOSX Contributors + * All rights reserved + * + * See top level LICENSE, COPYRIGHT, CONTRIBUTORS, NOTICE, and ACKNOWLEDGEMENTS files for details. + * ------------------------------------------------------------------------------------------------------------ + */ + +#include "FrictionDriver.hpp" + +#include "constitutive/ConstitutiveManager.hpp" +#include "constitutiveDrivers/LogLevelsInfo.hpp" +#include "constitutive/contact/FrictionBase.hpp" +#include "constitutive/contact/FrictionSelector.hpp" + +#include "functions/FunctionManager.hpp" +#include "functions/TableFunction.hpp" + +namespace geos +{ + +using namespace dataRepository; +using namespace constitutive; + +FrictionDriver::FrictionDriver( const string & name, Group * const parent ) + : ConstitutiveDriver( name, parent ) +{ + registerWrapper( viewKeyStruct::frictionNameString(), &m_frictionName ). + setRTTypeName( rtTypes::CustomTypes::groupNameRef ). + setInputFlag( InputFlags::REQUIRED ). + setDescription( "Friction model to test" ); + + registerWrapper( viewKeyStruct::numStepsString(), &m_numSteps ). + setInputFlag( InputFlags::REQUIRED ). + setDescription( "Number of sample step to take in both jumps and traction increments" ); + + registerWrapper( viewKeyStruct::jumpFunctionString(), &m_jumpFunctionName ). + setInputFlag( InputFlags::REQUIRED ). + setDescription( "Name of the input function representing jump function along world x-axis" ); + + registerWrapper( viewKeyStruct::dJumpFunctionString(), &m_dJumpFunctionName ). + setInputFlag( InputFlags::REQUIRED ). + setDescription( "Name of the input function representing deltaDisplacementJump function along world x-axis" ); + + registerWrapper( viewKeyStruct::tractionFunctionString(), &m_tractionFunctionName ). + setInputFlag( InputFlags::REQUIRED ). + setDescription( "Name of the input function representing traction function along world x-axis" ); + + registerWrapper( viewKeyStruct::thetaString(), &m_theta ). + setInputFlag( InputFlags::REQUIRED ). + setDescription( "x-Tilt angle in degree" ); + + registerWrapper( viewKeyStruct::phiString(), &m_phi ). + setInputFlag( InputFlags::INVALID ). + setDescription( "y-Tilt angle in degree" ); + + //first batch of parameters + registerWrapper( viewKeyStruct::normalDispTol(), &m_normalDispTol ). + setInputFlag( InputFlags::REQUIRED ). + setDescription( "normal Displacement Tolerance (scale as inverse of average Young modulus)." ); + + registerWrapper( viewKeyStruct::normalTractionTol(), &m_normalTracTol ). + setInputFlag( InputFlags::REQUIRED ). + setDescription( "normal Traction Tolerance" ); + + registerWrapper( viewKeyStruct::slidingTol(), &m_slidingTol ). + setInputFlag( InputFlags::REQUIRED ). + setDescription( "tangential Displacement Tolerance" ); + + registerWrapper( viewKeyStruct::simultaneous(), &m_isSimultaneous ). + setInputFlag( InputFlags::REQUIRED ). + setDescription( "isSimultaneous" ); + + + + addLogLevel< logInfo::LogOutput >(); +} + +void FrictionDriver::postInputInitialization() +{ + ConstitutiveDriver::postInputInitialization(); + + // Check that the functions exist + FunctionManager & functionManager = FunctionManager::getInstance(); + GEOS_ERROR_IF( !functionManager.hasGroup< TableFunction >( m_jumpFunctionName ), + GEOS_FMT( "Jump function with name '{}' not found", m_jumpFunctionName ), + getWrapperDataContext( viewKeyStruct::jumpFunctionString() ) ); + + GEOS_ERROR_IF( !functionManager.hasGroup< TableFunction >( m_dJumpFunctionName ), + GEOS_FMT( "dJump function with name '{}' not found", m_dJumpFunctionName ), + getWrapperDataContext( viewKeyStruct::dJumpFunctionString() ) ); + + GEOS_ERROR_IF( !functionManager.hasGroup< TableFunction >( m_tractionFunctionName ), + GEOS_FMT( "Traction function with name '{}' not found", m_tractionFunctionName ), + getWrapperDataContext( viewKeyStruct::tractionFunctionString() ) ); + + string_array columnNames; + getColumnNames( columnNames ); + integer const numCols = static_cast< integer >(columnNames.size()); + + // initialize functions + TableFunction & jumpFunction = functionManager.getGroup< TableFunction >( m_jumpFunctionName ); + TableFunction & dJumpFunction = functionManager.getGroup< TableFunction >( m_dJumpFunctionName ); + TableFunction & tractionFunction = functionManager.getGroup< TableFunction >( m_tractionFunctionName ); + + jumpFunction.initializeFunction(); + dJumpFunction.initializeFunction(); + tractionFunction.initializeFunction(); + + // TODO: Maybe we should take the maximum extent of jumpFunction and tractionFunction + ArrayOfArraysView< real64 > coordinates = jumpFunction.getCoordinates(); + real64 const minTime = coordinates[0][0]; + real64 const maxTime = coordinates[0][coordinates.sizeOfArray( 0 )-1]; + + // Allocate the data + allocateTable( numCols, minTime, maxTime ); + + // set input columns + initializeTable(); +} + +bool FrictionDriver::execute() +{ + FrictionBase & baseFriction = getFriction(); + + GEOS_LOG_LEVEL_RANK_0( logInfo::LogOutput, "Launching Friction Driver" ); + GEOS_LOG_LEVEL_RANK_0( logInfo::LogOutput, " Friction ............... " << m_frictionName ); + GEOS_LOG_LEVEL_RANK_0( logInfo::LogOutput, " Type ................... " << baseFriction.getCatalogName() ); + GEOS_LOG_LEVEL_RANK_0( logInfo::LogOutput, " Steps .................. " << m_numSteps ); + GEOS_LOG_LEVEL_RANK_0( logInfo::LogOutput, " Output ................. " << m_outputFile ); + + // create a dummy discretization with one quadrature point for + // storing constitutive data + conduit::Node node; + dataRepository::Group rootGroup( "root", node ); + dataRepository::Group discretization( "discretization", &rootGroup ); + + integer const numRows = m_table.size( 0 ); + discretization.resize( numRows ); // numRows elements + baseFriction.allocateConstitutiveData( discretization, 1 ); // one quadrature point + + constitutiveUpdatePassThru( baseFriction, [&]( auto & selectedFrictionModel ) + { + using FRICTION_TYPE = TYPEOFREF( selectedFrictionModel ); + runTest< FRICTION_TYPE >( selectedFrictionModel, m_table ); + } ); + + return false; +} + +void FrictionDriver::getColumnNames( string_array & columnNames ) const +{ + columnNames.emplace_back( "time" ); + columnNames.emplace_back( "traction,normal" ); + columnNames.emplace_back( "traction,tangent1" ); + columnNames.emplace_back( "traction,tangent2" ); + columnNames.emplace_back( "displacement jump,normal" ); + columnNames.emplace_back( "displacement jump,tangent1" ); + columnNames.emplace_back( "displacement jump,tangent2" ); + columnNames.emplace_back( "delta displacement jump,normal" ); + columnNames.emplace_back( "delta displacement jump,tangent1" ); + columnNames.emplace_back( "delta displacement jump,tangent2" ); + columnNames.emplace_back( "encoded constaint (0:converged, 1:stick & gn>0 (opening), 2: interpenetration, 3: stick & gt>lim (disp-sliding), 4: tau>taulim (trac-sliding) )" ); + columnNames.emplace_back( "fracture state (0:stick, 1:slip , 2: new slip, 3: open)" ); + columnNames.emplace_back( "newtraction,normal" ); + columnNames.emplace_back( "newtraction,tangent1" ); + columnNames.emplace_back( "newtraction,tangent2" ); + + if( dynamic_cast< CoulombFriction const * >(&getFriction()) != nullptr ) + { + columnNames.emplace_back( "tau limit" ); + } +} + +void FrictionDriver::initializeTable() +{ + integer const numRows = m_table.size( 0 ); + + FunctionManager & functionManager = FunctionManager::getInstance(); + TableFunction const & jumpFunction = functionManager.getGroup< TableFunction >( m_jumpFunctionName ); + TableFunction const & dJumpFunction = functionManager.getGroup< TableFunction >( m_dJumpFunctionName ); + TableFunction const & tractionFunction = functionManager.getGroup< TableFunction >( m_tractionFunctionName ); + + real64 const cos_theta = cos( m_theta * M_PI/180.0 ); + real64 const sin_theta = sin( m_theta * M_PI/180.0 ); + + real64 const cos_phi = cos( m_phi * M_PI/180.0 ); + real64 const sin_phi = sin( m_phi * M_PI/180.0 ); + + for( integer index = 0; index < numRows; ++index ) + { + real64 const time = m_table( index, TIME ); + + real64 const traction = tractionFunction.evaluate( &time ); + m_table( index, NTRAC ) = traction*sin_theta; + m_table( index, STRAC0 ) = traction*cos_theta*cos_phi; + m_table( index, STRAC1 ) = traction*cos_theta*sin_phi; + + real64 const jump = jumpFunction.evaluate( &time ); + m_table( index, NJUMP ) = jump*sin_theta; + m_table( index, SLIP0 ) = jump*cos_theta*cos_phi; + m_table( index, SLIP1 ) = jump*cos_theta*sin_phi; + + real64 const dJump = dJumpFunction.evaluate( &time ); + m_table( index, NDJUMP ) = dJump*sin_theta; + m_table( index, DSLIP0 ) = dJump*cos_theta*cos_phi; + m_table( index, DSLIP1 ) = dJump*cos_theta*sin_phi; + + m_table( index, CC ) = 0; + m_table( index, FS ) = fields::contact::FractureState::Stick; + } + + if( CoulombFriction const * coulombFriction = dynamic_cast< CoulombFriction const * >(&getFriction()) ) + { + real64 const cohesion = coulombFriction->getCohesion(); + real64 const frictionCoeff = coulombFriction->getFrictionCoeff(); + for( integer index = 0; index < numRows; ++index ) + { + real64 const normal_traction = m_table( index, NTRAC ); + m_table( index, TLIM ) = cohesion - normal_traction * frictionCoeff; + } + } +} + +FrictionBase & FrictionDriver::getFriction() +{ + return getConstitutiveManager().getGroup< FrictionBase >( m_frictionName ); +} + +FrictionBase const & FrictionDriver::getFriction() const +{ + return getConstitutiveManager().getGroup< FrictionBase >( m_frictionName ); +} + + +REGISTER_CATALOG_ENTRY( TaskBase, + FrictionDriver, + string const &, dataRepository::Group * const ) + +} diff --git a/src/coreComponents/constitutiveDrivers/contact/FrictionDriver.hpp b/src/coreComponents/constitutiveDrivers/contact/FrictionDriver.hpp new file mode 100644 index 00000000000..771cb411d36 --- /dev/null +++ b/src/coreComponents/constitutiveDrivers/contact/FrictionDriver.hpp @@ -0,0 +1,121 @@ +/* + * ------------------------------------------------------------------------------------------------------------ + * SPDX-License-Identifier: LGPL-2.1-only + * + * Copyright (c) 2016-2024 Lawrence Livermore National Security LLC + * Copyright (c) 2018-2024 TotalEnergies + * Copyright (c) 2018-2024 The Board of Trustees of the Leland Stanford Junior University + * Copyright (c) 2023-2024 Chevron + * Copyright (c) 2019- GEOS/GEOSX Contributors + * All rights reserved + * + * See top level LICENSE, COPYRIGHT, CONTRIBUTORS, NOTICE, and ACKNOWLEDGEMENTS files for details. + * ------------------------------------------------------------------------------------------------------------ + */ + +#ifndef GEOS_CONSTITUTIVEDRIVERS_CONTACT_FRICTIONDRIVER_HPP +#define GEOS_CONSTITUTIVEDRIVERS_CONTACT_FRICTIONDRIVER_HPP + +#include "constitutiveDrivers/ConstitutiveDriver.hpp" +#include "physicsSolvers/solidMechanics/contact/ContactSolverBase.hpp" +// #include "physicsSolvers/solidMechanics/contact/SolidMechanicsAugmentedLagrangianContact.hpp" + +namespace geos +{ +namespace constitutive +{ +class FrictionBase; +} + +class FrictionDriver : public ConstitutiveDriver +{ +public: + FrictionDriver( const string & name, + Group * const parent ); + + static string catalogName() + { return "FrictionDriver"; } + + void postInputInitialization() override; + + bool execute() override; + + void getColumnNames( string_array & columnNames ) const override; + + template< typename FRICTION_TYPE > + void + runTest( FRICTION_TYPE & friction, + const arrayView2d< real64, 1 > & table ); + +private: + /** + * @brief Get the friction model from the catalog + */ + constitutive::FrictionBase & getFriction(); + constitutive::FrictionBase const & getFriction() const; + + void initializeTable(); + + /** + * @struct viewKeyStruct holds char strings and viewKeys for fast lookup + */ + struct viewKeyStruct : ConstitutiveDriver::viewKeyStruct + { + constexpr static char const * frictionNameString() + { return "friction"; } + + constexpr static char const * contactNameString() + { return "contact"; } + + constexpr static char const * jumpFunctionString() + { return "jumpControl"; } + + constexpr static char const * dJumpFunctionString() + { return "dJumpControl"; } + + constexpr static char const * tractionFunctionString() + { return "tractionControl"; } + + constexpr static char const * thetaString() + { return "xTiltAngle";} + + constexpr static char const * phiString() + { return "yTiltAngle";} + + constexpr static char const * normalDispTol() + { return "tolJumpN"; } + + constexpr static char const * normalTractionTol() + { return "tolNormalTrac"; } + + constexpr static char const * slidingTol() + { return "tolJumpT"; } + + constexpr static char const * simultaneous() + { return "simultaneous"; } + + + }; + + // Time is defined in base class + enum columnKeys { NTRAC=1, STRAC0, STRAC1, NJUMP, SLIP0, SLIP1, NDJUMP, DSLIP0, DSLIP1, CC, FS, NEWTRAC, SNEWTRAC0, SNEWTRAC1, TLIM }; + + string m_jumpFunctionName; ///< + string m_dJumpFunctionName; ///< + string m_tractionFunctionName; ///< + + real64 m_theta{0.0}; ///< x-tilt of fault + real64 m_phi{0.0}; ///< y-tilt of fault + + real64 m_normalDispTol{1.e-8}; + real64 m_normalTracTol{100.}; + real64 m_slidingTol{1e-5}; + + integer m_isSimultaneous{1}; + + string m_frictionName; ///< frictionType identifier +}; + +} + +#endif //GEOS_CONSTITUTIVEDRIVERS_CONTACT_FRICTIONDRIVER_HPP diff --git a/src/coreComponents/constitutiveDrivers/contact/FrictionDriverRunTest.cpp b/src/coreComponents/constitutiveDrivers/contact/FrictionDriverRunTest.cpp new file mode 100644 index 00000000000..f06fe66436d --- /dev/null +++ b/src/coreComponents/constitutiveDrivers/contact/FrictionDriverRunTest.cpp @@ -0,0 +1,27 @@ +#include "FrictionDriverRunTest.hpp" +#include "constitutive/contact/CoulombFriction.hpp" +#include "constitutive/contact/FrictionlessContact.hpp" +#include "constitutive/contact/RateAndStateFriction.hpp" +#include + + +namespace geos +{ + +template +void +FrictionDriver::runTest( constitutive::CoulombFriction &, const arrayView2d< real64 > & ); + +template +void +FrictionDriver::runTest( constitutive::FrictionlessContact &, const arrayView2d< real64 > & ); + +template +void +FrictionDriver::runTest( constitutive::RateAndStateFriction< std::integral_constant< bool, true > > &, const arrayView2d< real64 > & ); + +template +void +FrictionDriver::runTest( constitutive::RateAndStateFriction< std::integral_constant< bool, false > > &, const arrayView2d< real64 > & ); + +} diff --git a/src/coreComponents/constitutiveDrivers/contact/FrictionDriverRunTest.hpp b/src/coreComponents/constitutiveDrivers/contact/FrictionDriverRunTest.hpp new file mode 100644 index 00000000000..9ba716392ef --- /dev/null +++ b/src/coreComponents/constitutiveDrivers/contact/FrictionDriverRunTest.hpp @@ -0,0 +1,117 @@ +/* + * ------------------------------------------------------------------------------------------------------------ + * SPDX-License-Identifier: LGPL-2.1-only + * + * Copyright (c) 2016-2024 Lawrence Livermore National Security LLC* + * Copyright (c) 2018-2024 TotalEnergies + * Copyright (c) 2018-2024 The Board of Trustees of the Leland Stanford Junior University + * Copyright (c) 2023-2024 Chevron + * Copyright (c) 2019- GEOS/GEOSX Contributors + * All rights reserved + * + * See top level LICENSE, COPYRIGHT, CONTRIBUTORS, NOTICE, and ACKNOWLEDGEMENTS files for details. + * ------------------------------------------------------------------------------------------------------------ + */ + +#ifndef GEOS_FRICTIONDRIVERRUNTEST_HPP_ +#define GEOS_FRICTIONDRIVERRUNTEST_HPP_ + +#include "constitutiveDrivers/contact/FrictionDriver.hpp" +#include "physicsSolvers/solidMechanics/contact/FractureState.hpp" +#include "physicsSolvers/solidMechanics/contact/SolidMechanicsAugmentedLagrangianContact.hpp" +#include "constitutive/solid/SolidFields.hpp" + +#include + +namespace geos +{ + +template< typename FRICTION_TYPE > +void +FrictionDriver::runTest( FRICTION_TYPE & friction, + const arrayView2d< real64 > & table ) +{ + + array1d< integer > const ghostRank( 1 ); ghostRank[0] = -1; + array1d< real64 > const normalDisplacementTol( 1 ); normalDisplacementTol[0]=m_normalDispTol;//normalDispTol should scale as 1/E + array1d< real64 > const normalTractionTol( 1 ); normalTractionTol[0]=m_normalTracTol; + array1d< real64 > const slidingTol( 1 ); slidingTol[0]=m_slidingTol; + + real64 kt = 1e7; + array2d< real64 > const iterPen( 1, 5 ); + iterPen[0][0] = 100*kt; + iterPen[0][1] = kt; + iterPen[0][2] = kt; iterPen[0][3] = kt; iterPen[0][4] = 0.; + array1d< integer > fractureState( 1 ); + + fractureState[0] = fields::contact::FractureState::Stick; + array2d< real64 > traction( 1, 3 ); + array2d< real64 > jump( 1, 3 ); + array2d< real64 > djump( 1, 3 ); + + bool isSimultaneous = m_isSimultaneous; + real64 const slidingCheckTol = .05; //default + + //TODO computeTolerance eleme to Elem + typename FRICTION_TYPE::KernelWrapper const kernelWrapper = friction.createKernelUpdates(); + + integer const numRows = m_table.size( 0 ); + forAll< parallelDevicePolicy<> >( numRows, + [&friction, &table, &kernelWrapper, + &ghostRank, &normalDisplacementTol, &normalTractionTol, &slidingTol, &iterPen, + &slidingCheckTol, &isSimultaneous, + &jump, &djump, + &fractureState, &traction ] + GEOS_HOST_DEVICE ( integer const ei ) + { + + GEOS_LOG_RANK( "[debug] Table Evaluation" ); + + jump[0][0] = table( ei, NJUMP ); + jump[0][1] = table( ei, SLIP0 ); + jump[0][2] = table( ei, SLIP1 ); + + djump[0][0] = table( ei, NDJUMP ); + djump[0][1] = table( ei, DSLIP0 ); + djump[0][2] = table( ei, DSLIP1 ); + + + traction[0][0] = table( ei, NTRAC ); + traction[0][1] = table( ei, STRAC0 ); + traction[0][2] = table( ei, STRAC1 ); + + kernelWrapper.updateFractureState( jump[0], + traction[0], + fractureState[0] ); + + auto [newTraction, condCov] = SolidMechanicsAugmentedLagrangianContact::updateTractionAndConstraintCheck( 1, + friction, + isSimultaneous, + slidingCheckTol, + normalDisplacementTol, + normalTractionTol, + slidingTol, + iterPen, + jump, + djump, + ghostRank, + fractureState.toView(), + traction.toView() + ); + + kernelWrapper.updateFractureState( jump[0], + newTraction[0], + fractureState[0] ); + + table( ei, CC ) = condCov[0]; + table( ei, FS ) = fractureState[0]; + table( ei, NEWTRAC ) = newTraction[0][0]; + table( ei, SNEWTRAC0 ) = newTraction[0][1]; + table( ei, SNEWTRAC1 ) = newTraction[0][2]; + + } ); +} + +} + +#endif //GEOS_FRICTIONDRIVERRUNTEST_HPP_ diff --git a/src/coreComponents/physicsSolvers/solidMechanics/contact/SolidMechanicsAugmentedLagrangianContact.cpp b/src/coreComponents/physicsSolvers/solidMechanics/contact/SolidMechanicsAugmentedLagrangianContact.cpp index 8dd1d59330e..d3717daadc1 100644 --- a/src/coreComponents/physicsSolvers/solidMechanics/contact/SolidMechanicsAugmentedLagrangianContact.cpp +++ b/src/coreComponents/physicsSolvers/solidMechanics/contact/SolidMechanicsAugmentedLagrangianContact.cpp @@ -40,7 +40,9 @@ #include "finiteElement/FiniteElementDiscretization.hpp" #include "mesh/DomainPartition.hpp" +#include #include +#include #if defined( GEOS_USE_CUDA ) #include @@ -1067,6 +1069,84 @@ void SolidMechanicsAugmentedLagrangianContact::updateState( DomainPartition & do GEOS_UNUSED_VAR( domain ); } + +std::tuple< array2d< real64 >, array1d< int > > SolidMechanicsAugmentedLagrangianContact::updateTractionAndConstraintCheck( + std::ptrdiff_t const rsize, + FrictionBase const & frictionLaw, + bool isSimultaneous, + real64 const slidingCheckTolerance, + arrayView1d< real64 const > const & normalDisplacementTolerance, + arrayView1d< real64 const > const & normalTractionTolerance, + arrayView1d< real64 const > const & slidingTolerance, + arrayView2d< real64 const > const & iterativePenalty, + arrayView2d< real64 const > const & dispJump, + arrayView2d< real64 const > const & deltaDispJump, + arrayView1d< integer const > const & ghostRank, + arrayView1d< integer const > const & fractureState, + arrayView2d< real64 > const & traction ) +{ + array2d< real64 > traction_new; + std::ptrdiff_t const sizes[2] = {rsize, 3}; + traction_new.resize( 2, sizes ); + array1d< int > condConv; + condConv.resize( rsize ); + + // Update the traction field based on the displacement results from the nonlinear solve + constitutiveUpdatePassThru( frictionLaw, [&] ( auto & castedFrictionLaw ) + { + using FrictionType = TYPEOFREF( castedFrictionLaw ); + typename FrictionType::KernelWrapper frictionWrapper = castedFrictionLaw.createKernelUpdates(); + + if( isSimultaneous ) + { + solidMechanicsALMKernels::ComputeTractionSimultaneousKernel:: + launch< parallelDevicePolicy<> >( rsize, + iterativePenalty, + traction, + dispJump, + deltaDispJump, + traction_new.toView() ); + } + else + { + solidMechanicsALMKernels::ComputeTractionKernel:: + launch< parallelDevicePolicy<> >( rsize, + frictionWrapper, + iterativePenalty, + traction, + dispJump, + deltaDispJump, + traction_new.toView() ); + } + } ); + + // real64 const slidingCheckTolerance = m_slidingCheckTolerance; + + constitutiveUpdatePassThru( frictionLaw, [&] ( auto & castedFrictionLaw ) + { + using FrictionType = TYPEOFREF( castedFrictionLaw ); + typename FrictionType::KernelWrapper frictionWrapper = castedFrictionLaw.createKernelUpdates(); + + solidMechanicsALMKernels::ConstraintCheckKernel:: + launch< parallelDevicePolicy<> >( rsize, + frictionWrapper, + ghostRank, + traction, + dispJump, + deltaDispJump, + normalTractionTolerance, + normalDisplacementTolerance, + slidingTolerance, + slidingCheckTolerance, + fractureState, + condConv.toView() ); + } ); + + + return std::make_tuple( traction_new, condConv ); +} + + bool SolidMechanicsAugmentedLagrangianContact::updateConfiguration( DomainPartition & domain, integer const GEOS_UNUSED_PARAM( configurationLoopIter ) ) { @@ -1107,61 +1187,27 @@ bool SolidMechanicsAugmentedLagrangianContact::updateConfiguration( DomainPartit std::ptrdiff_t const sizes[ 2 ] = {subRegion.size(), 3}; traction_new.resize( 2, sizes ); - arrayView2d< real64 > const traction_new_v = traction_new.toView(); - - condConv.resize( subRegion.size()); + // arrayView2d< real64 > const traction_new_v = traction_new.toView(); + + // condConv.resize( subRegion.size()); + + std::tie( traction_new, condConv ) = updateTractionAndConstraintCheck( + subRegion.size(), + frictionLaw, + m_simultaneous, + m_slidingCheckTolerance, + normalDisplacementTolerance, + normalTractionTolerance, + slidingTolerance, + iterativePenalty, + dispJump, + deltaDispJump, + ghostRank, + fractureState, + traction + ); arrayView1d< int > const condConv_v = condConv.toView(); - // Update the traction field based on the displacement results from the nonlinear solve - constitutiveUpdatePassThru( frictionLaw, [&] ( auto & castedFrictionLaw ) - { - using FrictionType = TYPEOFREF( castedFrictionLaw ); - typename FrictionType::KernelWrapper frictionWrapper = castedFrictionLaw.createKernelUpdates(); - - if( m_simultaneous ) - { - solidMechanicsALMKernels::ComputeTractionSimultaneousKernel:: - launch< parallelDevicePolicy<> >( subRegion.size(), - iterativePenalty, - traction, - dispJump, - deltaDispJump, - traction_new_v ); - } - else - { - solidMechanicsALMKernels::ComputeTractionKernel:: - launch< parallelDevicePolicy<> >( subRegion.size(), - frictionWrapper, - iterativePenalty, - traction, - dispJump, - deltaDispJump, - traction_new_v ); - } - } ); - - real64 const slidingCheckTolerance = m_slidingCheckTolerance; - - constitutiveUpdatePassThru( frictionLaw, [&] ( auto & castedFrictionLaw ) - { - using FrictionType = TYPEOFREF( castedFrictionLaw ); - typename FrictionType::KernelWrapper frictionWrapper = castedFrictionLaw.createKernelUpdates(); - - solidMechanicsALMKernels::ConstraintCheckKernel:: - launch< parallelDevicePolicy<> >( subRegion.size(), - frictionWrapper, - ghostRank, - traction, - dispJump, - deltaDispJump, - normalTractionTolerance, - normalDisplacementTolerance, - slidingTolerance, - slidingCheckTolerance, - fractureState, - condConv_v ); - } ); RAJA::ReduceSum< parallelDeviceReduce, localIndex > localSum[5] = { RAJA::ReduceSum< parallelDeviceReduce, localIndex >( 0 ), diff --git a/src/coreComponents/physicsSolvers/solidMechanics/contact/SolidMechanicsAugmentedLagrangianContact.hpp b/src/coreComponents/physicsSolvers/solidMechanics/contact/SolidMechanicsAugmentedLagrangianContact.hpp index 02af0550280..d4e411fdc39 100644 --- a/src/coreComponents/physicsSolvers/solidMechanics/contact/SolidMechanicsAugmentedLagrangianContact.hpp +++ b/src/coreComponents/physicsSolvers/solidMechanics/contact/SolidMechanicsAugmentedLagrangianContact.hpp @@ -22,6 +22,7 @@ #define GEOS_PHYSICSSOLVERS_CONTACT_SOLIDMECHANICSAUGMENTEDLAGRANGIANCONTACT_HPP_ #include "physicsSolvers/solidMechanics/contact/ContactSolverBase.hpp" +#include "constitutive/contact/FrictionSelector.hpp" namespace geos { @@ -215,8 +216,26 @@ class SolidMechanicsAugmentedLagrangianContact : public ContactSolverBase */ void createBubbleCellList( DomainPartition & domain ) const; + + static + std::tuple< array2d< real64 >, array1d< int > > updateTractionAndConstraintCheck( std::ptrdiff_t const rsize, + constitutive::FrictionBase const & frictionLaw, + bool isSimultaneous, + real64 const slidingCheckTolerance, + arrayView1d< real64 const > const & normalDisplacementTolerance, + arrayView1d< real64 const > const & normalTractionTolerance, + arrayView1d< real64 const > const & slidingTolerance, + arrayView2d< real64 const > const & iterativePenalty, + arrayView2d< real64 const > const & dispJump, + arrayView2d< real64 const > const & deltaDispJump, + arrayView1d< integer const > const & ghostRank, + arrayView1d< integer const > const & fractureState, + arrayView2d< real64 > const & traction ); + private: + + /** * @brief Validate that tetrahedral meshes use high-order quadrature rules * @param meshBodies the group containing the mesh bodies diff --git a/src/coreComponents/schema/schema.xsd b/src/coreComponents/schema/schema.xsd index bd776c2bac0..f21af84e624 100644 --- a/src/coreComponents/schema/schema.xsd +++ b/src/coreComponents/schema/schema.xsd @@ -569,6 +569,10 @@ + + + + @@ -6264,6 +6268,7 @@ When set to `all` output both convergence & iteration information to a csv.--> + @@ -6378,6 +6383,29 @@ Information output from lower logLevels is added with the desired log level + + + + + + + + + + + + + + + + + + + + + + + +