From 7c745adaec2ce9679b2893d622415c9aabf3929b Mon Sep 17 00:00:00 2001 From: jacques franc Date: Tue, 4 Nov 2025 18:08:17 +0100 Subject: [PATCH 1/4] starting implementation --- .../contact/FrictionDriver.cpp | 256 ++++++++++++++++++ .../contact/FrictionDriver.hpp | 122 +++++++++ .../contact/FrictionDriverRunTest.hpp | 106 ++++++++ .../contact/LogLevelsInfo.hpp | 52 ++++ 4 files changed, 536 insertions(+) create mode 100644 src/coreComponents/constitutiveDrivers/contact/FrictionDriver.cpp create mode 100644 src/coreComponents/constitutiveDrivers/contact/FrictionDriver.hpp create mode 100644 src/coreComponents/constitutiveDrivers/contact/FrictionDriverRunTest.hpp create mode 100644 src/coreComponents/constitutiveDrivers/contact/LogLevelsInfo.hpp diff --git a/src/coreComponents/constitutiveDrivers/contact/FrictionDriver.cpp b/src/coreComponents/constitutiveDrivers/contact/FrictionDriver.cpp new file mode 100644 index 00000000000..f51425517a9 --- /dev/null +++ b/src/coreComponents/constitutiveDrivers/contact/FrictionDriver.cpp @@ -0,0 +1,256 @@ +/* + * ------------------------------------------------------------------------------------------------------------ + * 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 "common/MpiWrapper.hpp" +#include "constitutive/ConstitutiveManager.hpp" +#include "constitutiveDrivers/contact/LogLevelsInfo.hpp" +#include "constitutive/contact/FrictionBase.hpp" +#include "constitutive/contact/FrictionSelector.hpp" + +#include "FrictionDriver.hpp" + +namespace geos { + +using namespace dataRepository; +using namespace constitutive; + +FrictionDriver::FrictionDriver(const string& name, Group * const parent) +: +TaskBase(name, parent) +{ + registerWrapper( viewKeyStruct::frictionNameString(), &m_frictionName ). + setRTTypeName( rtTypes::CustomTypes::groupNameRef ). + setInputFlag( InputFlags::REQUIRED ). + setDescription( "Relperm model to test" ); + + registerWrapper( viewKeyStruct::numStepsString(), &m_numSteps ). + setInputFlag( InputFlags::REQUIRED ). + setDescription( "Number of saturation steps to take" ); + + registerWrapper( viewKeyStruct::outputString(), &m_outputFile ). + setInputFlag( InputFlags::OPTIONAL ). + setApplyDefaultValue( "none" ). + setDescription( "Output file" ); + + registerWrapper( viewKeyStruct::baselineString(), &m_baselineFile ). + setInputFlag( InputFlags::OPTIONAL ). + setApplyDefaultValue( "none" ). + setDescription( "Baseline file" ); + + addLogLevel< logInfo::LogOutput >(); +} + +void FrictionDriver::outputResults() +{ + // TODO: improve file path output to grab command line -o directory + // for the moment, we just use the specified m_outputFile directly + + FILE * fp = fopen( m_outputFile.c_str(), "w" ); + + fprintf( fp, "# column 1 = time\n" ); + fprintf( fp, "# columns %d-%d = phase vol fractions\n", 2, 1 + m_numPhases ); + fprintf( fp, "# columns %d-%d = phase relperm\n", 2 + m_numPhases, 1 + 2 * m_numPhases ); + + if( ( m_numPhases == 2 && m_table.size( 1 ) > 5 ) || m_table.size( 1 ) > 7 ) + { + fprintf( fp, "# columns %d-%d = phase relperm (hyst)\n", 1 + 2 * m_numPhases, 1 + 3 * m_numPhases ); + } + + + for( integer n = 0; n < m_table.size( 0 ); ++n ) + { + for( integer col = 0; col < m_table.size( 1 ); ++col ) + { + fprintf( fp, "%.4e ", m_table( n, col ) ); + } + fprintf( fp, "\n" ); + } + fclose( fp ); + + +} + + +void FrictionDriver::postInputInitialization() +{ + ConstitutiveManager + & constitutiveManager = this->getGroupByPath< ConstitutiveManager >( "/Problem/domain/Constitutive" ); + FrictionBase& baseFriction = constitutiveManager.getGroup< FrictionBase >( m_frictionName ); + +// m_numPhases = baseFriction.numSubGroups(); + +} + + +bool FrictionDriver::execute( const geos::real64 GEOS_UNUSED_PARAM( time_n ), + const geos::real64 GEOS_UNUSED_PARAM( dt ), + const geos::integer GEOS_UNUSED_PARAM( cycleNumber ), + const geos::integer GEOS_UNUSED_PARAM( eventCounter ), + const geos::real64 GEOS_UNUSED_PARAM( eventProgress ), + geos::DomainPartition & + GEOS_UNUSED_PARAM( domain ) ) +{ + // this code only makes sense in serial + + GEOS_THROW_IF( MpiWrapper::commRank() > 0, "RelpermDriver should only be run in serial", std::runtime_error ); + + + ConstitutiveManager + & constitutiveManager = this->getGroupByPath< ConstitutiveManager >( "/Problem/domain/Constitutive" ); + FrictionBase + & baseFriction = constitutiveManager.getGroup< FrictionBase >( m_frictionName ); + + GEOS_LOG_LEVEL_RANK_0( logInfo::LogOutput, "Launching Relperm Driver" ); + GEOS_LOG_LEVEL_RANK_0( logInfo::LogOutput, " Relperm .................. " << m_frictionName ); + GEOS_LOG_LEVEL_RANK_0( logInfo::LogOutput, " Type ................... " << baseFriction.getCatalogName() ); +// GEOS_LOG_LEVEL_RANK_0( logInfo::LogOutput, " No. of Phases .......... " << m_numPhases ); + GEOS_LOG_LEVEL_RANK_0( logInfo::LogOutput, " Steps .................. " << m_numSteps ); + GEOS_LOG_LEVEL_RANK_0( logInfo::LogOutput, " Output ................. " << m_outputFile ); + GEOS_LOG_LEVEL_RANK_0( logInfo::LogOutput, " Baseline ............... " << m_baselineFile ); + + // 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 ); + + discretization.resize( 1 ); // one element + baseRelperm.allocateConstitutiveData( discretization, 1 ); // one quadrature point + + constitutiveUpdatePassThru( baseFriction, [&]( auto & selectedFrictionModel ) + { + using RELPERM_TYPE = TYPEOFREF( selectedFrictionModel ); + resizeTables< RELPERM_TYPE >(); + runTest< RELPERM_TYPE >( selectedFrictionModel, m_table ); + } ); + + // move table back to host for output + m_table.move( LvArray::MemorySpace::host ); + + if( m_outputFile != "none" ) + { + outputResults(); + } + + if( m_baselineFile != "none" ) + { + compareWithBaseline(); + } + + return false; +} + + +template< typename RELPERM_TYPE > +void FrictionDriver::resizeTables() +{ + ConstitutiveManager + & constitutiveManager = this->getGroupByPath< ConstitutiveManager >( "/Problem/domain/Constitutive" ); + FrictionBase + & baseRelperm = constitutiveManager.getGroup< FrictionBaseBase >( m_frictionName ); + + using PT = RelativePermeabilityBase::PhaseType; + integer const ipWater = baseRelperm.getPhaseOrder()[PT::WATER]; + integer const ipOil = baseRelperm.getPhaseOrder()[PT::OIL]; + integer const ipGas = baseRelperm.getPhaseOrder()[PT::GAS]; + + real64 minSw = 0., minSnw = 0.; + if( baseRelperm.numFluidPhases() > 2 ) + { + minSw = baseRelperm.getWettingPhaseMinVolumeFraction(); + minSnw = baseRelperm.getNonWettingMinVolumeFraction(); + } + else + { + if( ipWater < 0 )// a.k.a o/g + { + minSw = 0; + minSnw = baseRelperm.getNonWettingMinVolumeFraction(); + } + else if( ipGas < 0 || ipOil < 0 )// a.k.a w/o or w/g + { + minSnw = 0; + minSw = baseRelperm.getWettingPhaseMinVolumeFraction(); + } + } + + real64 const dSw = ( 1 - minSw - minSnw ) / m_numSteps; + // set input columns + + resizeTable< RELPERM_TYPE >(); + // 3-phase branch + if( m_numPhases > 2 ) + { + for( integer ni = 0; ni < m_numSteps + 1; ++ni ) + { + for( integer nj = 0; nj < m_numSteps + 1; ++nj ) + { + + integer index = ni * ( m_numSteps + 1 ) + nj; + m_table( index, TIME ) = minSw + index * dSw; + m_table( index, ipWater + 1 ) = minSw + nj * dSw; + m_table( index, ipGas + 1 ) = minSnw + ni * dSw; + m_table( index, ipOil + 1 ) = + 1. - m_table( index, ipWater + 1 ) - m_table( index, ipOil + 1 ); + } + } + } + else // 2-phase branch + { + for( integer ni = 0; ni < m_numSteps + 1; ++ni ) + { + integer index = ni; + m_table( index, TIME ) = minSw + index * dSw; + if( ipWater < 0 ) + { + m_table( index, ipGas + 1 ) = minSnw + ni * dSw; + m_table( index, ipOil + 1 ) = 1. - m_table( index, ipGas + 1 ); + } + else if( ipGas < 0 ) + { + m_table( index, ipWater + 1 ) = minSw + ni * dSw; + m_table( index, ipOil + 1 ) = 1. - m_table( index, ipWater + 1 ); + } + else if( ipOil < 0 ) + { + m_table( index, ipWater + 1 ) = minSw + ni * dSw; + m_table( index, ipGas + 1 ) = 1. - m_table( index, ipWater + 1 ); + } + } + + } + + +} + + +// template< typename RELPERM_TYPE > +// std::enable_if_t< std::is_same< TableRelativePermeabilityHysteresis, RELPERM_TYPE >::value, void > +// RelpermDriver::resizeTable() +// { +// if( m_numPhases > 2 ) +// { +// m_table.resize( ( m_numSteps + 1 ) * ( m_numSteps + 1 ), 1 + 3 * m_numPhases ); +// } +// else +// { +// m_table.resize( m_numSteps + 1, 1 + 3 * m_numPhases ); +// } + +// } + + +} \ No newline at end of file diff --git a/src/coreComponents/constitutiveDrivers/contact/FrictionDriver.hpp b/src/coreComponents/constitutiveDrivers/contact/FrictionDriver.hpp new file mode 100644 index 00000000000..ae02bb12a5c --- /dev/null +++ b/src/coreComponents/constitutiveDrivers/contact/FrictionDriver.hpp @@ -0,0 +1,122 @@ +/* + * ------------------------------------------------------------------------------------------------------------ + * 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_frictionDRIVER_HPP_ +#define GEOS_frictionDRIVER_HPP_ + +#include "events/tasks/TaskBase.hpp" + +namespace geos +{ + +class FrictionDriver : public TaskBase +{ + +public: + FrictionDriver( const string & name, + Group * const parent ); + + static string catalogName() + { return "frictionDriver"; } + + void postInputInitialization() override; + + virtual bool execute( real64 const GEOS_UNUSED_PARAM( time_n ), + real64 const GEOS_UNUSED_PARAM( dt ), + integer const GEOS_UNUSED_PARAM( cycleNumber ), + integer const GEOS_UNUSED_PARAM( eventCounter ), + real64 const GEOS_UNUSED_PARAM( eventProgress ), + DomainPartition & + GEOS_UNUSED_PARAM( domain ) ) override; + + // /** + // * @brief Run test using loading protocol in table + // * @param i friction constitutive model + // * @param table Table with input/output time history + // */ + // template< typename friction_TYPE > + // std::enable_if_t< std::is_same< constitutive::TableRelativePermeabilityHysteresis, friction_TYPE >::value, void > + // runTest( friction_TYPE & friction, + // const arrayView2d< real64, 1 > & table ); + + // template< typename friction_TYPE > + // std::enable_if_t< !std::is_same< constitutive::TableRelativePermeabilityHysteresis, friction_TYPE >::value, void > + // runTest( friction_TYPE & friction, + // const arrayView2d< real64, 1 > & table ); + + /** + * @brief Ouput table to file for easy plotting + */ + void outputResults(); + + /** + * @brief Read in a baseline table from file and compare with computed one (for unit testing purposes) + */ + void compareWithBaseline(); + +private: + + template< typename FRICTION_TYPE > + void resizeTables(); + + // template< typename friction_TYPE > + // std::enable_if_t< std::is_same< constitutive::TableRelativePermeabilityHysteresis, friction_TYPE >::value, void > + // resizeTable(); + + // template< typename friction_TYPE > + // std::enable_if_t< !std::is_same< constitutive::TableRelativePermeabilityHysteresis, friction_TYPE >::value, void > + // resizeTable(); + + /** + * @struct viewKeyStruct holds char strings and viewKeys for fast lookup + */ + struct viewKeyStruct + { + constexpr static char const * frictionNameString() + { return "friction"; } + + constexpr static char const * numStepsString() + { return "steps"; } + + constexpr static char const * outputString() + { return "output"; } + + constexpr static char const * baselineString() + { return "baseline"; } + }; + + integer m_numSteps; ///< Number of load steps + integer m_numColumns; ///< Number of columns in data table (depends on number of fluid phases) + integer m_numPhases; ///< Number of fluid phases + + string m_frictionName; ///< frictionType identifier + string m_outputFile; ///< Output file (optional, no output if not specified) + + array2d< real64 > m_table; ///< Table storing time-history of input/output + + Path m_baselineFile; ///< Baseline file (optional, for unit testing of solid models) + + enum columnKeys + { + TIME + }; ///< Enumeration of "input" column keys for readability + + static constexpr real64 m_baselineTol = 1e-3; ///< Comparison tolerance for baseline results +}; + + +} + +#endif //GEOS_FRICTIONDRIVER_HPP_ diff --git a/src/coreComponents/constitutiveDrivers/contact/FrictionDriverRunTest.hpp b/src/coreComponents/constitutiveDrivers/contact/FrictionDriverRunTest.hpp new file mode 100644 index 00000000000..eaaa5bae76a --- /dev/null +++ b/src/coreComponents/constitutiveDrivers/contact/FrictionDriverRunTest.hpp @@ -0,0 +1,106 @@ +/* + * ------------------------------------------------------------------------------------------------------------ + * 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 "constitutive/contact/FrictionDriver.hpp" +#include "physicsSolvers/solidMechanics/contact/FractureState.hpp" +#include "constitutive/solid/SolidFields.hpp" + + +namespace geos +{ + + +template< typename FRICTION_TYPE > +void +FrictionDriver::runTest( FRICTION_TYPE & relperm, + const arrayView2d< real64 > & table ) +{ + // get number of phases and components + + integer const numPhases = relperm.numFluidPhases(); + + // create kernel wrapper + + typename RELPERM_TYPE::KernelWrapper const kernelWrapper = relperm.createKernelWrapper(); + + // set saturation to user specified feed + // it is more convenient to provide input in molar, so perform molar to mass conversion here + + array2d< real64, compflow::LAYOUT_PHASE > saturationValues; + if( numPhases > 2 ) + { + saturationValues.resize(( m_numSteps + 1 ) * ( m_numSteps + 1 ), numPhases ); + } + else + { + saturationValues.resize( m_numSteps + 1, numPhases ); + } + integer const ipWater = relperm.getPhaseOrder()[PT::WATER]; + integer const ipOil = relperm.getPhaseOrder()[PT::OIL]; + integer const ipGas = relperm.getPhaseOrder()[PT::GAS]; + const localIndex offset = std::max( std::max( ipOil, ipWater ), std::max( ipOil, ipGas ) ) + 1; + + for( integer n = 0; n < table.size( 0 ); ++n ) + { + + + if( m_numPhases > 2 ) + { + saturationValues[n][ipWater] = table( n, ipWater + 1 ); + saturationValues[n][ipOil] = table( n, ipOil + 1 ); + saturationValues[n][ipGas] = table( n, ipGas + 1 ); + } + else//two-phase + { + if( ipWater < 0 ) + { + saturationValues[n][ipOil] = table( n, ipOil + 1 ); + saturationValues[n][ipGas] = table( n, ipGas + 1 ); + } + else if( ipGas < 0 ) + { + saturationValues[n][ipWater] = table( n, ipWater + 1 ); + saturationValues[n][ipOil] = table( n, ipOil + 1 ); + } + } + + } + + arrayView2d< real64 const, compflow::USD_PHASE > const saturation = saturationValues.toViewConst(); + + // perform relperm update using table (Swet,Snonwet) and save resulting total density, etc. + // note: column indexing should be kept consistent with output file header below. + + forAll< parallelDevicePolicy<> >( saturation.size( 0 ), + [numPhases, kernelWrapper, saturation, table, + offset] GEOS_HOST_DEVICE ( integer const n ) + { + kernelWrapper.update( 0, 0, saturation[n] ); + for( integer p = 0; p < numPhases; ++p ) + { + table( n, offset + 1 + p ) = kernelWrapper.relperm()( 0, 0, p ); + } + } ); + +} + + +} + + +#endif //GEOS_FRICTIONDRIVERRUNTEST_HPP_ diff --git a/src/coreComponents/constitutiveDrivers/contact/LogLevelsInfo.hpp b/src/coreComponents/constitutiveDrivers/contact/LogLevelsInfo.hpp new file mode 100644 index 00000000000..1a0319a7db4 --- /dev/null +++ b/src/coreComponents/constitutiveDrivers/contact/LogLevelsInfo.hpp @@ -0,0 +1,52 @@ +/* + * ------------------------------------------------------------------------------------------------------------ + * 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. + * ------------------------------------------------------------------------------------------------------------ + */ + +/** + * @file LogLevelsInfo.hpp + * This file contains common log level informations for PVT driver + */ + +#ifndef GEOS_CONSTITUTIVEDRIVERS_CONTACT_LOGLEVELSINFO_HPP_ +#define GEOS_CONSTITUTIVEDRIVERS_CONTACT_LOGLEVELSINFO_HPP_ + +#include "common/DataTypes.hpp" + +namespace geos +{ + +namespace logInfo +{ + +/** + * @name Common LogLevels info structures. They must comply with the `is_log_level_info` trait. + */ +///@{ + +/// @cond DO_NOT_DOCUMENT + +struct LogOutput +{ + static constexpr int getMinLogLevel() { return 1; } + static constexpr std::string_view getDescription() { return "Enable log output"; } +}; + +/// @endcond +///@} + +} + +} + +#endif // GEOS_CONSTITUTIVEDRIVERS_CONTACT_LOGLEVELSINFO_HPP_ From 67b1df99800130c68793fa6268473a488802cd80 Mon Sep 17 00:00:00 2001 From: jacques franc Date: Wed, 5 Nov 2025 17:45:17 +0100 Subject: [PATCH 2/4] first working v --- .../frictionDriver/frictionDriver_Coulomb.xml | 20 ++ .../frictionDriver/frictionDriver_base.xml | 62 +++++ inputFiles/frictionDriver/tables/jumps.geos | 6 + inputFiles/frictionDriver/tables/time.geos | 6 + .../frictionDriver/tables/tractions.geos | 6 + .../constitutiveDrivers/CMakeLists.txt | 4 + .../contact/FrictionDriver.cpp | 233 ++++++++++-------- .../contact/FrictionDriver.hpp | 43 ++-- .../contact/FrictionDriverRunTest.cpp | 26 ++ .../contact/FrictionDriverRunTest.hpp | 80 ++---- 10 files changed, 317 insertions(+), 169 deletions(-) create mode 100644 inputFiles/frictionDriver/frictionDriver_Coulomb.xml create mode 100644 inputFiles/frictionDriver/frictionDriver_base.xml create mode 100644 inputFiles/frictionDriver/tables/jumps.geos create mode 100644 inputFiles/frictionDriver/tables/time.geos create mode 100644 inputFiles/frictionDriver/tables/tractions.geos create mode 100644 src/coreComponents/constitutiveDrivers/contact/FrictionDriverRunTest.cpp diff --git a/inputFiles/frictionDriver/frictionDriver_Coulomb.xml b/inputFiles/frictionDriver/frictionDriver_Coulomb.xml new file mode 100644 index 00000000000..3a1c0a5b630 --- /dev/null +++ b/inputFiles/frictionDriver/frictionDriver_Coulomb.xml @@ -0,0 +1,20 @@ + + + + + + + + + + + + + diff --git a/inputFiles/frictionDriver/frictionDriver_base.xml b/inputFiles/frictionDriver/frictionDriver_base.xml new file mode 100644 index 00000000000..df5ef688301 --- /dev/null +++ b/inputFiles/frictionDriver/frictionDriver_base.xml @@ -0,0 +1,62 @@ + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + diff --git a/inputFiles/frictionDriver/tables/jumps.geos b/inputFiles/frictionDriver/tables/jumps.geos new file mode 100644 index 00000000000..3c922732f7a --- /dev/null +++ b/inputFiles/frictionDriver/tables/jumps.geos @@ -0,0 +1,6 @@ +0 +1e-3 +2e-3 +3e-3 +4e-3 +5e-3 diff --git a/inputFiles/frictionDriver/tables/time.geos b/inputFiles/frictionDriver/tables/time.geos new file mode 100644 index 00000000000..e8371f00609 --- /dev/null +++ b/inputFiles/frictionDriver/tables/time.geos @@ -0,0 +1,6 @@ +0 +1 +2 +3 +4 +5 diff --git a/inputFiles/frictionDriver/tables/tractions.geos b/inputFiles/frictionDriver/tables/tractions.geos new file mode 100644 index 00000000000..c80b1f6d1bd --- /dev/null +++ b/inputFiles/frictionDriver/tables/tractions.geos @@ -0,0 +1,6 @@ +0 +1e6 +2e6 +3e6 +4e6 +5e6 diff --git a/src/coreComponents/constitutiveDrivers/CMakeLists.txt b/src/coreComponents/constitutiveDrivers/CMakeLists.txt index de519ea2fdb..b64c154dd1d 100644 --- a/src/coreComponents/constitutiveDrivers/CMakeLists.txt +++ b/src/coreComponents/constitutiveDrivers/CMakeLists.txt @@ -29,6 +29,8 @@ set( constitutiveDrivers_headers relativePermeability/RelpermDriver.hpp relativePermeability/RelpermDriverRunTest.hpp solid/TriaxialDriver.hpp + contact/FrictionDriver.hpp + contact/FrictionDriverRunTest.hpp ) # # Specify all sources @@ -57,6 +59,8 @@ set( constitutiveDrivers_sources relativePermeability/RelpermDriverTableRelativeRunTest.cpp relativePermeability/RelpermDriverTableRelativeHysteresisRunTest.cpp solid/TriaxialDriver.cpp + contact/FrictionDriver.cpp + contact/FrictionDriverRunTest.cpp ) set( dependencyList ${parallelDeps} constitutive events ) diff --git a/src/coreComponents/constitutiveDrivers/contact/FrictionDriver.cpp b/src/coreComponents/constitutiveDrivers/contact/FrictionDriver.cpp index f51425517a9..0239141cddc 100644 --- a/src/coreComponents/constitutiveDrivers/contact/FrictionDriver.cpp +++ b/src/coreComponents/constitutiveDrivers/contact/FrictionDriver.cpp @@ -19,6 +19,10 @@ #include "constitutive/contact/FrictionBase.hpp" #include "constitutive/contact/FrictionSelector.hpp" +#include "functions/FunctionManager.hpp" +#include "functions/TableFunction.hpp" +#include + #include "FrictionDriver.hpp" namespace geos { @@ -33,11 +37,27 @@ TaskBase(name, parent) registerWrapper( viewKeyStruct::frictionNameString(), &m_frictionName ). setRTTypeName( rtTypes::CustomTypes::groupNameRef ). setInputFlag( InputFlags::REQUIRED ). - setDescription( "Relperm model to test" ); + setDescription( "Friction model to test" ); registerWrapper( viewKeyStruct::numStepsString(), &m_numSteps ). setInputFlag( InputFlags::REQUIRED ). - setDescription( "Number of saturation steps to take" ); + 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::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( "Number of increment step to take in both jumps and traction increments" ); + + registerWrapper( viewKeyStruct::phiString(), &m_phi ). + setInputFlag( InputFlags::INVALID). + setDescription( "Number of increment step to take in both jumps and traction increments" ); registerWrapper( viewKeyStruct::outputString(), &m_outputFile ). setInputFlag( InputFlags::OPTIONAL ). @@ -52,22 +72,67 @@ TaskBase(name, parent) addLogLevel< logInfo::LogOutput >(); } -void FrictionDriver::outputResults() +void FrictionDriver::compareWithBaseline() { - // TODO: improve file path output to grab command line -o directory - // for the moment, we just use the specified m_outputFile directly + // open baseline file - FILE * fp = fopen( m_outputFile.c_str(), "w" ); + std::ifstream file( m_baselineFile.c_str() ); + GEOS_THROW_IF( !file.is_open(), "Can't seem to open the baseline file " << m_baselineFile, InputError ); - fprintf( fp, "# column 1 = time\n" ); - fprintf( fp, "# columns %d-%d = phase vol fractions\n", 2, 1 + m_numPhases ); - fprintf( fp, "# columns %d-%d = phase relperm\n", 2 + m_numPhases, 1 + 2 * m_numPhases ); + // discard file header - if( ( m_numPhases == 2 && m_table.size( 1 ) > 5 ) || m_table.size( 1 ) > 7 ) + string line; + for( integer row=0; row < m_numColumns; ++row ) { - fprintf( fp, "# columns %d-%d = phase relperm (hyst)\n", 1 + 2 * m_numPhases, 1 + 3 * m_numPhases ); + getline( file, line ); } + // read data block. we assume the file size is consistent with m_table, + // but check for a premature end-of-file. we then compare results value by value. + // we ignore the newton iteration and residual columns, as those may be platform + // specific. + + real64 value; + real64 error; + + for( integer row=0; row < m_table.size( 0 ); ++row ) + { + for( integer col=0; col < m_table.size( 1 ); ++col ) + { + GEOS_THROW_IF( file.eof(), "Baseline file appears shorter than internal results", std::runtime_error ); + file >> value; + + error = fabs( m_table[row][col]-value ) / ( fabs( value )+1 ); + GEOS_THROW_IF( error > m_baselineTol, "Results do not match baseline at data row " << row+1 + << " (row " << row+10 << " with header)" + << " and column " << col+1, std::runtime_error ); + + } + } + + // check we actually reached the end of the baseline file + + file >> value; + GEOS_THROW_IF( !file.eof(), "Baseline file appears longer than internal results", std::runtime_error ); + + // success + + GEOS_LOG_LEVEL_RANK_0( logInfo::LogOutput, " Comparison ........ Internal results consistent with baseline." ); + + file.close(); +} + +void FrictionDriver::outputResults() +{ + // TODO: improve file path output to grab command line -o directory + // for the moment, we just use the specified m_outputFile directly + + FILE * fp = fopen( m_outputFile.c_str(), "w" ); + + fprintf( fp, "# column 1-3 = normal and in-plane displacement jump\n" ); + fprintf( fp, "# columns 4-6 = normal and in-place tractions\n"); + fprintf( fp, "# columns 7 = fracture state (0:Stick,1-2:[new]Slip,3:Open)\n"); + fprintf( fp, "# columns 8 = tau lim\n"); for( integer n = 0; n < m_table.size( 0 ); ++n ) { @@ -79,7 +144,6 @@ void FrictionDriver::outputResults() } fclose( fp ); - } @@ -104,7 +168,7 @@ bool FrictionDriver::execute( const geos::real64 GEOS_UNUSED_PARAM( time_n ), { // this code only makes sense in serial - GEOS_THROW_IF( MpiWrapper::commRank() > 0, "RelpermDriver should only be run in serial", std::runtime_error ); + GEOS_THROW_IF( MpiWrapper::commRank() > 0, "FrictionDriver should only be run in serial", std::runtime_error ); ConstitutiveManager @@ -112,8 +176,8 @@ bool FrictionDriver::execute( const geos::real64 GEOS_UNUSED_PARAM( time_n ), FrictionBase & baseFriction = constitutiveManager.getGroup< FrictionBase >( m_frictionName ); - GEOS_LOG_LEVEL_RANK_0( logInfo::LogOutput, "Launching Relperm Driver" ); - GEOS_LOG_LEVEL_RANK_0( logInfo::LogOutput, " Relperm .................. " << m_frictionName ); + 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, " No. of Phases .......... " << m_numPhases ); GEOS_LOG_LEVEL_RANK_0( logInfo::LogOutput, " Steps .................. " << m_numSteps ); @@ -128,13 +192,13 @@ bool FrictionDriver::execute( const geos::real64 GEOS_UNUSED_PARAM( time_n ), dataRepository::Group discretization( "discretization", &rootGroup ); discretization.resize( 1 ); // one element - baseRelperm.allocateConstitutiveData( discretization, 1 ); // one quadrature point + baseFriction.allocateConstitutiveData( discretization, 1 ); // one quadrature point constitutiveUpdatePassThru( baseFriction, [&]( auto & selectedFrictionModel ) { - using RELPERM_TYPE = TYPEOFREF( selectedFrictionModel ); - resizeTables< RELPERM_TYPE >(); - runTest< RELPERM_TYPE >( selectedFrictionModel, m_table ); + using FRICTION_TYPE = TYPEOFREF( selectedFrictionModel ); + resizeTables< FRICTION_TYPE >(); + runTest< FRICTION_TYPE >( selectedFrictionModel, m_table ); } ); // move table back to host for output @@ -154,103 +218,76 @@ bool FrictionDriver::execute( const geos::real64 GEOS_UNUSED_PARAM( time_n ), } -template< typename RELPERM_TYPE > +template< typename FRICTION_TYPE > void FrictionDriver::resizeTables() { ConstitutiveManager & constitutiveManager = this->getGroupByPath< ConstitutiveManager >( "/Problem/domain/Constitutive" ); FrictionBase - & baseRelperm = constitutiveManager.getGroup< FrictionBaseBase >( m_frictionName ); + & baseFriction = constitutiveManager.getGroup< FrictionBase >( m_frictionName ); - using PT = RelativePermeabilityBase::PhaseType; - integer const ipWater = baseRelperm.getPhaseOrder()[PT::WATER]; - integer const ipOil = baseRelperm.getPhaseOrder()[PT::OIL]; - integer const ipGas = baseRelperm.getPhaseOrder()[PT::GAS]; + // initialize table functions + FunctionManager & functionManager = FunctionManager::getInstance(); - real64 minSw = 0., minSnw = 0.; - if( baseRelperm.numFluidPhases() > 2 ) - { - minSw = baseRelperm.getWettingPhaseMinVolumeFraction(); - minSnw = baseRelperm.getNonWettingMinVolumeFraction(); - } - else - { - if( ipWater < 0 )// a.k.a o/g - { - minSw = 0; - minSnw = baseRelperm.getNonWettingMinVolumeFraction(); - } - else if( ipGas < 0 || ipOil < 0 )// a.k.a w/o or w/g - { - minSnw = 0; - minSw = baseRelperm.getWettingPhaseMinVolumeFraction(); - } - } + TableFunction & jumpFunction = functionManager.getGroup< TableFunction >( m_jumpFunctionName ); + TableFunction & tractionFunction = functionManager.getGroup< TableFunction >( m_tractionFunctionName ); - real64 const dSw = ( 1 - minSw - minSnw ) / m_numSteps; + jumpFunction.initializeFunction(); + tractionFunction.initializeFunction(); + + ArrayOfArraysView< real64 > Jcoordinates = jumpFunction.getCoordinates(); + real64 const minJump = Jcoordinates[0][0]; + real64 const maxJump = Jcoordinates[0][Jcoordinates.sizeOfArray( 0 )-1]; + real64 const dJ = (maxJump-minJump) / m_numSteps; + ArrayOfArraysView< real64 > Tcoordinates = jumpFunction.getCoordinates(); + real64 const minTrac = Tcoordinates[0][0]; + real64 const maxTrac = Tcoordinates[0][Tcoordinates.sizeOfArray( 0 )-1]; + real64 const dT = (maxTrac-minTrac) / m_numSteps; // set input columns + resizeTable< FRICTION_TYPE >(); - resizeTable< RELPERM_TYPE >(); - // 3-phase branch - if( m_numPhases > 2 ) - { - for( integer ni = 0; ni < m_numSteps + 1; ++ni ) - { - for( integer nj = 0; nj < m_numSteps + 1; ++nj ) - { - - integer index = ni * ( m_numSteps + 1 ) + nj; - m_table( index, TIME ) = minSw + index * dSw; - m_table( index, ipWater + 1 ) = minSw + nj * dSw; - m_table( index, ipGas + 1 ) = minSnw + ni * dSw; - m_table( index, ipOil + 1 ) = - 1. - m_table( index, ipWater + 1 ) - m_table( index, ipOil + 1 ); - } - } - } - else // 2-phase branch + // real64 const cohesion = baseFriction.getWrapper("cohesion"); + // real64 const frictionCoeff = baseFriction.getGroup("frictionCoefficient"); + + //All variation + for( integer nt = 0; nt < m_numSteps+1; ++nt ) { - for( integer ni = 0; ni < m_numSteps + 1; ++ni ) + + + for( integer nj = 0; nj < m_numSteps+1; ++nj ) { - integer index = ni; - m_table( index, TIME ) = minSw + index * dSw; - if( ipWater < 0 ) - { - m_table( index, ipGas + 1 ) = minSnw + ni * dSw; - m_table( index, ipOil + 1 ) = 1. - m_table( index, ipGas + 1 ); - } - else if( ipGas < 0 ) - { - m_table( index, ipWater + 1 ) = minSw + ni * dSw; - m_table( index, ipOil + 1 ) = 1. - m_table( index, ipWater + 1 ); - } - else if( ipOil < 0 ) - { - m_table( index, ipWater + 1 ) = minSw + ni * dSw; - m_table( index, ipGas + 1 ) = 1. - m_table( index, ipWater + 1 ); - } - } - } + integer index = nt * (m_numSteps+1) + nj; + m_table( index, NTRAC ) = (minTrac + nt*dT)*cos(m_theta * 180/M_PI); + m_table( index, STRAC0 ) = (minTrac + nt*dT)*sin(m_theta * 180/M_PI); + m_table( index, STRAC1 ) = (minTrac + nt*dT)*sin(m_theta * 180/M_PI); + + m_table( index, NJUMP ) = (minJump + nj*dJ)*cos(m_theta * 180/M_PI); + m_table( index, SLIP0 ) = (minJump + nj*dJ)*sin(m_theta * 180/M_PI); + m_table( index, SLIP1 ) = (minJump + nj*dJ)*sin(m_theta * 180/M_PI); + m_table( index, FS ) = fields::contact::FractureState::Stick; -} + m_table( index, TLIM ) = -1; + //Only for Coulomb + // m_table( index, TLIM ) = cohesion - m_table(index, NTRAC) * frictionCoeff; -// template< typename RELPERM_TYPE > -// std::enable_if_t< std::is_same< TableRelativePermeabilityHysteresis, RELPERM_TYPE >::value, void > -// RelpermDriver::resizeTable() -// { -// if( m_numPhases > 2 ) -// { -// m_table.resize( ( m_numSteps + 1 ) * ( m_numSteps + 1 ), 1 + 3 * m_numPhases ); -// } -// else -// { -// m_table.resize( m_numSteps + 1, 1 + 3 * m_numPhases ); -// } + } + } + +} -// } +template< typename FRICTION_TYPE > +void +FrictionDriver::resizeTable() +{ + m_table.resize((m_numSteps + 1)*(m_numSteps + 1), m_numColumns ); +} + +REGISTER_CATALOG_ENTRY( TaskBase, + FrictionDriver, + string const &, dataRepository::Group * const ) } \ No newline at end of file diff --git a/src/coreComponents/constitutiveDrivers/contact/FrictionDriver.hpp b/src/coreComponents/constitutiveDrivers/contact/FrictionDriver.hpp index ae02bb12a5c..c7e5f4ce6c2 100644 --- a/src/coreComponents/constitutiveDrivers/contact/FrictionDriver.hpp +++ b/src/coreComponents/constitutiveDrivers/contact/FrictionDriver.hpp @@ -29,7 +29,7 @@ class FrictionDriver : public TaskBase Group * const parent ); static string catalogName() - { return "frictionDriver"; } + { return "FrictionDriver"; } void postInputInitialization() override; @@ -51,10 +51,10 @@ class FrictionDriver : public TaskBase // runTest( friction_TYPE & friction, // const arrayView2d< real64, 1 > & table ); - // template< typename friction_TYPE > - // std::enable_if_t< !std::is_same< constitutive::TableRelativePermeabilityHysteresis, friction_TYPE >::value, void > - // runTest( friction_TYPE & friction, - // const arrayView2d< real64, 1 > & table ); + template< typename FRICTION_TYPE > + void + runTest( FRICTION_TYPE & friction, + const arrayView2d< real64, 1 > & table ); /** * @brief Ouput table to file for easy plotting @@ -69,11 +69,10 @@ class FrictionDriver : public TaskBase private: template< typename FRICTION_TYPE > - void resizeTables(); + void resizeTable(); - // template< typename friction_TYPE > - // std::enable_if_t< std::is_same< constitutive::TableRelativePermeabilityHysteresis, friction_TYPE >::value, void > - // resizeTable(); + template< typename FRICTION_TYPE > + void resizeTables(); // template< typename friction_TYPE > // std::enable_if_t< !std::is_same< constitutive::TableRelativePermeabilityHysteresis, friction_TYPE >::value, void > @@ -90,16 +89,34 @@ class FrictionDriver : public TaskBase constexpr static char const * numStepsString() { return "steps"; } + constexpr static char const * jumpFunctionString() + { return "jumpControl"; } + + 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 * outputString() { return "output"; } constexpr static char const * baselineString() { return "baseline"; } + }; integer m_numSteps; ///< Number of load steps - integer m_numColumns; ///< Number of columns in data table (depends on number of fluid phases) - integer m_numPhases; ///< Number of fluid phases + static integer const m_numColumns = 8; ///< Number of columns in dat + enum columnKeys { NJUMP, SLIP0, SLIP1, NTRAC, STRAC0, STRAC1, FS, TLIM }; + + string m_jumpFunctionName; ///< + string m_tractionFunctionName; ///< + + float m_theta, m_phi;///< x- and y-tilt of fault string m_frictionName; ///< frictionType identifier string m_outputFile; ///< Output file (optional, no output if not specified) @@ -108,10 +125,6 @@ class FrictionDriver : public TaskBase Path m_baselineFile; ///< Baseline file (optional, for unit testing of solid models) - enum columnKeys - { - TIME - }; ///< Enumeration of "input" column keys for readability static constexpr real64 m_baselineTol = 1e-3; ///< Comparison tolerance for baseline results }; diff --git a/src/coreComponents/constitutiveDrivers/contact/FrictionDriverRunTest.cpp b/src/coreComponents/constitutiveDrivers/contact/FrictionDriverRunTest.cpp new file mode 100644 index 00000000000..2d40f34235f --- /dev/null +++ b/src/coreComponents/constitutiveDrivers/contact/FrictionDriverRunTest.cpp @@ -0,0 +1,26 @@ +#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 &); + +template +void +FrictionDriver::runTest(constitutive::FrictionlessContact&, const arrayView2d &); + +template +void +FrictionDriver::runTest(constitutive::RateAndStateFriction>&, const arrayView2d &); + +template +void +FrictionDriver::runTest(constitutive::RateAndStateFriction>&, const arrayView2d &); + +} diff --git a/src/coreComponents/constitutiveDrivers/contact/FrictionDriverRunTest.hpp b/src/coreComponents/constitutiveDrivers/contact/FrictionDriverRunTest.hpp index eaaa5bae76a..cde385ff032 100644 --- a/src/coreComponents/constitutiveDrivers/contact/FrictionDriverRunTest.hpp +++ b/src/coreComponents/constitutiveDrivers/contact/FrictionDriverRunTest.hpp @@ -16,7 +16,7 @@ #ifndef GEOS_FRICTIONDRIVERRUNTEST_HPP_ #define GEOS_FRICTIONDRIVERRUNTEST_HPP_ -#include "constitutive/contact/FrictionDriver.hpp" +#include "constitutiveDrivers/contact/FrictionDriver.hpp" #include "physicsSolvers/solidMechanics/contact/FractureState.hpp" #include "constitutive/solid/SolidFields.hpp" @@ -27,79 +27,47 @@ namespace geos template< typename FRICTION_TYPE > void -FrictionDriver::runTest( FRICTION_TYPE & relperm, +FrictionDriver::runTest( FRICTION_TYPE & friction, const arrayView2d< real64 > & table ) { - // get number of phases and components - integer const numPhases = relperm.numFluidPhases(); + array2d< real64 > jumps, tractions; + jumps.resize(table.size(0),3); + tractions.resize(table.size(0),3); - // create kernel wrapper - - typename RELPERM_TYPE::KernelWrapper const kernelWrapper = relperm.createKernelWrapper(); - - // set saturation to user specified feed - // it is more convenient to provide input in molar, so perform molar to mass conversion here - - array2d< real64, compflow::LAYOUT_PHASE > saturationValues; - if( numPhases > 2 ) - { - saturationValues.resize(( m_numSteps + 1 ) * ( m_numSteps + 1 ), numPhases ); - } - else - { - saturationValues.resize( m_numSteps + 1, numPhases ); - } - integer const ipWater = relperm.getPhaseOrder()[PT::WATER]; - integer const ipOil = relperm.getPhaseOrder()[PT::OIL]; - integer const ipGas = relperm.getPhaseOrder()[PT::GAS]; - const localIndex offset = std::max( std::max( ipOil, ipWater ), std::max( ipOil, ipGas ) ) + 1; - - for( integer n = 0; n < table.size( 0 ); ++n ) + for( integer n = 0; n < table.size(0); ++n) { + jumps[n][0] = table(n,NJUMP); + jumps[n][1] = table(n,SLIP0); + jumps[n][2] = table(n,SLIP1); - - if( m_numPhases > 2 ) - { - saturationValues[n][ipWater] = table( n, ipWater + 1 ); - saturationValues[n][ipOil] = table( n, ipOil + 1 ); - saturationValues[n][ipGas] = table( n, ipGas + 1 ); - } - else//two-phase - { - if( ipWater < 0 ) - { - saturationValues[n][ipOil] = table( n, ipOil + 1 ); - saturationValues[n][ipGas] = table( n, ipGas + 1 ); - } - else if( ipGas < 0 ) - { - saturationValues[n][ipWater] = table( n, ipWater + 1 ); - saturationValues[n][ipOil] = table( n, ipOil + 1 ); - } - } + tractions[n][0] = table(n,NTRAC); + tractions[n][1] = table(n,STRAC0); + tractions[n][2] = table(n,STRAC1); } - arrayView2d< real64 const, compflow::USD_PHASE > const saturation = saturationValues.toViewConst(); - // perform relperm update using table (Swet,Snonwet) and save resulting total density, etc. - // note: column indexing should be kept consistent with output file header below. + // create kernel wrapper + typename FRICTION_TYPE::KernelWrapper const kernelWrapper = friction.createKernelUpdates(); - forAll< parallelDevicePolicy<> >( saturation.size( 0 ), - [numPhases, kernelWrapper, saturation, table, - offset] GEOS_HOST_DEVICE ( integer const n ) + forAll< parallelDevicePolicy<> >( 1, + [ kernelWrapper, table, jumps, tractions ] + GEOS_HOST_DEVICE ( integer const ei ) { - kernelWrapper.update( 0, 0, saturation[n] ); - for( integer p = 0; p < numPhases; ++p ) + for( integer i = 1; i < table.size(0) ; ++i ) { - table( n, offset + 1 + p ) = kernelWrapper.relperm()( 0, 0, p ); + integer fs = fields::contact::FractureState::Stick; + kernelWrapper.updateFractureState( jumps[i], + tractions[i], + fs ); + + table(i,FS) = fs; } } ); } - } From 833f25aea0e91a7e95aad8f5c00aedc10e4d7201 Mon Sep 17 00:00:00 2001 From: jacques franc Date: Wed, 5 Nov 2025 18:17:45 +0100 Subject: [PATCH 3/4] Add tau lim --- .../constitutive/contact/CoulombFriction.hpp | 8 ++++++++ .../contact/FrictionDriver.cpp | 16 +++++++++------- 2 files changed, 17 insertions(+), 7 deletions(-) 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/contact/FrictionDriver.cpp b/src/coreComponents/constitutiveDrivers/contact/FrictionDriver.cpp index 0239141cddc..1388b8ecc6c 100644 --- a/src/coreComponents/constitutiveDrivers/contact/FrictionDriver.cpp +++ b/src/coreComponents/constitutiveDrivers/contact/FrictionDriver.cpp @@ -22,6 +22,7 @@ #include "functions/FunctionManager.hpp" #include "functions/TableFunction.hpp" #include +#include #include "FrictionDriver.hpp" @@ -246,14 +247,17 @@ void FrictionDriver::resizeTables() // set input columns resizeTable< FRICTION_TYPE >(); - // real64 const cohesion = baseFriction.getWrapper("cohesion"); - // real64 const frictionCoeff = baseFriction.getGroup("frictionCoefficient"); + real64 cohesion = 0.; + real64 frictionCoeff = 0.; + if constexpr ( std::is_same_v< CoulombFriction, FRICTION_TYPE> ) { + + cohesion = dynamic_cast(&baseFriction)->getCohesion(); + frictionCoeff = dynamic_cast(&baseFriction)->getFrictionCoeff(); + } //All variation for( integer nt = 0; nt < m_numSteps+1; ++nt ) { - - for( integer nj = 0; nj < m_numSteps+1; ++nj ) { @@ -268,10 +272,8 @@ void FrictionDriver::resizeTables() m_table( index, FS ) = fields::contact::FractureState::Stick; - - m_table( index, TLIM ) = -1; //Only for Coulomb - // m_table( index, TLIM ) = cohesion - m_table(index, NTRAC) * frictionCoeff; + m_table( index, TLIM ) = cohesion - m_table(index, NTRAC) * frictionCoeff; } } From 49116212d639040b0c68953d4e90665f08d11b10 Mon Sep 17 00:00:00 2001 From: jacques franc Date: Thu, 6 Nov 2025 14:54:34 +0100 Subject: [PATCH 4/4] fixes --- .../contact/FrictionDriver.cpp | 42 +++++++++++-------- .../contact/FrictionDriver.hpp | 4 +- 2 files changed, 26 insertions(+), 20 deletions(-) diff --git a/src/coreComponents/constitutiveDrivers/contact/FrictionDriver.cpp b/src/coreComponents/constitutiveDrivers/contact/FrictionDriver.cpp index 1388b8ecc6c..1f475f27a45 100644 --- a/src/coreComponents/constitutiveDrivers/contact/FrictionDriver.cpp +++ b/src/coreComponents/constitutiveDrivers/contact/FrictionDriver.cpp @@ -130,10 +130,11 @@ void FrictionDriver::outputResults() FILE * fp = fopen( m_outputFile.c_str(), "w" ); - fprintf( fp, "# column 1-3 = normal and in-plane displacement jump\n" ); - fprintf( fp, "# columns 4-6 = normal and in-place tractions\n"); - fprintf( fp, "# columns 7 = fracture state (0:Stick,1-2:[new]Slip,3:Open)\n"); - fprintf( fp, "# columns 8 = tau lim\n"); + fprintf( fp, "# column 1 = index \n" ); + fprintf( fp, "# column 2-3 = normal and in-plane displacement jump\n" ); + fprintf( fp, "# columns 5-7 = normal and in-place tractions\n"); + fprintf( fp, "# columns 8 = fracture state (0:Stick,1-2:[new]Slip,3:Open)\n"); + fprintf( fp, "# columns 9 = tau lim\n"); for( integer n = 0; n < m_table.size( 0 ); ++n ) { @@ -236,16 +237,21 @@ void FrictionDriver::resizeTables() jumpFunction.initializeFunction(); tractionFunction.initializeFunction(); - ArrayOfArraysView< real64 > Jcoordinates = jumpFunction.getCoordinates(); - real64 const minJump = Jcoordinates[0][0]; - real64 const maxJump = Jcoordinates[0][Jcoordinates.sizeOfArray( 0 )-1]; - real64 const dJ = (maxJump-minJump) / m_numSteps; - ArrayOfArraysView< real64 > Tcoordinates = jumpFunction.getCoordinates(); - real64 const minTrac = Tcoordinates[0][0]; - real64 const maxTrac = Tcoordinates[0][Tcoordinates.sizeOfArray( 0 )-1]; - real64 const dT = (maxTrac-minTrac) / m_numSteps; + ArrayOfArraysView< real64 > coordinates = jumpFunction.getCoordinates(); + real64 const minTime = coordinates[0][0]; + real64 const maxTime = coordinates[0][coordinates.sizeOfArray( 0 )-1]; + real64 const dt = (maxTime-minTime) / m_numSteps; + // set input columns resizeTable< FRICTION_TYPE >(); + + // set time column + for( integer k=0; k