From d10121df8f56963686c742dcba74470f41fbf8d0 Mon Sep 17 00:00:00 2001 From: LilianeRA Date: Fri, 9 Oct 2026 18:25:14 -0300 Subject: [PATCH 1/2] Internal computation of the hydraulic aperture from the mechanical aperture using JRC; Adding unit test --- .../ParallelPlatesPermeability.cpp | 17 ++ .../ParallelPlatesPermeability.hpp | 33 +++- .../constitutive/unitTests/CMakeLists.txt | 3 +- .../testJointRoughnessCoefficient.cpp | 163 ++++++++++++++++++ 4 files changed, 210 insertions(+), 6 deletions(-) create mode 100644 src/coreComponents/constitutive/unitTests/testJointRoughnessCoefficient.cpp diff --git a/src/coreComponents/constitutive/permeability/ParallelPlatesPermeability.cpp b/src/coreComponents/constitutive/permeability/ParallelPlatesPermeability.cpp index 0d4b8923210..0947d023c10 100644 --- a/src/coreComponents/constitutive/permeability/ParallelPlatesPermeability.cpp +++ b/src/coreComponents/constitutive/permeability/ParallelPlatesPermeability.cpp @@ -40,6 +40,13 @@ ParallelPlatesPermeability::ParallelPlatesPermeability( string const & name, Gro setSizedFromParent( 0 ). setDescription( "Default value of the permeability normal to the surface. If not specified the permeability is updated using the cubic law. " ); + // JRC is valid for values between 0 and 20, but we will limit the range to 1-20 to avoid divsion by 0. + registerWrapper( viewKeyStruct::jointRoughnessCoefficientString(), &m_jointRoughnessCoefficient ). + setApplyDefaultValue( -1 ). + setInputFlag( InputFlags::OPTIONAL ). + setSizedFromParent( 0 ). + setDescription( "Default value of the joint roughness coefficient." ); + registerField< fields::permeability::dPerm_dDispJump >( &m_dPerm_dDispJump ); } @@ -60,6 +67,16 @@ void ParallelPlatesPermeability::postInputInitialization() { m_updateTransversalComponent = false; } + + // if the JRC value is not -1 (the flag for not using it), then check for validity. + GEOS_LOG_RANK_0(std::abs( m_jointRoughnessCoefficient + 1.0)); + if( std::abs( m_jointRoughnessCoefficient + 1.0) > 1e-10 ) + { + GEOS_THROW_IF( m_jointRoughnessCoefficient < 0 || m_jointRoughnessCoefficient > 20, + GEOS_FMT( "Joint Roughness Coefficient must be between 0 and 20. Provided value: {}", + m_jointRoughnessCoefficient ), + InputError, getDataContext() ); + } } void ParallelPlatesPermeability::initializeState() const diff --git a/src/coreComponents/constitutive/permeability/ParallelPlatesPermeability.hpp b/src/coreComponents/constitutive/permeability/ParallelPlatesPermeability.hpp index 80b5901714b..bf1ecc3d4ee 100644 --- a/src/coreComponents/constitutive/permeability/ParallelPlatesPermeability.hpp +++ b/src/coreComponents/constitutive/permeability/ParallelPlatesPermeability.hpp @@ -35,10 +35,13 @@ class ParallelPlatesPermeabilityUpdate : public PermeabilityBaseUpdate ParallelPlatesPermeabilityUpdate( arrayView3d< real64 > const & permeability, arrayView3d< real64 > const & dPerm_dPressure, arrayView4d< real64 > const & dPerm_dDispJump, - bool const updateTransversalComponent ) + bool const updateTransversalComponent, + real64 const jointRoughnessCoefficient ) : PermeabilityBaseUpdate( permeability, dPerm_dPressure ), m_dPerm_dDispJump( dPerm_dDispJump ), - m_numDimensionsToUpdate( 3 ) + m_numDimensionsToUpdate( 3 ), + m_jointRoughnessCoefficient( jointRoughnessCoefficient ), + m_printWarning(true) { m_numDimensionsToUpdate = updateTransversalComponent ? 3 : 2; } @@ -52,8 +55,23 @@ class ParallelPlatesPermeabilityUpdate : public PermeabilityBaseUpdate { GEOS_UNUSED_VAR( oldHydraulicAperture ); - real64 const perm = newHydraulicAperture*newHydraulicAperture / 12.0; - real64 const dPerm_dHydraulicAperture = newHydraulicAperture / 6.0; + real64 aperture = newHydraulicAperture; + if ( m_jointRoughnessCoefficient >= 1.0 ) + { + if( std::pow( newHydraulicAperture, 0.4 ) <= m_jointRoughnessCoefficient ) + { + aperture = (newHydraulicAperture*newHydraulicAperture * 1e6) / std::pow( m_jointRoughnessCoefficient, 2.5 ); + if ( aperture > newHydraulicAperture ) // invalid region of the plot + { + aperture = newHydraulicAperture; + GEOS_LOG_RANK_0_IF( m_printWarning, "Warning: computed aperture with JRC is greater than the hydraulic aperture. This is an invalid region of the plot. Using the original aperture value instead."); + m_printWarning = false; + } + } + } + + real64 const perm = aperture*aperture / 12.0; + real64 const dPerm_dHydraulicAperture = aperture / 6.0; for( int dim=0; dim < m_numDimensionsToUpdate; dim++ ) { @@ -99,6 +117,8 @@ class ParallelPlatesPermeabilityUpdate : public PermeabilityBaseUpdate arrayView4d< real64 > m_dPerm_dDispJump; int m_numDimensionsToUpdate; + real64 m_jointRoughnessCoefficient; + bool mutable m_printWarning; }; @@ -129,13 +149,15 @@ class ParallelPlatesPermeability : public PermeabilityBase return KernelWrapper( m_permeability, m_dPerm_dPressure, m_dPerm_dDispJump, - m_updateTransversalComponent ); + m_updateTransversalComponent, + m_jointRoughnessCoefficient ); } struct viewKeyStruct : public PermeabilityBase::viewKeyStruct { static constexpr char const * transversalPermeabilityString() { return "transversalPermeability"; } + static constexpr char const * jointRoughnessCoefficientString() { return "jointRoughnessCoefficient"; } }; protected: @@ -150,6 +172,7 @@ class ParallelPlatesPermeability : public PermeabilityBase array4d< real64 > m_dPerm_dDispJump; real64 m_transversalPermeability; + real64 m_jointRoughnessCoefficient; bool m_updateTransversalComponent; diff --git a/src/coreComponents/constitutive/unitTests/CMakeLists.txt b/src/coreComponents/constitutive/unitTests/CMakeLists.txt index 94276e254da..ef68eec13af 100644 --- a/src/coreComponents/constitutive/unitTests/CMakeLists.txt +++ b/src/coreComponents/constitutive/unitTests/CMakeLists.txt @@ -40,7 +40,8 @@ set( gtest_geosx_tests testStabilityTest2Comp.cpp testStabilityTest9Comp.cpp testRachfordRice.cpp - testTwoPhaseImmiscibleFluid.cpp ) + testTwoPhaseImmiscibleFluid.cpp + testJointRoughnessCoefficient.cpp ) set( dependencyList gtest blas lapack constitutive ${parallelDeps} ) diff --git a/src/coreComponents/constitutive/unitTests/testJointRoughnessCoefficient.cpp b/src/coreComponents/constitutive/unitTests/testJointRoughnessCoefficient.cpp new file mode 100644 index 00000000000..a618069e90d --- /dev/null +++ b/src/coreComponents/constitutive/unitTests/testJointRoughnessCoefficient.cpp @@ -0,0 +1,163 @@ +/* + * ------------------------------------------------------------------------------------------------------------ + * 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. + * ------------------------------------------------------------------------------------------------------------ + */ + +// Source includes +#include "constitutive/ConstitutiveManager.hpp" +#include "constitutive/permeability/ParallelPlatesPermeability.hpp" +#include "constitutive/permeability/PermeabilityFields.hpp" +#include "dataRepository/xmlWrapper.hpp" + +// TPL includes +#include +#include + +using namespace geos; +using namespace ::geos::constitutive; + + +TEST( JointRoughnessCoefficientTests, testJRC5 ) +{ + conduit::Node node; + dataRepository::Group rootGroup( "root", node ); + ConstitutiveManager constitutiveManager( "constitutive", &rootGroup ); + + real64 constexpr defaultAperture = 1.0e-4; // mechanical aperture defined in SurfaceElementRegion + + string const inputStream = + "" + " " + ""; + + xmlWrapper::xmlDocument xmlDocument; + xmlWrapper::xmlResult xmlResult = xmlDocument.loadString( inputStream ); + if( !xmlResult ) + { + GEOS_LOG_RANK_0( "XML parsed with errors!" ); + GEOS_LOG_RANK_0( "Error description: " << xmlResult.description()); + GEOS_LOG_RANK_0( "Error offset: " << xmlResult.offset ); + } + + xmlWrapper::xmlNode xmlConstitutiveNode = xmlDocument.getChild( "Constitutive" ); + constitutiveManager.processInputFileRecursive( xmlDocument, xmlConstitutiveNode ); + constitutiveManager.postInputInitializationRecursive(); + + ParallelPlatesPermeability & cm = constitutiveManager.getConstitutiveRelation< ParallelPlatesPermeability >( "fracturePerm" ); + + ParallelPlatesPermeability::KernelWrapper cmw = cm.createKernelWrapper(); + + { + // When JRC = 5 in this case, the new aperture should be in the invalid region. + // So the new aperture should be equal to the default aperture. + real64 newHydraulicAperture = cmw.computeApertureUsingJRC( defaultAperture ); + + EXPECT_DOUBLE_EQ( newHydraulicAperture, defaultAperture ); + } +} + + +TEST( JointRoughnessCoefficientTests, testJRC10 ) +{ + conduit::Node node; + dataRepository::Group rootGroup( "root", node ); + ConstitutiveManager constitutiveManager( "constitutive", &rootGroup ); + + real64 constexpr defaultAperture = 1.0e-4; // mechanical aperture defined in SurfaceElementRegion + + string const inputStream = + "" + " " + ""; + + xmlWrapper::xmlDocument xmlDocument; + xmlWrapper::xmlResult xmlResult = xmlDocument.loadString( inputStream ); + if( !xmlResult ) + { + GEOS_LOG_RANK_0( "XML parsed with errors!" ); + GEOS_LOG_RANK_0( "Error description: " << xmlResult.description()); + GEOS_LOG_RANK_0( "Error offset: " << xmlResult.offset ); + } + + xmlWrapper::xmlNode xmlConstitutiveNode = xmlDocument.getChild( "Constitutive" ); + constitutiveManager.processInputFileRecursive( xmlDocument, xmlConstitutiveNode ); + constitutiveManager.postInputInitializationRecursive(); + + ParallelPlatesPermeability & cm = constitutiveManager.getConstitutiveRelation< ParallelPlatesPermeability >( "fracturePerm" ); + + ParallelPlatesPermeability::KernelWrapper cmw = cm.createKernelWrapper(); + + { + // When JRC = 5 in this case, the new aperture should be in the invalid region. + // So the new aperture should be equal to the default aperture. + real64 newHydraulicAperture = cmw.computeApertureUsingJRC( defaultAperture ); + + EXPECT_DOUBLE_EQ( newHydraulicAperture, 3.162277660168379e-05 ); + } +} +/* +TEST( BartonBandisStressPathDrivenTests, testPressure ) +{ + conduit::Node node; + dataRepository::Group rootGroup( "root", node ); + ConstitutiveManager constitutiveManager( "constitutive", &rootGroup ); + + real64 constexpr referenceAperture = 1.0e-4; // in-situ + std::stringstream ss; + ss << std::scientific << std::setprecision(4) << referenceAperture; + + string const inputStream = + "" + " " + ""; + + xmlWrapper::xmlDocument xmlDocument; + xmlWrapper::xmlResult xmlResult = xmlDocument.loadString( inputStream ); + if( !xmlResult ) + { + GEOS_LOG_RANK_0( "XML parsed with errors!" ); + GEOS_LOG_RANK_0( "Error description: " << xmlResult.description()); + GEOS_LOG_RANK_0( "Error offset: " << xmlResult.offset ); + } + + xmlWrapper::xmlNode xmlConstitutiveNode = xmlDocument.getChild( "Constitutive" ); + constitutiveManager.processInputFileRecursive( xmlDocument, xmlConstitutiveNode ); + constitutiveManager.postInputInitializationRecursive(); + + BartonBandisStressPathDriven & cm = constitutiveManager.getConstitutiveRelation< BartonBandisStressPathDriven >( "hApertureModel" ); + + BartonBandisStressPathDriven::KernelWrapper cmw = cm.createKernelWrapper(); + + { + array1d < real64 > normal(3); + normal[0] = 0.0; + normal[1] = 0.0; + normal[2] = 1.0; + + real64 const newAperture = cmw.computeHydraulicAperture(13842265.4230671, normal); + EXPECT_DOUBLE_EQ( newAperture, 0.00022328650392488669 ); + } + +} + */ From 29b9a2d856a0779543d6043fdcdf96fcd8472d58 Mon Sep 17 00:00:00 2001 From: LilianeRA Date: Fri, 9 Oct 2026 18:38:23 -0300 Subject: [PATCH 2/2] Adding sample files. --- .../buckleyLeverett_JRC10.xml | 260 +++++++++++++++++ .../single_phase_edfm_10x1x1.xml | 265 ++++++++++++++++++ 2 files changed, 525 insertions(+) create mode 100644 inputFiles/jointRoughnessCoefficient/buckleyLeverett_JRC10.xml create mode 100644 inputFiles/jointRoughnessCoefficient/single_phase_edfm_10x1x1.xml diff --git a/inputFiles/jointRoughnessCoefficient/buckleyLeverett_JRC10.xml b/inputFiles/jointRoughnessCoefficient/buckleyLeverett_JRC10.xml new file mode 100644 index 00000000000..4071d5e5348 --- /dev/null +++ b/inputFiles/jointRoughnessCoefficient/buckleyLeverett_JRC10.xml @@ -0,0 +1,260 @@ + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + diff --git a/inputFiles/jointRoughnessCoefficient/single_phase_edfm_10x1x1.xml b/inputFiles/jointRoughnessCoefficient/single_phase_edfm_10x1x1.xml new file mode 100644 index 00000000000..391a61f75b4 --- /dev/null +++ b/inputFiles/jointRoughnessCoefficient/single_phase_edfm_10x1x1.xml @@ -0,0 +1,265 @@ + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + +