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 @@
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
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 );
+ }
+
+}
+ */