diff --git a/host-configs/ubuntu.cmake b/host-configs/ubuntu.cmake
new file mode 100644
index 00000000000..fabff86e1a7
--- /dev/null
+++ b/host-configs/ubuntu.cmake
@@ -0,0 +1,32 @@
+set( CONFIG_NAME "ubuntu" )
+
+# Set compilers path
+set(CMAKE_C_COMPILER "/usr/bin/gcc" CACHE PATH "") # This is typically something like /usr/bin/gcc ... or clang
+set(CMAKE_CXX_COMPILER "/usr/bin/g++" CACHE PATH "") # This is typically something like /usr/bin/g++ ... or clang++
+set(ENABLE_FORTRAN OFF CACHE BOOL "" FORCE)
+
+# Set paths to mpi
+set(ENABLE_MPI ON CACHE PATH "")
+set(MPI_C_COMPILER "/usr/bin/mpicc" CACHE PATH "") # This is typically something like /usr/bin/mpicc
+set(MPI_CXX_COMPILER "/usr/bin/mpicxx" CACHE PATH "") # This is typically something like /usr/bin/mpicxx
+set(MPIEXEC "/usr/bin/mpirun" CACHE PATH "") # This is typically something like /usr/bin/mpirun
+
+# Set paths to blas and lapack
+set( BLAS_LIBRARIES "/usr/lib/x86_64-linux-gnu/libblas.so" CACHE PATH "" FORCE ) # This is typically something like /usr/lib64/libblas.so
+set( LAPACK_LIBRARIES "/usr/lib/x86_64-linux-gnu/liblapack.so" CACHE PATH "" FORCE ) # This is typically something like /usr/lib64/liblapack.so
+
+# Cuda and openMP
+set( ENABLE_CUDA OFF CACHE PATH "" FORCE )
+set( ENABLE_OPENMP OFF CACHE PATH "" FORCE )
+
+# TPLs
+set( ENABLE_TRILINOS OFF CACHE PATH "" FORCE )
+set( ENABLE_CALIPER OFF CACHE PATH "" FORCE )
+set( ENABLE_DOXYGEN OFF CACHE BOOL "" FORCE)
+set( ENABLE_MATHPRESSO OFF CACHE BOOL "" FORCE )
+
+if(NOT ( EXISTS "${GEOS_TPL_DIR}" AND IS_DIRECTORY "${GEOS_TPL_DIR}" ) )
+ set(GEOS_TPL_DIR "${CMAKE_SOURCE_DIR}/../../../thirdPartyLibs/install-${CONFIG_NAME}-release" CACHE PATH "" FORCE )
+endif()
+
+include(${CMAKE_CURRENT_LIST_DIR}/tpls.cmake)
diff --git a/inputFiles/stressPathDrivenGeomechanics/fractureMatrixFlow_edfm_BBP_base.xml b/inputFiles/stressPathDrivenGeomechanics/fractureMatrixFlow_edfm_BBP_base.xml
new file mode 100644
index 00000000000..f50d208d13f
--- /dev/null
+++ b/inputFiles/stressPathDrivenGeomechanics/fractureMatrixFlow_edfm_BBP_base.xml
@@ -0,0 +1,158 @@
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
diff --git a/inputFiles/stressPathDrivenGeomechanics/fractureMatrixFlow_edfm_horizontalFrac_BBP_smoke.xml b/inputFiles/stressPathDrivenGeomechanics/fractureMatrixFlow_edfm_horizontalFrac_BBP_smoke.xml
new file mode 100644
index 00000000000..198e1abdcd4
--- /dev/null
+++ b/inputFiles/stressPathDrivenGeomechanics/fractureMatrixFlow_edfm_horizontalFrac_BBP_smoke.xml
@@ -0,0 +1,122 @@
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
diff --git a/src/coreComponents/constitutive/CMakeLists.txt b/src/coreComponents/constitutive/CMakeLists.txt
index 7606cfe25f6..23e2671fd1c 100644
--- a/src/coreComponents/constitutive/CMakeLists.txt
+++ b/src/coreComponents/constitutive/CMakeLists.txt
@@ -151,6 +151,7 @@ set( constitutive_headers
permeability/DamagePermeability.hpp
permeability/ExponentialDecayPermeability.hpp
permeability/ParallelPlatesPermeability.hpp
+ permeability/BartonBandisPermeability.hpp
permeability/PermeabilityBase.hpp
permeability/PermeabilityFields.hpp
permeability/PressurePermeability.hpp
@@ -303,6 +304,7 @@ set( constitutive_sources
permeability/DamagePermeability.cpp
permeability/ExponentialDecayPermeability.cpp
permeability/ParallelPlatesPermeability.cpp
+ permeability/BartonBandisPermeability.cpp
permeability/PermeabilityBase.cpp
permeability/PressurePermeability.cpp
permeability/ProppantPermeability.cpp
diff --git a/src/coreComponents/constitutive/ConstitutivePassThru.hpp b/src/coreComponents/constitutive/ConstitutivePassThru.hpp
index d78d7bb1bf8..ccad125f20d 100644
--- a/src/coreComponents/constitutive/ConstitutivePassThru.hpp
+++ b/src/coreComponents/constitutive/ConstitutivePassThru.hpp
@@ -49,6 +49,7 @@
#include "permeability/CarmanKozenyPermeability.hpp"
#include "permeability/ExponentialDecayPermeability.hpp"
#include "permeability/ParallelPlatesPermeability.hpp"
+#include "permeability/BartonBandisPermeability.hpp"
#include "permeability/PressurePermeability.hpp"
#include "permeability/ProppantPermeability.hpp"
#include "permeability/SlipDependentPermeability.hpp"
@@ -375,6 +376,7 @@ struct ConstitutivePassThru< CompressibleSolidBase >
CompressibleSolid< PressurePorosity, CarmanKozenyPermeability >,
CompressibleSolid< PressurePorosity, ExponentialDecayPermeability >,
CompressibleSolid< PressurePorosity, ParallelPlatesPermeability >,
+ CompressibleSolid< PressurePorosity, BartonBandisPermeability >,
CompressibleSolid< PressurePorosity, PressurePermeability >,
CompressibleSolid< PressurePorosity, SlipDependentPermeability >,
CompressibleSolid< PressurePorosity, WillisRichardsPermeability >
@@ -389,6 +391,7 @@ struct ConstitutivePassThru< CompressibleSolidBase >
CompressibleSolid< PressurePorosity, CarmanKozenyPermeability >,
CompressibleSolid< PressurePorosity, ExponentialDecayPermeability >,
CompressibleSolid< PressurePorosity, ParallelPlatesPermeability >,
+ CompressibleSolid< PressurePorosity, BartonBandisPermeability >,
CompressibleSolid< PressurePorosity, PressurePermeability >,
CompressibleSolid< PressurePorosity, SlipDependentPermeability >,
CompressibleSolid< PressurePorosity, WillisRichardsPermeability >
@@ -462,6 +465,7 @@ struct ConstitutivePassThru< CoupledSolidBase >
CompressibleSolid< PressurePorosity, CarmanKozenyPermeability >,
CompressibleSolid< PressurePorosity, ExponentialDecayPermeability >,
CompressibleSolid< PressurePorosity, ParallelPlatesPermeability >,
+ CompressibleSolid< PressurePorosity, BartonBandisPermeability >,
CompressibleSolid< PressurePorosity, PressurePermeability >,
CompressibleSolid< PressurePorosity, SlipDependentPermeability >,
CompressibleSolid< PressurePorosity, WillisRichardsPermeability >,
@@ -503,6 +507,7 @@ struct ConstitutivePassThru< CoupledSolidBase >
CompressibleSolid< PressurePorosity, CarmanKozenyPermeability >,
CompressibleSolid< PressurePorosity, ExponentialDecayPermeability >,
CompressibleSolid< PressurePorosity, ParallelPlatesPermeability >,
+ CompressibleSolid< PressurePorosity, BartonBandisPermeability >,
CompressibleSolid< PressurePorosity, PressurePermeability >,
CompressibleSolid< PressurePorosity, SlipDependentPermeability >,
CompressibleSolid< PressurePorosity, WillisRichardsPermeability >,
diff --git a/src/coreComponents/constitutive/permeability/BartonBandisPermeability.cpp b/src/coreComponents/constitutive/permeability/BartonBandisPermeability.cpp
new file mode 100644
index 00000000000..be39b892c00
--- /dev/null
+++ b/src/coreComponents/constitutive/permeability/BartonBandisPermeability.cpp
@@ -0,0 +1,102 @@
+/*
+ * ------------------------------------------------------------------------------------------------------------
+ * 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 BartonBandisPermeability.cpp
+ */
+
+#include "BartonBandisPermeability.hpp"
+
+#include "constitutive/permeability/PermeabilityFields.hpp"
+
+namespace geos
+{
+
+using namespace dataRepository;
+
+namespace constitutive
+{
+
+
+BartonBandisPermeability::BartonBandisPermeability( string const & name, Group * const parent ):
+ PermeabilityBase( name, parent ),
+ m_updateTransversalComponent( true )
+{
+ registerWrapper( viewKeyStruct::transversalPermeabilityString(), &m_transversalPermeability ).
+ setApplyDefaultValue( -1 ).
+ setInputFlag( InputFlags::OPTIONAL ).
+ setSizedFromParent( 0 ).
+ setDescription( "Default value of the permeability normal to the surface. If not specified the permeability is updated using the cubic law. " );
+
+ /// TODO: must become a required parameter.
+ registerWrapper( viewKeyStruct::apertureZeroString(), &m_aperture0 ).
+ setInputFlag( dataRepository::InputFlags::OPTIONAL ).
+ setApplyDefaultValue( 1e-6 ).
+ setDescription( "Reference hydraulic aperture. It is the aperture at zero normal stress." );
+
+ registerWrapper( viewKeyStruct::biotCoefficientString(), &m_biotCoefficient ).
+ setApplyDefaultValue( 1.0 ).
+ setInputFlag( InputFlags::REQUIRED ).
+ setDescription( "Biot coefficient." );
+ registerWrapper( viewKeyStruct::poissonRatioString(), &m_poissonRatio ).
+ setApplyDefaultValue( 0.3 ).
+ setInputFlag( InputFlags::REQUIRED ).
+ setDescription( "Poisson ratio." );
+ registerWrapper( viewKeyStruct::normalStiffnessString(), &m_normalStiffness ).
+ setInputFlag( InputFlags::REQUIRED ).
+ setDescription( "Normal stiffness: Kni." );
+ registerWrapper( viewKeyStruct::referencePressureString(), &m_referencePressure ).
+ setApplyDefaultValue( 1.0e5 ).
+ setInputFlag( InputFlags::REQUIRED ).
+ setDescription( "Reference pressure: p_0." );
+ registerWrapper( viewKeyStruct::referenceTotalStressString(), &m_referenceTotalStress ).
+ setApplyDefaultValue( { 85.0e6, 85.0e6, 105.0e6 } ).
+ setInputFlag( InputFlags::REQUIRED ).
+ setDescription( "Total stress at reference state: sigmaT_0." );
+}
+
+
+void BartonBandisPermeability::postInputInitialization()
+{
+ PermeabilityBase::postInputInitialization();
+
+ if( m_transversalPermeability > -1 )
+ {
+ m_updateTransversalComponent = false;
+ }
+}
+
+void BartonBandisPermeability::initializeState() const
+{
+ localIndex const numE = m_permeability.size( 0 );
+ integer constexpr numQuad = 1; // NOTE: enforcing 1 quadrature point
+
+ auto permView = m_permeability.toView();
+
+ real64 const transversalPerm = m_transversalPermeability;
+
+ forAll< parallelDevicePolicy<> >( numE, [=] GEOS_HOST_DEVICE ( localIndex const ei )
+ {
+ for( localIndex q = 0; q < numQuad; ++q )
+ {
+ permView[ei][q][2] = transversalPerm;
+ }
+ } );
+}
+
+REGISTER_CATALOG_ENTRY( ConstitutiveBase, BartonBandisPermeability, string const &, Group * const )
+
+}
+} /* namespace geos */
diff --git a/src/coreComponents/constitutive/permeability/BartonBandisPermeability.hpp b/src/coreComponents/constitutive/permeability/BartonBandisPermeability.hpp
new file mode 100644
index 00000000000..b66eb9d473f
--- /dev/null
+++ b/src/coreComponents/constitutive/permeability/BartonBandisPermeability.hpp
@@ -0,0 +1,257 @@
+/*
+ * ------------------------------------------------------------------------------------------------------------
+ * 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 BartonBandisPermeability.hpp
+ */
+
+#ifndef GEOS_CONSTITUTIVE_PERMEABILITY_BARTONBANDISPERMEABILITY_HPP_
+#define GEOS_CONSTITUTIVE_PERMEABILITY_BARTONBANDISPERMEABILITY_HPP_
+
+#include "constitutive/permeability/PermeabilityBase.hpp"
+#include "LvArray/src/genericTensorOps.hpp"
+
+
+namespace geos
+{
+namespace constitutive
+{
+
+class BartonBandisPermeabilityUpdate : public PermeabilityBaseUpdate
+{
+public:
+
+ BartonBandisPermeabilityUpdate( arrayView3d< real64 > const & permeability,
+ arrayView3d< real64 > const & dPerm_dPressure,
+ bool const updateTransversalComponent,
+ real64 const aperture0,
+ real64 const biot,
+ real64 const poisson,
+ real64 const normalStiffness,
+ real64 const referencePressure,
+ R1Tensor const &referenceTotalStress )
+ : PermeabilityBaseUpdate( permeability, dPerm_dPressure ),
+ m_numDimensionsToUpdate( 3 ),
+ m_aperture0( aperture0 ),
+ m_biotCoefficient( biot ),
+ m_poissonRatio( poisson ),
+ m_normalStiffness( normalStiffness ), // Kni
+ m_referencePressure( referencePressure ),
+ m_referenceTotalStress( referenceTotalStress )
+ {
+ m_numDimensionsToUpdate = updateTransversalComponent ? 3 : 2;
+
+ m_referenceEffectiveStress[0] = m_referenceTotalStress[0] - m_biotCoefficient*m_referencePressure;
+ m_referenceEffectiveStress[1] = m_referenceTotalStress[1] - m_biotCoefficient*m_referencePressure;
+ m_referenceEffectiveStress[2] = m_referenceTotalStress[2] - m_biotCoefficient*m_referencePressure;
+ }
+
+ GEOS_HOST_DEVICE
+ void compute( real64 const & oldHydraulicAperture,
+ real64 const & newHydraulicAperture,
+ arraySlice1d< real64 > const & permeability,
+ real64 & dPerm_dHydraulicAperture ) const
+ {
+ GEOS_UNUSED_VAR( oldHydraulicAperture );
+
+ real64 const perm = newHydraulicAperture*newHydraulicAperture / 12.0;
+ dPerm_dHydraulicAperture = newHydraulicAperture / 6.0;
+
+ for( int dim=0; dim < m_numDimensionsToUpdate; dim++ )
+ {
+ permeability[dim] = perm;
+ }
+ }
+
+ GEOS_HOST_DEVICE
+ virtual void updateFromPressureApertureAndNormal( localIndex const k,
+ localIndex const q,
+ real64 const & pressure,
+ real64 const & oldHydraulicAperture,
+ real64 const & newHydraulicAperture,
+ arraySlice1d< real64 const > const & normal,
+ real64 const & dHydraulicAperture_dNormalJump ) const override final
+ {
+ GEOS_UNUSED_VAR( q, dHydraulicAperture_dNormalJump);
+ // compute effective normal stress on the fracture
+ real64 dStress_dPressure = -1.0;
+ real64 const fractureStress = computeFractureStress(pressure, normal, dStress_dPressure);
+ // compute new aperture using Barton Bandis model
+ real64 dAperture_dStress = -1.0;
+ real64 hydraulicAperture = computeHydraulicAperture(pressure, fractureStress, normal, dAperture_dStress, k);
+
+ real64 dPerm_dHydraulicAperture = -1.0;
+ compute( oldHydraulicAperture,
+ hydraulicAperture,
+ m_permeability[k][0],
+ dPerm_dHydraulicAperture );
+
+ real64 const dPerm_dPressure = dPerm_dHydraulicAperture * dAperture_dStress * dStress_dPressure;
+ for( localIndex i=0; i < m_permeability[k][0].size(); i++ ) // size = 3
+ {
+ m_dPerm_dPressure[k][0][i] = dPerm_dPressure;
+ }
+ }
+
+
+
+ GEOS_HOST_DEVICE
+ GEOS_FORCE_INLINE
+ real64 computeHydraulicAperture( real64 const pressure, real64 const normalComponentOfStressOnFracture,
+ arraySlice1d< real64 const > const & normal, real64 & dAperture_dStress, int k ) const
+ {
+ real64 const referenceTotalStress[ 3 ] = LVARRAY_TENSOROPS_INIT_LOCAL_3 (m_referenceTotalStress);
+ real64 const biot_pressure = m_biotCoefficient * m_referencePressure; // biot is alpha in the equations
+
+ // Computation of maximum fracture closure (Barton-Bandis parameter)
+ // Fracture traction via Terzaghi's Principle
+ real64 sigma_c0[3] = {0.0};
+ LvArray::tensorOps::hadamardProduct< 3 >( sigma_c0, referenceTotalStress, normal );
+ LvArray::tensorOps::scaledAdd< 3 >(sigma_c0, normal, -biot_pressure);
+
+ real64 const sigma_n0 = LvArray::tensorOps::AiBi< 3 >( sigma_c0, normal );
+
+ // \dfrac{-K_{ni}a_0 + \sqrt{(K_{ni}a_0)^2 + 4K_{ni}a_0\sigma_{n0}}}{2K_{ni}}.
+ real64 const normalStiffApertureProduct = m_normalStiffness * m_aperture0;
+ real64 const sqrtTerm = normalStiffApertureProduct * (normalStiffApertureProduct + 4.0*sigma_n0) ;
+ real64 const g0 = ( -normalStiffApertureProduct + std::sqrt(sqrtTerm) ) / (2.0 * m_normalStiffness);
+
+ real64 const maximumFractureClosure = g0 + m_aperture0; // V_m or a_m -> aperture at stress-free state
+
+ // Normal effective stress on the fracture
+ // g_n(\sigma_n) = \dfrac{\sigma_n * V_m}{K_{ni} * V_m + \sigma_n}
+ real64 const fractureClosure = (normalComponentOfStressOnFracture*maximumFractureClosure) / (m_normalStiffness*maximumFractureClosure + normalComponentOfStressOnFracture);
+
+ real64 const newHydraulicAperture = maximumFractureClosure - fractureClosure;
+
+ // derivative
+ // \frac{da}{d\sigma}(\sigma) = -\dfrac{K_{ni} V_m^2}{(K_{ni} V_m + \sigma)^2}
+ real64 const normalStiffMaxClosureProduct = m_normalStiffness*maximumFractureClosure;
+ real64 const denom = normalStiffMaxClosureProduct + normalComponentOfStressOnFracture;
+ dAperture_dStress = -(normalStiffMaxClosureProduct*maximumFractureClosure) / (denom * denom);
+
+ return newHydraulicAperture;
+ }
+
+private:
+
+ int m_numDimensionsToUpdate;
+
+
+ real64 m_aperture0;
+
+ real64 m_biotCoefficient;
+ real64 m_poissonRatio;
+ real64 m_normalStiffness; // Kni
+ real64 m_referencePressure; // p_0
+
+ R1Tensor m_referenceTotalStress; // sigmaT_0 computed analytically
+ R1Tensor m_referenceEffectiveStress; // sigma_0
+
+ GEOS_HOST_DEVICE
+ GEOS_FORCE_INLINE
+ real64 computeFractureStress( real64 const pressure, arraySlice1d< real64 const > const & normal, real64 & dStress_dPressure ) const
+ {
+ real64 const deltaSigmaZ = m_biotCoefficient * (pressure - m_referencePressure);
+ real64 const poisson_deltaSigma = deltaSigmaZ * m_poissonRatio/(1.0 - m_poissonRatio);
+ // sigma: matrix diagonal
+ real64 effectiveStress[3] = { m_referenceEffectiveStress[0] - poisson_deltaSigma,
+ m_referenceEffectiveStress[1] - poisson_deltaSigma,
+ m_referenceEffectiveStress[2] - deltaSigmaZ };
+ real64 effectiveStressOnFracture[3] = {0.0}; // sigma_c
+ LvArray::tensorOps::hadamardProduct< 3 >( effectiveStressOnFracture, normal, effectiveStress );
+ real64 const normalComponentOfStressOnFracture = LvArray::tensorOps::AiBi< 3 >(effectiveStressOnFracture, normal); // sigmaN_N
+
+ // derivative
+ dStress_dPressure = -m_biotCoefficient;
+
+ return normalComponentOfStressOnFracture;
+ }
+
+};
+
+
+class BartonBandisPermeability : public PermeabilityBase
+{
+public:
+
+ BartonBandisPermeability( string const & name, dataRepository::Group * const parent );
+
+ static string catalogName() { return "BartonBandisPermeability"; }
+
+ virtual string getCatalogName() const override { return catalogName(); }
+
+ virtual void initializeState() const override final;
+
+ /// Type of kernel wrapper for in-kernel update
+ using KernelWrapper = BartonBandisPermeabilityUpdate;
+
+ /**
+ * @brief Create an update kernel wrapper.
+ * @return the wrapper
+ */
+ KernelWrapper createKernelWrapper() const
+ {
+ return KernelWrapper( m_permeability,
+ m_dPerm_dPressure,
+ m_updateTransversalComponent,
+ m_aperture0,
+ m_biotCoefficient,
+ m_poissonRatio,
+ m_normalStiffness,
+ m_referencePressure,
+ m_referenceTotalStress);
+ }
+
+
+ struct viewKeyStruct : public PermeabilityBase::viewKeyStruct
+ {
+ static constexpr char const * transversalPermeabilityString() { return "transversalPermeability"; }
+
+ /// string/key for aperture under zero normal stress
+ static constexpr char const * apertureZeroString() { return "referenceAperture"; }
+ static constexpr char const * biotCoefficientString() { return "biotCoefficient"; }
+ static constexpr char const * poissonRatioString() { return "poissonRatio"; }
+ static constexpr char const * normalStiffnessString() { return "normalStiffness"; }
+ static constexpr char const * referencePressureString() { return "referencePressure"; }
+ static constexpr char const * referenceTotalStressString() { return "referenceTotalStress"; }
+ };
+
+protected:
+
+ virtual void postInputInitialization() override;
+
+private:
+
+ real64 m_transversalPermeability;
+
+ bool m_updateTransversalComponent;
+
+ /// Reference hydraulic aperture. Aperture at zero normal stress
+ real64 m_aperture0; /// TODO: this will replace what is currently called defaultAperture.
+ real64 m_biotCoefficient;
+ real64 m_poissonRatio;
+ real64 m_normalStiffness; // Kni
+ real64 m_referencePressure; // p_0
+
+ R1Tensor m_referenceTotalStress; // sigmaT_0 computed analytically
+};
+
+}/* namespace constitutive */
+
+} /* namespace geos */
+
+
+#endif //GEOS_CONSTITUTIVE_PERMEABILITY_BARTONBANDISPERMEABILITY_HPP_
diff --git a/src/coreComponents/constitutive/permeability/PermeabilityBase.hpp b/src/coreComponents/constitutive/permeability/PermeabilityBase.hpp
index 190b8f1b4f6..dfcde6ddba7 100644
--- a/src/coreComponents/constitutive/permeability/PermeabilityBase.hpp
+++ b/src/coreComponents/constitutive/permeability/PermeabilityBase.hpp
@@ -66,6 +66,18 @@ class PermeabilityBaseUpdate
GEOS_UNUSED_VAR( k, q, oldHydraulicAperture, newHydraulicAperture, dHydraulicAperture_dNormalJump );
}
+ GEOS_HOST_DEVICE
+ virtual void updateFromPressureApertureAndNormal( localIndex const k,
+ localIndex const q,
+ real64 const & pressure,
+ real64 const & oldHydraulicAperture,
+ real64 const & newHydraulicAperture,
+ arraySlice1d< real64 const > const & normal,
+ real64 const & dHydraulicAperture_dNormalJump ) const
+ {
+ GEOS_UNUSED_VAR( k, q, pressure, oldHydraulicAperture, newHydraulicAperture, normal, dHydraulicAperture_dNormalJump );
+ }
+
GEOS_HOST_DEVICE
virtual void updateFromApertureAndShearDisplacement( localIndex const k,
localIndex const q,
diff --git a/src/coreComponents/constitutive/solid/CompressibleSolid.cpp b/src/coreComponents/constitutive/solid/CompressibleSolid.cpp
index 5519d77543a..04791b16c44 100644
--- a/src/coreComponents/constitutive/solid/CompressibleSolid.cpp
+++ b/src/coreComponents/constitutive/solid/CompressibleSolid.cpp
@@ -24,6 +24,7 @@
#include "constitutive/permeability/CarmanKozenyPermeability.hpp"
#include "constitutive/permeability/ExponentialDecayPermeability.hpp"
#include "constitutive/permeability/ParallelPlatesPermeability.hpp"
+#include "constitutive/permeability/BartonBandisPermeability.hpp"
#include "constitutive/permeability/PressurePermeability.hpp"
#include "constitutive/permeability/SlipDependentPermeability.hpp"
#include "constitutive/permeability/WillisRichardsPermeability.hpp"
@@ -48,6 +49,7 @@ typedef CompressibleSolid< PressurePorosity, CarmanKozenyPermeability > Compress
typedef CompressibleSolid< PressurePorosity, PressurePermeability > CompressibleRockPressurePerm;
typedef CompressibleSolid< PressurePorosity, ExponentialDecayPermeability > FaultED;
typedef CompressibleSolid< PressurePorosity, ParallelPlatesPermeability > FractureRock;
+typedef CompressibleSolid< PressurePorosity, BartonBandisPermeability > FractureRockBB;
typedef CompressibleSolid< PressurePorosity, SlipDependentPermeability > Fault;
typedef CompressibleSolid< PressurePorosity, WillisRichardsPermeability > FaultWR;
@@ -55,6 +57,7 @@ REGISTER_CATALOG_ENTRY( ConstitutiveBase, CompressibleRockConstant, string const
REGISTER_CATALOG_ENTRY( ConstitutiveBase, CompressibleRockCK, string const &, Group * const )
REGISTER_CATALOG_ENTRY( ConstitutiveBase, CompressibleRockPressurePerm, string const &, Group * const )
REGISTER_CATALOG_ENTRY( ConstitutiveBase, FractureRock, string const &, Group * const )
+REGISTER_CATALOG_ENTRY( ConstitutiveBase, FractureRockBB, string const &, Group * const )
REGISTER_CATALOG_ENTRY( ConstitutiveBase, FaultED, string const &, Group * const )
REGISTER_CATALOG_ENTRY( ConstitutiveBase, Fault, string const &, Group * const )
REGISTER_CATALOG_ENTRY( ConstitutiveBase, FaultWR, string const &, Group * const )
diff --git a/src/coreComponents/constitutive/solid/CompressibleSolid.hpp b/src/coreComponents/constitutive/solid/CompressibleSolid.hpp
index fb1832b2804..5028feaaa9c 100644
--- a/src/coreComponents/constitutive/solid/CompressibleSolid.hpp
+++ b/src/coreComponents/constitutive/solid/CompressibleSolid.hpp
@@ -75,6 +75,20 @@ class CompressibleSolidUpdates : public CoupledSolidUpdates< NullModel, PORO_TYP
m_permUpdate.updateFromAperture( k, q, oldHydraulicAperture, newHydraulicAperture, dHydraulicAperture_dNormalJump );
}
+ GEOS_HOST_DEVICE
+ void updateStateFromPressureApertureAndNormal( localIndex const k,
+ localIndex const q,
+ real64 const & pressure,
+ real64 const & oldHydraulicAperture,
+ real64 const & newHydraulicAperture,
+ array1d< real64 > const & normal ) const
+ {
+ real64 const temperature = 0;
+ real64 const dHydraulicAperture_dNormalJump = 1.0;
+ m_porosityUpdate.updateFromPressureAndTemperature( k, q, pressure, temperature );
+ m_permUpdate.updateFromPressureApertureAndNormal( k, q, pressure, oldHydraulicAperture, newHydraulicAperture, normal, dHydraulicAperture_dNormalJump);
+ }
+
GEOS_HOST_DEVICE
void updateStateFromPressureApertureJumpAndTraction( localIndex const k,
localIndex const q,
diff --git a/src/coreComponents/physicsSolvers/fluidFlow/FlowSolverBase.cpp b/src/coreComponents/physicsSolvers/fluidFlow/FlowSolverBase.cpp
index 3a948b6fee9..235d780cb57 100644
--- a/src/coreComponents/physicsSolvers/fluidFlow/FlowSolverBase.cpp
+++ b/src/coreComponents/physicsSolvers/fluidFlow/FlowSolverBase.cpp
@@ -113,6 +113,31 @@ void updatePorosityAndPermeabilityFromPressureAndAperture( POROUSWRAPPER_TYPE po
} );
}
+template< typename POROUSWRAPPER_TYPE >
+void updatePorosityAndPermeabilityFromPressureApertureAndNormal( POROUSWRAPPER_TYPE porousWrapper,
+ SurfaceElementSubRegion & subRegion,
+ arrayView1d< real64 const > const & pressure,
+ arrayView1d< real64 const > const & oldHydraulicAperture,
+ arrayView1d< real64 const > const & newHydraulicAperture,
+ arrayView2d< real64 const > const & normalVector)
+{
+ forAll< parallelDevicePolicy<> >( subRegion.size(), [=] GEOS_DEVICE ( localIndex const k )
+ {
+ array1d < real64 > normal(3);
+ normal[0] = normalVector[k][0];
+ normal[1] = normalVector[k][1];
+ normal[2] = normalVector[k][2];
+ for( localIndex q = 0; q < porousWrapper.numGauss(); ++q )
+ {
+ porousWrapper.updateStateFromPressureApertureAndNormal( k, q,
+ pressure[k],
+ oldHydraulicAperture[k],
+ newHydraulicAperture[k],
+ normal );
+ }
+ } );
+}
+
template< typename POROUSWRAPPER_TYPE >
void updatePorosityAndPermeabilityFromPressureApertureJumpAndTraction( POROUSWRAPPER_TYPE porousWrapper,
SurfaceElementSubRegion & subRegion,
@@ -677,6 +702,7 @@ void FlowSolverBase::updatePorosityAndPermeability( SurfaceElementSubRegion & su
GEOS_MARK_FUNCTION;
arrayView1d< real64 const > const & pressure = subRegion.getField< flow::pressure >();
+ arrayView2d< real64 const > const & normalVector = subRegion.getField< fields::normalVector >(); // mesh/MeshFields.hp
arrayView1d< real64 const > const newHydraulicAperture = subRegion.getField< flow::hydraulicAperture >();
arrayView1d< real64 const > const oldHydraulicAperture = subRegion.getField< flow::aperture0 >();
@@ -707,12 +733,14 @@ void FlowSolverBase::updatePorosityAndPermeability( SurfaceElementSubRegion & su
}
else
{
- updatePorosityAndPermeabilityFromPressureAndAperture( porousWrapper, subRegion, pressure, oldHydraulicAperture, newHydraulicAperture );
+ // updatePorosityAndPermeabilityFromPressureAndAperture( porousWrapper, subRegion, pressure, oldHydraulicAperture, newHydraulicAperture );
+ updatePorosityAndPermeabilityFromPressureApertureAndNormal( porousWrapper, subRegion, pressure, oldHydraulicAperture, newHydraulicAperture, normalVector );
}
}
else
{
- updatePorosityAndPermeabilityFromPressureAndAperture( porousWrapper, subRegion, pressure, oldHydraulicAperture, newHydraulicAperture );
+ // updatePorosityAndPermeabilityFromPressureAndAperture( porousWrapper, subRegion, pressure, oldHydraulicAperture, newHydraulicAperture );
+ updatePorosityAndPermeabilityFromPressureApertureAndNormal( porousWrapper, subRegion, pressure, oldHydraulicAperture, newHydraulicAperture, normalVector );
}
} );
diff --git a/src/coreComponents/schema/schema.xsd b/src/coreComponents/schema/schema.xsd
index e78160e9ce0..bc1c457712d 100644
--- a/src/coreComponents/schema/schema.xsd
+++ b/src/coreComponents/schema/schema.xsd
@@ -222,18 +222,10 @@
-
-
-
-
-
-
-
-
@@ -667,6 +659,10 @@
+
+
+
+
@@ -743,6 +739,10 @@
+
+
+
+
@@ -1620,26 +1620,10 @@ stress - traction is applied to the faces as specified by the inner product of i
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
@@ -1648,18 +1632,6 @@ stress - traction is applied to the faces as specified by the inner product of i
-
-
-
-
-
-
-
-
-
-
-
-
@@ -6765,6 +6737,7 @@ Information output from lower logLevels is added with the desired log level
+
@@ -6784,6 +6757,7 @@ Information output from lower logLevels is added with the desired log level
+
@@ -6881,6 +6855,24 @@ Information output from lower logLevels is added with the desired log level
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
@@ -7459,6 +7451,18 @@ if the table is requested to be output in the log, and it is too large, a CSV fi
+
+
+
+
+
+
+
+
+
+
+
+
diff --git a/src/coreComponents/schema/schema.xsd.other b/src/coreComponents/schema/schema.xsd.other
index 12b9eae3e1d..41a98cbcad8 100644
--- a/src/coreComponents/schema/schema.xsd.other
+++ b/src/coreComponents/schema/schema.xsd.other
@@ -351,15 +351,11 @@ A field can represent a physical variable. (pressure, temperature, global compos
-
-
-
-
@@ -526,7 +522,7 @@ A field can represent a physical variable. (pressure, temperature, global compos
-
+
@@ -1609,13 +1605,14 @@ A field can represent a physical variable. (pressure, temperature, global compos
-
+
+
@@ -1635,6 +1632,7 @@ A field can represent a physical variable. (pressure, temperature, global compos
+
@@ -1725,6 +1723,12 @@ A field can represent a physical variable. (pressure, temperature, global compos
+
+
+
+
+
+
@@ -2421,6 +2425,7 @@ A field can represent a physical variable. (pressure, temperature, global compos
+
diff --git a/src/docs/sphinx/basicExamples/stressPathDriven/Example.rst b/src/docs/sphinx/basicExamples/stressPathDriven/Example.rst
new file mode 100644
index 00000000000..e1bfb2b7731
--- /dev/null
+++ b/src/docs/sphinx/basicExamples/stressPathDriven/Example.rst
@@ -0,0 +1,172 @@
+.. _stressPathDrivenExample:
+
+
+#############################################################
+Stress-Path Driven Simulation
+#############################################################
+
+
+**Context**
+
+In this example, we demonstrate how to set up a single-phase flow simulation driven by the oedometric stress path. The domain contains a single fracture, whose aperture is updated according to the Barton–Bandis model.
+
+**Objectives**
+
+At the end of this example, you will know:
+
+ - how to set up the Barton-Bandis parameters for the oedometric stress path simulation
+ - how to update the fracture aperture without a mechanics solver
+
+
+**Input file**
+
+The XML file for this test case is located at:
+
+.. code-block:: console
+
+ inputFiles/stressPathDrivenGeomechanics/fractureMatrixFlow_edfm_horizontalFrac_BBP_smoke.xml
+
+
+------------------------------------------------------------------
+Description of the case
+------------------------------------------------------------------
+
+This example is part of the research about integrating geomechanical effects into the Embedded Discrete Fracture Model (EDFM). Here, the EDFM formulation is coupled with analytical poromechanics under an Oedometric Stress Path assumption (flat reservoir with low contrast between it and adjacent formation stiffness), enabling a direct relationship between pore-pressure variations and changes in effective normal stress. These stresses, in turn, govern aperture evolution via the Barton–Bandis constitutive law. Also, we assume:
+
+ - In-situ stress data are available
+ - Andersonian stress regime (no lateral stress, constant vertical total stress)
+ - Permeability computed by Parallel Plates model
+
+The simulation presented below is a simplified case designed to validate the implementation against the analytical solution detailed in the Section ..... The simplification include the use of a single-phase flow solver with TPFA as discretization method, an impermeable matrix, and a through-going fracture.
+
+------------------------------
+Mesh and geometry
+------------------------------
+
+The hexahedral mesh is generated internally and consists of 11 cells aligned along the Y-axis, each measuring 1x1x1 m. The fracture is represented as a plane that crosses the entire domain, dividing it into two halves. A pressure gradient is applied at the boundary cells of the domain and directly at the fracture tips.
+
+.. figure:: stressPathDrivenExampleParaview.png
+ :align: center
+ :width: 500
+ :figclass: align-center
+
+ Visualization of the first simulation step.
+
+------------------------------------------------------------------
+Initial Conditions (Input file setup)
+------------------------------------------------------------------
+
+This example needs the following input values to work properly:
+
+- Surface Element Region's default aperture and Stress-path driven reference aperture must be the same.
+- Initial pressure in the rock and fracture must be the same.
+- Stress-path driven's reference pressure must equal the initial pressure.
+
+
+------------------------------
+Initial Conditions
+------------------------------
+
+This example needs the following input values to work properly:
+
+- Surface Element Region's default aperture and Stress-path driven reference aperture must be the same.
+- Initial pressure in the rock and fracture must be the same.
+- Stress-path driven reference pressure must equal the initial pressure.
+
+
+------------------------------
+Constitutive law
+------------------------------
+
+The new class ``BartonBandisPermeability`` follows a similar implementation to the class ``ParallelPlatesPermeability``. Both are defined as permeability model managed by ``CompressibleSolid`` class, and each contains a private kernel responsible for computing the updated hydraulic aperture. However, because the proposed formulation of the Barton–Bandis model is stress-dependent, the implementation of the new class’s kernel differs and requires the following inputs:
+
+.. literalinclude:: ../../../../../inputFiles/stressPathDrivenGeomechanics/fractureMatrixFlow_edfm_BBP_base.xml
+ :language: xml
+ :start-after:
+ :end-before:
+
+In summary, by assuming the oedometric conditions, it is possible to relate pore pressure and total stresses by
+
+.. math::
+ \Delta \sigma_V = \alpha\, \Delta p, \qquad \Delta \sigma_H = \Delta \sigma_h = \frac{\nu}{1 - \nu}\, \Delta \sigma_V
+
+where :math:`\alpha` is the Biot coefficient and :math:`\nu` is Poisson's ratio.
+
+Given the in-situ effective stress,
+
+.. math::
+ \sigma_0 = \sigma_{\text{ref}} - \alpha\, p_0\, I,
+
+where :math:`p_0` is the reference pore pressure, we compute the normal stress acting on the fracture plane by projecting the updated total stress tensor onto the fracture normal direction:
+
+.. math::
+ \sigma_n = \left( \sigma_0 - \sigma_T \right) : (\vec{n} \otimes \vec{n}),
+
+where :math:`\vec{n}` is the unit normal vector to the fracture plane.
+
+Once the normal effective stress is obtained, the fracture closure is computed using the Barton--Bandis hyperbolic law, which relates normal stress to updated mechanical aperture. The fracture closure is given by
+
+.. math::
+ g_n(\sigma_n) = \frac{\sigma_n\, V_m}{K_{ni}\, V_m + \sigma_n},
+
+where :math:`K_{ni}` is the initial normal stiffness and :math:`V_m` is the maximum possible closure.
+The resulting closure is then used to update the aperture (:math:`a`), from which the intrinsic fracture permeability (:math:`K_f`) is obtained using the parallel-plates law
+
+.. math::
+ K_f = \frac{a^2}{12}.
+
+Both quantities — the updated aperture and the resulting permeability — affect the EDFM transmissibilities.
+
+------------------------------------------------------------------
+Single-phase flow solver with fracture aperture update
+------------------------------------------------------------------
+
+The single-phase flow solver definition in the input file remains the same, with the exception of one optional flag. The flag ``computePrescribedStressPath`` enables the computation of the updated fracture aperture using the new ``BartonBandisPermeability`` class. If this constitutive law is not defined in the input file while the flag is set to 1, an error will be raised. On the other hand, if this flag is set to 0, the fracture aperture will not be updated, even if the constitutive law is defined.
+
+.. literalinclude:: ../../../../../inputFiles/stressPathDrivenGeomechanics/fractureMatrixFlow_edfm_BBP_base.xml
+ :language: xml
+ :start-after:
+ :end-before:
+
+The new class integration is illustrated in the diagram below. Since the proposed permeability class needs the current pressure and fracture normal vector for the effective stress computation, the function's signature in ``CompressibleSolid`` and in ``PermeabilityBase`` class were modified to reflect this requirement.
+
+
+.. figure:: bartonBandisPermeabilityDiagram.png
+ :align: center
+ :width: 800
+ :figclass: align-center
+
+ Diagram of the new implementation to incorporate Barton-Bandis model to flow solvers.
+
+------------------------
+Events
+------------------------
+
+The code below is to save the pressure, aperture and permeability of the fracture at each time step for later validation against the analytical solution.
+
+.. literalinclude:: ../../../../../inputFiles/stressPathDrivenGeomechanics/fractureMatrixFlow_edfm_horizontalFrac_BBP_smoke.xml
+ :language: xml
+ :start-after:
+ :end-before:
+
+------------------------------
+Unit test
+------------------------------
+
+The normal effective stress and aperture update computation can be tested by providing the in-situ (reference) values of pressure and aperture. The resulting aperture should equal the in-situ aperture.
+
+
+------------------------------
+Integration test
+------------------------------
+
+...
+
+------------------------------------------------------------------
+To go further
+------------------------------------------------------------------
+
+
+**Feedback on this example**
+
+For any feedback on this example, please submit a `GitHub issue on the project's GitHub page `_.
diff --git a/src/docs/sphinx/basicExamples/stressPathDriven/bartonBandisPermeabilityDiagram.png b/src/docs/sphinx/basicExamples/stressPathDriven/bartonBandisPermeabilityDiagram.png
new file mode 100644
index 00000000000..fb4d15b279e
Binary files /dev/null and b/src/docs/sphinx/basicExamples/stressPathDriven/bartonBandisPermeabilityDiagram.png differ
diff --git a/src/docs/sphinx/basicExamples/stressPathDriven/bartonBandisSPDDiagram.png b/src/docs/sphinx/basicExamples/stressPathDriven/bartonBandisSPDDiagram.png
new file mode 100644
index 00000000000..435eb647117
Binary files /dev/null and b/src/docs/sphinx/basicExamples/stressPathDriven/bartonBandisSPDDiagram.png differ
diff --git a/src/docs/sphinx/basicExamples/stressPathDriven/stressPathDrivenExampleParaview.png b/src/docs/sphinx/basicExamples/stressPathDriven/stressPathDrivenExampleParaview.png
new file mode 100644
index 00000000000..5da55276b52
Binary files /dev/null and b/src/docs/sphinx/basicExamples/stressPathDriven/stressPathDrivenExampleParaview.png differ
diff --git a/src/docs/sphinx/basicExamples/stressPathDriven/stressPathDrivenExampleParaviewShortFrac.png b/src/docs/sphinx/basicExamples/stressPathDriven/stressPathDrivenExampleParaviewShortFrac.png
new file mode 100644
index 00000000000..45037257521
Binary files /dev/null and b/src/docs/sphinx/basicExamples/stressPathDriven/stressPathDrivenExampleParaviewShortFrac.png differ