diff --git a/examples/functions/function_examples/function_examples.xml b/examples/functions/function_examples/function_examples.xml
index 0cc84199b48..a5b2958eaa6 100644
--- a/examples/functions/function_examples/function_examples.xml
+++ b/examples/functions/function_examples/function_examples.xml
@@ -5,7 +5,7 @@
-
@@ -13,7 +13,7 @@
+ expression="sqrt((x*x)+(y*y)+(z*z))*t"/>-->
diff --git a/host-configs/ubuntu.cmake b/host-configs/ubuntu.cmake
new file mode 100644
index 00000000000..73542768822
--- /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_SPD_base.xml b/inputFiles/stressPathDrivenGeomechanics/fractureMatrixFlow_edfm_SPD_base.xml
new file mode 100644
index 00000000000..6fe717ff2ce
--- /dev/null
+++ b/inputFiles/stressPathDrivenGeomechanics/fractureMatrixFlow_edfm_SPD_base.xml
@@ -0,0 +1,163 @@
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
diff --git a/inputFiles/stressPathDrivenGeomechanics/fractureMatrixFlow_edfm_base.xml b/inputFiles/stressPathDrivenGeomechanics/fractureMatrixFlow_edfm_base.xml
new file mode 100644
index 00000000000..82138056ec3
--- /dev/null
+++ b/inputFiles/stressPathDrivenGeomechanics/fractureMatrixFlow_edfm_base.xml
@@ -0,0 +1,139 @@
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
diff --git a/inputFiles/stressPathDrivenGeomechanics/fractureMatrixFlow_edfm_horizontalFrac_SPD_smoke.xml b/inputFiles/stressPathDrivenGeomechanics/fractureMatrixFlow_edfm_horizontalFrac_SPD_smoke.xml
new file mode 100644
index 00000000000..fc7f8d523aa
--- /dev/null
+++ b/inputFiles/stressPathDrivenGeomechanics/fractureMatrixFlow_edfm_horizontalFrac_SPD_smoke.xml
@@ -0,0 +1,121 @@
+
+
+
+
+
+
+
+
+
+
+
+
+
+ cb1
+
+
+ { -1, 0, 0 }
+
+
+ 1e-8
+
+
+
+
+
+ cb1
+
+
+ { 1, 0, 0 }
+
+
+ 1e-8
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
diff --git a/inputFiles/stressPathDrivenGeomechanics/fractureMatrixFlow_edfm_horizontalFrac_smoke.xml b/inputFiles/stressPathDrivenGeomechanics/fractureMatrixFlow_edfm_horizontalFrac_smoke.xml
new file mode 100644
index 00000000000..3f38b974830
--- /dev/null
+++ b/inputFiles/stressPathDrivenGeomechanics/fractureMatrixFlow_edfm_horizontalFrac_smoke.xml
@@ -0,0 +1,111 @@
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
diff --git a/src/coreComponents/constitutive/CMakeLists.txt b/src/coreComponents/constitutive/CMakeLists.txt
index c835c57966f..7ad8c16e274 100644
--- a/src/coreComponents/constitutive/CMakeLists.txt
+++ b/src/coreComponents/constitutive/CMakeLists.txt
@@ -37,6 +37,7 @@ set( constitutive_headers
capillaryPressure/CapillaryPressureSelector.hpp
capillaryPressure/Layouts.hpp
contact/BartonBandis.hpp
+ contact/BartonBandisStressPathDriven.hpp
contact/CoulombFriction.hpp
contact/FrictionSelector.hpp
contact/FrictionBase.hpp
@@ -232,6 +233,7 @@ set( constitutive_sources
capillaryPressure/TableCapillaryPressureHelpers.cpp
capillaryPressure/VanGenuchtenCapillaryPressure.cpp
contact/BartonBandis.cpp
+ contact/BartonBandisStressPathDriven.cpp
contact/CoulombFriction.cpp
contact/FrictionBase.cpp
contact/FrictionlessContact.cpp
diff --git a/src/coreComponents/constitutive/contact/BartonBandisStressPathDriven.cpp b/src/coreComponents/constitutive/contact/BartonBandisStressPathDriven.cpp
new file mode 100644
index 00000000000..a82dffce671
--- /dev/null
+++ b/src/coreComponents/constitutive/contact/BartonBandisStressPathDriven.cpp
@@ -0,0 +1,86 @@
+/*
+ * ------------------------------------------------------------------------------------------------------------
+ * 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 BartonBandis.cpp
+ */
+
+
+#include "BartonBandisStressPathDriven.hpp"
+
+
+namespace geos
+{
+
+namespace constitutive
+{
+
+using namespace dataRepository;
+
+BartonBandisStressPathDriven::BartonBandisStressPathDriven( string const & name, Group * const parent ):
+ HydraulicApertureBase( name, parent ),
+ m_referenceNormalStress( 0.0 )
+{
+ setInputFlags( InputFlags::REQUIRED );
+
+ registerWrapper( viewKeyStruct::referenceNormalStressString(), &m_referenceNormalStress ).
+ setInputFlag( InputFlags::REQUIRED ).
+ setDescription( "Reference normal stress." );
+ registerWrapper( viewKeyStruct::biotString(), &m_biot ).
+ setApplyDefaultValue( 1.0 ).
+ setInputFlag( InputFlags::REQUIRED ).
+ setDescription( "Biot coefficient." );
+ registerWrapper( viewKeyStruct::poissonString(), &m_poisson ).
+ 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 BartonBandisStressPathDriven::postInputInitialization()
+{
+ GEOS_THROW_IF( m_referenceNormalStress <= 0.0,
+ getFullName() << ": The provided reference stress is zero or negative. Value: " << m_referenceNormalStress,
+ InputError );
+
+}
+
+BartonBandisStressPathDrivenUpdates BartonBandisStressPathDriven::createKernelWrapper() const
+{
+ return KernelWrapper( m_aperture0,
+ m_referenceNormalStress,
+ m_biot,
+ m_poisson,
+ m_normalStiffness,
+ m_referencePressure,
+ m_referenceTotalStress);
+}
+
+REGISTER_CATALOG_ENTRY( ConstitutiveBase, BartonBandisStressPathDriven, string const &, Group * const )
+
+} /* namespace constitutive */
+
+} /* namespace geos */
diff --git a/src/coreComponents/constitutive/contact/BartonBandisStressPathDriven.hpp b/src/coreComponents/constitutive/contact/BartonBandisStressPathDriven.hpp
new file mode 100644
index 00000000000..9bee8095abc
--- /dev/null
+++ b/src/coreComponents/constitutive/contact/BartonBandisStressPathDriven.hpp
@@ -0,0 +1,221 @@
+/*
+ * ------------------------------------------------------------------------------------------------------------
+ * 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 BartonBandis.hpp
+ */
+
+#ifndef GEOS_CONSTITUTIVE_CONTACT_BARTONBANDISSTRESSPATHDRIVEN_HPP_
+#define GEOS_CONSTITUTIVE_CONTACT_BARTONBANDISSTRESSPATHDRIVEN_HPP_
+
+#include "constitutive/contact/HydraulicApertureBase.hpp"
+#include "functions/TableFunction.hpp"
+#include "physicsSolvers/solidMechanics/contact/FractureState.hpp"
+
+namespace geos
+{
+
+namespace constitutive
+{
+
+/**
+ * @class BartonBandisUpdates
+ *
+ * This class is used for in-kernel contact relation updates
+ */
+class BartonBandisStressPathDrivenUpdates
+{
+public:
+
+ BartonBandisStressPathDrivenUpdates( real64 const aperture0,
+ real64 const referenceNormalStress,
+ real64 const biot,
+ real64 const poisson,
+ real64 const normalStiffness,
+ real64 const referencePressure,
+ R1Tensor const &referenceTotalStress)
+ : m_aperture0( aperture0 ),
+ m_referenceNormalStress( referenceNormalStress ),
+ m_biot( biot ),
+ m_poisson( poisson ),
+ m_normalStiffness( normalStiffness ), // Kni
+ m_referencePressure( referencePressure ),
+ m_referenceTotalStress( referenceTotalStress )
+ {
+ m_referenceEffectiveStress[0] = m_referenceTotalStress[0] - m_biot*m_referencePressure;
+ m_referenceEffectiveStress[1] = m_referenceTotalStress[1] - m_biot*m_referencePressure;
+ m_referenceEffectiveStress[2] = m_referenceTotalStress[2] - m_biot*m_referencePressure;
+ }
+
+ /// Default copy constructor
+ BartonBandisStressPathDrivenUpdates( BartonBandisStressPathDrivenUpdates const & ) = default;
+
+ /// Default move constructor
+ BartonBandisStressPathDrivenUpdates( BartonBandisStressPathDrivenUpdates && ) = default;
+
+ /// Deleted default constructor
+ BartonBandisStressPathDrivenUpdates() = default;
+
+ /// Deleted copy assignment operator
+ BartonBandisStressPathDrivenUpdates & operator=( BartonBandisStressPathDrivenUpdates const & ) = delete;
+
+ /// Deleted move assignment operator
+ BartonBandisStressPathDrivenUpdates & operator=( BartonBandisStressPathDrivenUpdates && ) = delete;
+
+ /**
+ * @brief Evaluate the effective aperture, and its derivative wrt aperture
+ * @param[in] aperture the model aperture/gap
+ * @param[out] dHydraulicAperture_dAperture the derivative of the effective aperture wrt aperture
+ * @return The hydraulic aperture that is always > 0
+ */
+ GEOS_HOST_DEVICE
+ real64 computeHydraulicAperture( real64 const pressure,
+ array1d< real64 > const & normal) const;
+
+
+private:
+ real64 m_aperture0;
+ real64 m_referenceNormalStress;
+
+ real64 m_biot;
+ real64 m_poisson;
+ real64 m_normalStiffness; // Kni
+ real64 m_referencePressure; // p_0
+
+ R1Tensor m_referenceTotalStress; // sigmaT_0 computed analytically
+ R1Tensor m_referenceEffectiveStress; // sigma_0
+
+ real64 computeFractureStress( real64 const pressure, array1d< real64 > const & normal ) const;
+};
+
+
+/**
+ * @class BartonBandis
+ *
+ * This class serves as the interface for implementing contact enforcement constitutive relations.
+ * This does not include the actual enforcement algorithm, but only the constitutive relations that
+ * govern the behavior of the contact. So things like penalty, or friction, or kinematic constraint.
+ */
+class BartonBandisStressPathDriven : public HydraulicApertureBase
+{
+public:
+
+ /**
+ * @brief The standard data repository constructor
+ * @param name The name of the relation in the data repository
+ * @param parent The name of the parent Group that holds this relation object.
+ */
+ BartonBandisStressPathDriven( string const & name,
+ Group * const parent );
+
+ static string catalogName() { return "BartonBandisStressPathDriven"; }
+
+ virtual string getCatalogName() const override { return catalogName(); }
+
+
+ /// Type of kernel wrapper for in-kernel update
+ using KernelWrapper = BartonBandisStressPathDrivenUpdates;
+
+ /**
+ * @brief Create an update kernel wrapper.
+ * @return the wrapper
+ */
+ KernelWrapper createKernelWrapper() const;
+
+ struct viewKeyStruct : public HydraulicApertureBase::viewKeyStruct
+ {
+ static constexpr char const * referenceNormalStressString() { return "referenceNormalStress"; }
+ static constexpr char const * biotString() { return "biot"; }
+ static constexpr char const * poissonString() { return "poisson"; }
+ 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_referenceNormalStress;
+
+ real64 m_biot;
+ real64 m_poisson;
+ real64 m_normalStiffness; // Kni
+ real64 m_referencePressure; // p_0
+
+ R1Tensor m_referenceTotalStress; // sigmaT_0 computed analytically
+};
+
+GEOS_HOST_DEVICE
+GEOS_FORCE_INLINE
+real64 BartonBandisStressPathDrivenUpdates::computeHydraulicAperture( real64 const pressure,
+ array1d< real64 > const & normal ) const
+{
+ real64 const biot_pressure = m_biot * m_referencePressure; // biot is alpha in the equations
+
+ // Computation of maximum fracture closure (Barton-Bandis parameter)
+ // Fracture traction via Terzaghi's Principle
+ real64 const sigma_c0[3] = { m_referenceTotalStress[0] * normal[0] - biot_pressure * normal[0],
+ m_referenceTotalStress[1] * normal[1] - biot_pressure * normal[1],
+ m_referenceTotalStress[2] * normal[2] - biot_pressure * normal[2] };
+ real64 const sigma_n0 = sigma_c0[0]*normal[0] +
+ sigma_c0[1]*normal[1] +
+ sigma_c0[2]*normal[2];
+ real64 const g0 = (-m_normalStiffness*m_aperture0 +
+ std::sqrt((m_normalStiffness*m_aperture0)*
+ (m_normalStiffness*m_aperture0) +
+ 4.0*m_normalStiffness*sigma_n0*m_aperture0)) / (2.0*m_normalStiffness);
+ real64 const maximumFractureClosure = g0 + m_aperture0; // Vm
+
+ // Normal effective stress on the fracture
+ real64 const sigmaN_N = computeFractureStress( pressure, normal );
+ real64 const fractureClosure = sigmaN_N*maximumFractureClosure/(m_normalStiffness*maximumFractureClosure + sigmaN_N); // gn_BB
+
+ // Compute the new aperture which is equal to the aperture at the free-stress state
+ // minus the closure from the free-stress state to the current state
+ real64 const newHydraulicAperture = maximumFractureClosure - fractureClosure;
+
+ return newHydraulicAperture;
+}
+
+GEOS_HOST_DEVICE
+GEOS_FORCE_INLINE
+real64 BartonBandisStressPathDrivenUpdates::computeFractureStress( real64 const pressure,
+ array1d< real64 > const & normal ) const
+{
+ auto matmul = [](real64 const (&u)[3], array1d< real64 > const &v, array1d< real64 > &r) -> void
+ {
+ r[0] = u[0]*v[0];
+ r[1] = u[1]*v[1];
+ r[2] = u[2]*v[2];
+ };
+
+ real64 const deltaSigmaZ = m_biot * (pressure - m_referencePressure);
+ real64 const poisson_deltaSigma = deltaSigmaZ * m_poisson/(1.0 - m_poisson);
+ // sigma: matrix diagonal
+ real64 effectiveStress[3] = { m_referenceEffectiveStress[0] - poisson_deltaSigma,
+ m_referenceEffectiveStress[1] - poisson_deltaSigma,
+ m_referenceEffectiveStress[2] - deltaSigmaZ };
+ array1d< real64 > effectiveStressOnFracture(3); // sigma_c
+ matmul(effectiveStress, normal, effectiveStressOnFracture);
+ real64 normalComponentOfStressOnFracture = dot(effectiveStressOnFracture, normal); // sigmaN_N
+ return normalComponentOfStressOnFracture;
+}
+} /* namespace constitutive */
+
+} /* namespace geos */
+
+#endif /* GEOS_CONSTITUTIVE_CONTACT_BARTONBANDISSTRESSPATHDRIVEN_HPP_ */
diff --git a/src/coreComponents/constitutive/contact/HydraulicApertureRelationSelector.hpp b/src/coreComponents/constitutive/contact/HydraulicApertureRelationSelector.hpp
index f8675f1dc1e..5e1930fa7ea 100644
--- a/src/coreComponents/constitutive/contact/HydraulicApertureRelationSelector.hpp
+++ b/src/coreComponents/constitutive/contact/HydraulicApertureRelationSelector.hpp
@@ -23,6 +23,7 @@
#include "constitutive/ConstitutivePassThruHandler.hpp"
#include "constitutive/contact/HydraulicApertureTable.hpp"
#include "constitutive/contact/BartonBandis.hpp"
+#include "constitutive/contact/BartonBandisStressPathDriven.hpp"
namespace geos
{
@@ -46,6 +47,19 @@ void constitutiveUpdatePassThru( HydraulicApertureBase & contact,
BartonBandis >::execute( contact, std::forward< LAMBDA >( lambda ) );
}
+template< typename LAMBDA >
+void constitutiveUpdatePassThru( BartonBandisStressPathDriven const & contact,
+ LAMBDA && lambda )
+{
+ ConstitutivePassThruHandler< BartonBandisStressPathDriven >::execute( contact, std::forward< LAMBDA >( lambda ) );
+}
+
+template< typename LAMBDA >
+void constitutiveUpdatePassThru( BartonBandisStressPathDriven & contact,
+ LAMBDA && lambda )
+{
+ ConstitutivePassThruHandler< BartonBandisStressPathDriven >::execute( contact, std::forward< LAMBDA >( lambda ) );
+}
} /* namespace constitutive */
} /* namespace geos */
diff --git a/src/coreComponents/constitutive/permeability/PressurePermeability.hpp b/src/coreComponents/constitutive/permeability/PressurePermeability.hpp
index 9d654048a0e..6e4ed374512 100644
--- a/src/coreComponents/constitutive/permeability/PressurePermeability.hpp
+++ b/src/coreComponents/constitutive/permeability/PressurePermeability.hpp
@@ -87,7 +87,7 @@ class PressurePermeabilityUpdate : public PermeabilityBaseUpdate
referencePermeability[0] = m_referencePermeability[k][0][0];
referencePermeability[1] = m_referencePermeability[k][0][1];
referencePermeability[2] = m_referencePermeability[k][0][2];
-
+
switch( m_presModelType )
{
case PressureModelType::Exponential:
diff --git a/src/coreComponents/physicsSolvers/fluidFlow/FlowSolverBase.cpp b/src/coreComponents/physicsSolvers/fluidFlow/FlowSolverBase.cpp
index 8a0474bf79c..2eb4de61b06 100644
--- a/src/coreComponents/physicsSolvers/fluidFlow/FlowSolverBase.cpp
+++ b/src/coreComponents/physicsSolvers/fluidFlow/FlowSolverBase.cpp
@@ -36,6 +36,8 @@
#include "physicsSolvers/fluidFlow/kernels/MinPoreVolumeMaxPorosityKernel.hpp"
#include "physicsSolvers/fluidFlow/kernels/StencilWeightsUpdateKernel.hpp"
+#include "constitutive/contact/HydraulicApertureRelationSelector.hpp"
+
namespace geos
{
@@ -109,6 +111,7 @@ FlowSolverBase::FlowSolverBase( string const & name,
PhysicsSolverBase( name, parent ),
m_numDofPerCell( 0 ),
m_isThermal( 0 ),
+ m_computePrescribedStressPath( 0 ),
m_keepVariablesConstantDuringInitStep( false ),
m_isFixedStressPoromechanicsUpdate( false ),
m_isJumpStabilized( false ),
@@ -119,6 +122,11 @@ FlowSolverBase::FlowSolverBase( string const & name,
setInputFlag( InputFlags::OPTIONAL ).
setDescription( "Flag indicating whether the problem is thermal or not." );
+ this->registerWrapper( viewKeyStruct::computesPrescribedStressPathString(), &m_computePrescribedStressPath ).
+ setApplyDefaultValue( 0 ).
+ setInputFlag( InputFlags::OPTIONAL ).
+ setDescription( "Flag to determine whether or not this simulation computes the precribed stress path." );
+
this->registerWrapper( viewKeyStruct::allowNegativePressureString(), &m_allowNegativePressure ).
setApplyDefaultValue( 0 ). // negative pressure is not allowed by default
setInputFlag( InputFlags::OPTIONAL ).
@@ -625,6 +633,7 @@ void FlowSolverBase::updatePorosityAndPermeability( SurfaceElementSubRegion & su
{
GEOS_MARK_FUNCTION;
+
arrayView1d< real64 const > const & pressure = subRegion.getField< flow::pressure >();
arrayView1d< real64 const > const newHydraulicAperture = subRegion.getField< flow::hydraulicAperture >();
@@ -642,6 +651,42 @@ void FlowSolverBase::updatePorosityAndPermeability( SurfaceElementSubRegion & su
} );
}
+void FlowSolverBase::updateHydraulicAperture( SurfaceElementSubRegion & subRegion ) const
+{
+ GEOS_MARK_FUNCTION;
+
+ arrayView1d< real64 const > const & pressure = subRegion.getField< fields::flow::pressure >();
+ arrayView1d< real64 > const & newHydraulicAperture = subRegion.getField< fields::flow::hydraulicAperture >();
+ arrayView2d< real64 const > const & normalVector = subRegion.getField< fields::normalVector >(); // mesh/MeshFields.hpp
+
+ string const & hydraulicApertureRelationName =
+ subRegion.template getReference< string >( viewKeyStruct::hydraulicApertureRelationNameString() );
+ BartonBandisStressPathDriven const & hydraulicApertureModel =
+ this->template getConstitutiveModel< BartonBandisStressPathDriven >( subRegion, hydraulicApertureRelationName );
+
+ BartonBandisStressPathDrivenUpdates hydraulicApertureWrapper = hydraulicApertureModel.createKernelWrapper();
+
+ real64 sumAperture = 0.0;
+ 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];
+
+ newHydraulicAperture[k] = hydraulicApertureWrapper.computeHydraulicAperture( pressure[k], normal );
+
+ sumAperture += newHydraulicAperture[k];
+ } );
+
+ real64 const averageAperture = sumAperture / subRegion.size();
+
+ forAll< parallelDevicePolicy<> >( subRegion.size(), [&newHydraulicAperture, averageAperture] GEOS_DEVICE ( localIndex const k )
+ {
+ newHydraulicAperture[k] = averageAperture;
+ } );
+}
+
void FlowSolverBase::findMinMaxElevationInEquilibriumTarget( DomainPartition & domain, // cannot be const...
std::map< string, localIndex > const & equilNameToEquilId,
diff --git a/src/coreComponents/physicsSolvers/fluidFlow/FlowSolverBase.hpp b/src/coreComponents/physicsSolvers/fluidFlow/FlowSolverBase.hpp
index 110c27f293e..453f050eb13 100644
--- a/src/coreComponents/physicsSolvers/fluidFlow/FlowSolverBase.hpp
+++ b/src/coreComponents/physicsSolvers/fluidFlow/FlowSolverBase.hpp
@@ -82,9 +82,13 @@ class FlowSolverBase : public PhysicsSolverBase
static constexpr char const * fluidNamesString() { return "fluidNames"; }
static constexpr char const * solidNamesString() { return "solidNames"; }
+ static constexpr char const * stressPathNamesString() { return "stressPathNames"; }
static constexpr char const * permeabilityNamesString() { return "permeabilityNames"; }
static constexpr char const * solidInternalEnergyNamesString() { return "solidInternalEnergyNames"; }
static constexpr char const * thermalConductivityNamesString() { return "thermalConductivityNames"; }
+
+ static constexpr char const * computesPrescribedStressPathString() { return "computesPrescribedStressPath"; }
+ static constexpr char const * hydraulicApertureRelationNameString() { return "hydraulicApertureRelationName"; }
};
/**
@@ -109,6 +113,8 @@ class FlowSolverBase : public PhysicsSolverBase
virtual void updatePorosityAndPermeability( SurfaceElementSubRegion & subRegion ) const;
+ void updateHydraulicAperture( SurfaceElementSubRegion & subRegion ) const;
+
/**
* @brief Utility function to save the iteration state (useful for sequential simulations)
* @param[in] domain the domain partition
@@ -258,6 +264,9 @@ class FlowSolverBase : public PhysicsSolverBase
/// flag to determine whether or not this is a thermal simulation
integer m_isThermal;
+ /// flag to determine whether or not this simulation computes the precribed stress path
+ integer m_computePrescribedStressPath;
+
/// the input temperature
real64 m_inputTemperature;
diff --git a/src/coreComponents/physicsSolvers/fluidFlow/SinglePhaseBase.cpp b/src/coreComponents/physicsSolvers/fluidFlow/SinglePhaseBase.cpp
index bb4d66c6122..f636ee2239c 100644
--- a/src/coreComponents/physicsSolvers/fluidFlow/SinglePhaseBase.cpp
+++ b/src/coreComponents/physicsSolvers/fluidFlow/SinglePhaseBase.cpp
@@ -47,6 +47,7 @@
#include "physicsSolvers/fluidFlow/kernels/singlePhase/FluidUpdateKernel.hpp"
#include "physicsSolvers/fluidFlow/kernels/singlePhase/SolidInternalEnergyUpdateKernel.hpp"
+#include "constitutive/contact/HydraulicApertureRelationSelector.hpp"
namespace geos
{
@@ -131,6 +132,15 @@ void SinglePhaseBase::setConstitutiveNames( ElementSubRegionBase & subRegion ) c
{
setConstitutiveName< SinglePhaseThermalConductivityBase >( subRegion, viewKeyStruct::thermalConductivityNamesString(), "singlephase thermal conductivity" );
}
+
+ if( m_computePrescribedStressPath )
+ {
+ // Only for fractures
+ if( dynamic_cast< SurfaceElementSubRegion * >( &subRegion ) )
+ {
+ setConstitutiveName< constitutive::BartonBandisStressPathDriven >( subRegion, viewKeyStruct::hydraulicApertureRelationNameString(), "hydraulic aperture" );
+ }
+ }
}
void SinglePhaseBase::initializeAquiferBC() const
@@ -1184,6 +1194,20 @@ void SinglePhaseBase::updateState( DomainPartition & domain )
{
GEOS_MARK_FUNCTION;
+ if(m_computePrescribedStressPath)
+ {
+ forDiscretizationOnMeshTargets( domain.getMeshBodies(), [&]( string const &,
+ MeshLevel & mesh,
+ string_array const & regionNames )
+ {
+ mesh.getElemManager().forElementSubRegions< SurfaceElementSubRegion >( regionNames, [&]( localIndex const,
+ auto & subRegion )
+ {
+ updateHydraulicAperture( subRegion );
+ } );
+ } );
+ }
+
forDiscretizationOnMeshTargets( domain.getMeshBodies(), [&]( string const &,
MeshLevel & mesh,
string_array const & regionNames )
diff --git a/src/coreComponents/physicsSolvers/fluidFlow/SinglePhaseBase.hpp b/src/coreComponents/physicsSolvers/fluidFlow/SinglePhaseBase.hpp
index a9e0af8e961..e72bf051256 100644
--- a/src/coreComponents/physicsSolvers/fluidFlow/SinglePhaseBase.hpp
+++ b/src/coreComponents/physicsSolvers/fluidFlow/SinglePhaseBase.hpp
@@ -266,7 +266,7 @@ class SinglePhaseBase : public FlowSolverBase
arrayView1d< real64 > const & localRhs ) const = 0;
virtual void
- updateState ( DomainPartition & domain ) override final;
+ updateState ( DomainPartition & domain ) override;
/**
* @brief Function to update all constitutive state and dependent variables
diff --git a/src/coreComponents/physicsSolvers/fluidFlow/SinglePhaseStressPathDrivenFVM.cpp b/src/coreComponents/physicsSolvers/fluidFlow/SinglePhaseStressPathDrivenFVM.cpp
new file mode 100644
index 00000000000..7a3a6dc3a24
--- /dev/null
+++ b/src/coreComponents/physicsSolvers/fluidFlow/SinglePhaseStressPathDrivenFVM.cpp
@@ -0,0 +1,127 @@
+/*
+ * ------------------------------------------------------------------------------------------------------------
+ * 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 SinglePhaseFVM.cpp
+ */
+
+#include "SinglePhaseStressPathDrivenFVM.hpp"
+
+#include "physicsSolvers/LogLevelsInfo.hpp"
+#include "constitutive/contact/HydraulicApertureRelationSelector.hpp"
+#include "mesh/DomainPartition.hpp"
+#include "constitutive/ConstitutivePassThru.hpp"
+#include "constitutive/solid/CoupledSolidBase.hpp"
+#include "physicsSolvers/solidMechanics/contact/FractureState.hpp"
+
+/**
+ * @namespace the geos namespace that encapsulates the majority of the code
+ */
+namespace geos
+{
+
+using namespace constitutive;
+
+SinglePhaseStressPathDrivenFVM::SinglePhaseStressPathDrivenFVM( const string & name,
+ dataRepository::Group * const parent ):
+ SinglePhaseFVM< SinglePhaseBase >( name, parent )
+{
+ //SinglePhaseFVM< SinglePhaseBase >::template addLogLevel< logInfo::Convergence >();
+}
+
+void SinglePhaseStressPathDrivenFVM::setConstitutiveNamesCallSuper( ElementSubRegionBase & subRegion ) const
+{
+ FlowSolverBase::setConstitutiveNamesCallSuper( subRegion );
+
+ // To make Barton Bandis constitutive law mandatory
+ if( dynamic_cast< SurfaceElementSubRegion * >( &subRegion ) )
+ {
+ this->template setConstitutiveName< constitutive::BartonBandisStressPathDriven >( subRegion,
+ viewKeyStruct::hydraulicApertureRelationNameString(), "hydraulic aperture" );
+ }
+}
+
+void SinglePhaseStressPathDrivenFVM::updateState( DomainPartition & domain )
+{
+ GEOS_MARK_FUNCTION;
+ // As in PoromechanicsConformingFractures.hpp
+
+ // remove the contribution of the hydraulic aperture from the stencil weights
+ prepareStencilWeights( domain );
+
+ updateHydraulicAperture( domain );
+
+ // update the stencil weights using the updated hydraulic aperture
+ updateStencilWeights( domain );
+
+ SinglePhaseBase::updateState( domain ); // updates the permeability of the cell and the fracture
+}
+
+void SinglePhaseStressPathDrivenFVM::updateHydraulicAperture( DomainPartition & domain ) const
+{
+ GEOS_MARK_FUNCTION;
+
+ forDiscretizationOnMeshTargets( domain.getMeshBodies(), [&]( string const &,
+ MeshLevel & mesh,
+ string_array const & regionNames )
+ {
+ mesh.getElemManager().forElementSubRegions< SurfaceElementSubRegion >( regionNames, [&]( localIndex const,
+ auto & subRegion )
+ {
+ updateFractureAperture( subRegion );
+ } );
+ } );
+}
+
+void SinglePhaseStressPathDrivenFVM::updateFractureAperture( SurfaceElementSubRegion & subRegion ) const
+{
+ GEOS_MARK_FUNCTION;
+
+ arrayView1d< real64 const > const & pressure = subRegion.getField< fields::flow::pressure >();
+ arrayView1d< real64 > const & newHydraulicAperture = subRegion.getField< fields::flow::hydraulicAperture >();
+ arrayView2d< real64 const > const & normalVector = subRegion.getField< fields::normalVector >(); // mesh/MeshFields.hpp
+
+ string const & hydraulicApertureRelationName =
+ subRegion.template getReference< string >( viewKeyStruct::hydraulicApertureRelationNameString() );
+ BartonBandisStressPathDriven const & hydraulicApertureModel =
+ this->template getConstitutiveModel< BartonBandisStressPathDriven >( subRegion, hydraulicApertureRelationName );
+
+ BartonBandisStressPathDrivenUpdates hydraulicApertureWrapper = hydraulicApertureModel.createKernelWrapper();
+
+ real64 sumAperture = 0.0;
+ 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];
+
+ newHydraulicAperture[k] = hydraulicApertureWrapper.computeHydraulicAperture( pressure[k], normal );
+
+ sumAperture += newHydraulicAperture[k];
+ } );
+
+ real64 const averageAperture = sumAperture / subRegion.size();
+
+ forAll< parallelDevicePolicy<> >( subRegion.size(), [&newHydraulicAperture, averageAperture] GEOS_DEVICE ( localIndex const k )
+ {
+ newHydraulicAperture[k] = averageAperture;
+ } );
+}
+
+
+REGISTER_CATALOG_ENTRY( PhysicsSolverBase, SinglePhaseStressPathDrivenFVM, string const &, dataRepository::Group * const )
+
+} /* namespace geos */
diff --git a/src/coreComponents/physicsSolvers/fluidFlow/SinglePhaseStressPathDrivenFVM.hpp b/src/coreComponents/physicsSolvers/fluidFlow/SinglePhaseStressPathDrivenFVM.hpp
new file mode 100644
index 00000000000..e771d6b123a
--- /dev/null
+++ b/src/coreComponents/physicsSolvers/fluidFlow/SinglePhaseStressPathDrivenFVM.hpp
@@ -0,0 +1,107 @@
+/*
+ * ------------------------------------------------------------------------------------------------------------
+ * 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 SinglePhaseStressPathDrivenFVM.hpp
+ */
+
+#ifndef GEOS_PHYSICSSOLVERS_FLUIDFLOW_SINGLEPHASESTRESSPATHDRIVENFVM_HPP_
+#define GEOS_PHYSICSSOLVERS_FLUIDFLOW_SINGLEPHASESTRESSPATHDRIVENFVM_HPP_
+
+#include "physicsSolvers/fluidFlow/SinglePhaseFVM.hpp"
+
+namespace geos
+{
+
+
+/**
+ * @class SinglePhaseStressPathDrivenFVM
+ *
+ * class to perform a single phase finite volume solve
+ * using only cell-centered variables
+ * works with both TPFA and MPFA
+ */
+class SinglePhaseStressPathDrivenFVM : public SinglePhaseFVM< SinglePhaseBase >
+{
+public:
+
+ /**
+ * @brief main constructor for Group Objects
+ * @param name the name of this instantiation of Group in the repository
+ * @param parent the parent group of this instantiation of Group
+ */
+ SinglePhaseStressPathDrivenFVM( const string & name,
+ dataRepository::Group * const parent );
+
+
+ /// deleted default constructor
+ SinglePhaseStressPathDrivenFVM() = delete;
+
+ /// deleted copy constructor
+ SinglePhaseStressPathDrivenFVM( SinglePhaseStressPathDrivenFVM const & ) = delete;
+
+ /// default move constructor
+ SinglePhaseStressPathDrivenFVM( SinglePhaseStressPathDrivenFVM && ) = default;
+
+ /// deleted assignment operator
+ SinglePhaseStressPathDrivenFVM & operator=( SinglePhaseStressPathDrivenFVM const & ) = delete;
+
+ /// deleted move operator
+ SinglePhaseStressPathDrivenFVM & operator=( SinglePhaseStressPathDrivenFVM && ) = delete;
+
+ /**
+ * @brief default destructor
+ */
+ virtual ~SinglePhaseStressPathDrivenFVM() override = default;
+
+ /**
+ * @brief name of the node manager in the object catalog
+ * @return string that contains the catalog name to generate a new NodeManager object through the object catalog.
+ */
+ static string catalogName()
+ {
+ return "SinglePhaseStressPathDrivenFVM";
+ }
+
+ /**
+ * @copydoc PhysicsSolverBase::getCatalogName()
+ */
+ string getCatalogName() const override { return catalogName(); }
+
+ struct viewKeyStruct : SinglePhaseBase::viewKeyStruct
+ {
+ /// Name of the hydraulicApertureRelationName
+ static constexpr char const * hydraulicApertureRelationNameString() { return "hydraulicApertureRelationName"; }
+ };
+
+ // The same as in PoromechanicsSolver.hpp
+ virtual void setConstitutiveNamesCallSuper( ElementSubRegionBase & subRegion ) const override;
+
+ virtual void updateState ( DomainPartition & domain ) override final;
+
+
+private:
+
+ void updateHydraulicAperture( DomainPartition & domain ) const;
+
+ void updateFractureAperture( SurfaceElementSubRegion & subRegion ) const;
+
+ real64 m_apertureInSitu;
+
+};
+
+} /* namespace geos */
+
+#endif //GEOS_PHYSICSSOLVERS_FLUIDFLOW_SINGLEPHASESTRESSPATHDRIVENFVM_HPP_
diff --git a/src/coreComponents/schema/schema.xsd b/src/coreComponents/schema/schema.xsd
index b4480011c41..b9655ab7bf9 100644
--- a/src/coreComponents/schema/schema.xsd
+++ b/src/coreComponents/schema/schema.xsd
@@ -222,18 +222,10 @@
-
-
-
-
-
-
-
-
@@ -835,42 +827,82 @@
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
@@ -1433,26 +1465,10 @@ stress - traction is applied to the faces as specified by the inner product of i
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
@@ -1461,18 +1477,6 @@ stress - traction is applied to the faces as specified by the inner product of i
-
-
-
-
-
-
-
-
-
-
-
-
@@ -3058,6 +3062,8 @@ When set to 2 also output convergence information to a csv-->
+
+
@@ -3186,6 +3192,8 @@ When set to 2 also output convergence information to a csv-->
+
+
@@ -3927,6 +3935,8 @@ When set to 2 also output convergence information to a csv-->
+
+
@@ -4444,6 +4454,8 @@ When set to 2 also output convergence information to a csv-->
+
+
@@ -4509,6 +4521,8 @@ When set to 2 also output convergence information to a csv-->
+
+
@@ -4615,6 +4629,8 @@ When set to 2 also output convergence information to a csv-->
+
+
@@ -4669,6 +4685,8 @@ When set to 2 also output convergence information to a csv-->
+
+
@@ -4962,6 +4980,8 @@ When set to 2 also output convergence information to a csv-->
+
+
@@ -6024,6 +6044,7 @@ Information output from lower logLevels is added with the desired log level
+
@@ -6081,15 +6102,25 @@ Information output from lower logLevels is added with the desired log level
+
+
+
+
+
+
+
+
+
+
@@ -6125,6 +6156,22 @@ Information output from lower logLevels is added with the desired log level
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
@@ -7315,6 +7362,18 @@ If you want to do a three-phase simulation, please use instead wettingIntermedia
+
+
+
+
+
+
+
+
+
+
+
+
@@ -7327,6 +7386,18 @@ If you want to do a three-phase simulation, please use instead wettingIntermedia
+
+
+
+
+
+
+
+
+
+
+
+
@@ -7339,6 +7410,18 @@ If you want to do a three-phase simulation, please use instead wettingIntermedia
+
+
+
+
+
+
+
+
+
+
+
+
@@ -7351,6 +7434,18 @@ If you want to do a three-phase simulation, please use instead wettingIntermedia
+
+
+
+
+
+
+
+
+
+
+
+
@@ -7363,6 +7458,18 @@ If you want to do a three-phase simulation, please use instead wettingIntermedia
+
+
+
+
+
+
+
+
+
+
+
+
@@ -7375,6 +7482,18 @@ If you want to do a three-phase simulation, please use instead wettingIntermedia
+
+
+
+
+
+
+
+
+
+
+
+
@@ -7387,6 +7506,18 @@ If you want to do a three-phase simulation, please use instead wettingIntermedia
+
+
+
+
+
+
+
+
+
+
+
+
@@ -7399,6 +7530,18 @@ If you want to do a three-phase simulation, please use instead wettingIntermedia
+
+
+
+
+
+
+
+
+
+
+
+
@@ -7411,6 +7554,18 @@ If you want to do a three-phase simulation, please use instead wettingIntermedia
+
+
+
+
+
+
+
+
+
+
+
+
@@ -7423,6 +7578,18 @@ If you want to do a three-phase simulation, please use instead wettingIntermedia
+
+
+
+
+
+
+
+
+
+
+
+
diff --git a/src/coreComponents/schema/schema.xsd.other b/src/coreComponents/schema/schema.xsd.other
index 38ec3095e83..f032ffbf5a4 100644
--- a/src/coreComponents/schema/schema.xsd.other
+++ b/src/coreComponents/schema/schema.xsd.other
@@ -341,15 +341,11 @@
-
-
-
-
@@ -516,7 +512,7 @@
-
+
@@ -1502,13 +1498,14 @@
-
+
+
@@ -1566,15 +1563,25 @@
+
+
+
+
+
+
+
+
+
+
@@ -1603,6 +1610,7 @@
+
@@ -2741,15 +2749,25 @@
+
+
+
+
+
+
+
+
+
+