From 92aa5c1de65d74272f0bcae19ba7f247109a2b04 Mon Sep 17 00:00:00 2001 From: jacques franc Date: Mon, 16 Mar 2026 16:57:24 +0100 Subject: [PATCH 1/8] draft/ enable disp update in SEQ ALM --- .../constitutive/solid/CoupledSolid.hpp | 4 ++ .../fluidFlow/FlowSolverBase.cpp | 58 ++++++++++++++++++- 2 files changed, 61 insertions(+), 1 deletion(-) diff --git a/src/coreComponents/constitutive/solid/CoupledSolid.hpp b/src/coreComponents/constitutive/solid/CoupledSolid.hpp index ddab1b97a19..da74fb42ce8 100644 --- a/src/coreComponents/constitutive/solid/CoupledSolid.hpp +++ b/src/coreComponents/constitutive/solid/CoupledSolid.hpp @@ -160,8 +160,12 @@ class CoupledSolid : public CoupledSolidBase getPermModel() ); } +public: + typedef PERM_TYPE PermType; + //START_SPHINX_INCLUDE_01 protected: + SOLID_TYPE const & getSolidModel() const { return this->getParent().template getGroup< SOLID_TYPE >( m_solidModelName ); } diff --git a/src/coreComponents/physicsSolvers/fluidFlow/FlowSolverBase.cpp b/src/coreComponents/physicsSolvers/fluidFlow/FlowSolverBase.cpp index bbc5a8466c3..6dea8e1506f 100644 --- a/src/coreComponents/physicsSolvers/fluidFlow/FlowSolverBase.cpp +++ b/src/coreComponents/physicsSolvers/fluidFlow/FlowSolverBase.cpp @@ -33,6 +33,12 @@ #include "mesh/DomainPartition.hpp" #include "physicsSolvers/LogLevelsInfo.hpp" #include "physicsSolvers/fluidFlow/FlowSolverBaseFields.hpp" + +#include "physicsSolvers/solidMechanics/contact/ContactFields.hpp" //should not be here -- testing only design flaw +#include "constitutive/permeability/SlipDependentPermeability.hpp" +#include "constitutive/permeability/WillisRichardsPermeability.hpp" +#include "constitutive/permeability/ParallelPlatesPermeability.hpp" + #include "physicsSolvers/fluidFlow/kernels/MinPoreVolumeMaxPorosityKernel.hpp" #include "physicsSolvers/fluidFlow/kernels/StencilWeightsUpdateKernel.hpp" @@ -104,6 +110,35 @@ void updatePorosityAndPermeabilityFromPressureAndAperture( POROUSWRAPPER_TYPE po } ); } +template< typename POROUSWRAPPER_TYPE > +void updatePorosityAndPermeabilityFromPressurApertureJumpAndTraction( POROUSWRAPPER_TYPE porousWrapper, + SurfaceElementSubRegion & subRegion, + arrayView1d< real64 const > const & pressure, + arrayView1d< real64 const > const & oldHydraulicAperture, + arrayView1d< real64 const > const & newHydraulicAperture, + arrayView1d< real64 const > const & dHydraulicAperture_dNormalJump, + arrayView2d< real64 const > const & dispJump, + arrayView2d< real64 const > const & fracTraction) +{ + forAll< parallelDevicePolicy<> >( subRegion.size(), [=] GEOS_DEVICE ( localIndex const k ) + { + for( localIndex q = 0; q < porousWrapper.numGauss(); ++q ) + { + real64 const jump[3] = LVARRAY_TENSOROPS_INIT_LOCAL_3 ( dispJump[k] ); + real64 const traction[3] = LVARRAY_TENSOROPS_INIT_LOCAL_3 ( fracTraction[k] ); + porousWrapper.updateStateFromPressureApertureJumpAndTraction( k, q, + pressure[k], + oldHydraulicAperture[k], + newHydraulicAperture[k], + dHydraulicAperture_dNormalJump[k], + jump, + traction + ); + } + } ); + +} + FlowSolverBase::FlowSolverBase( string const & name, Group * const parent ): PhysicsSolverBase( name, parent ), @@ -632,13 +667,34 @@ void FlowSolverBase::updatePorosityAndPermeability( SurfaceElementSubRegion & su arrayView1d< real64 const > const newHydraulicAperture = subRegion.getField< flow::hydraulicAperture >(); arrayView1d< real64 const > const oldHydraulicAperture = subRegion.getField< flow::aperture0 >(); + arrayView2d< real64 const > const dispJump = subRegion.getField< fields::contact::dispJump >(); + arrayView2d< real64 const > const fractureTraction = subRegion.getField< fields::contact::traction >(); + string const & solidName = subRegion.getReference< string >( viewKeyStruct::solidNamesString() ); CoupledSolidBase & porousSolid = subRegion.getConstitutiveModel< CoupledSolidBase >( solidName ); constitutive::ConstitutivePassThru< CompressibleSolidBase >::execute( porousSolid, [=, &subRegion] ( auto & castedPorousSolid ) { typename TYPEOFREF( castedPorousSolid ) ::KernelWrapper porousWrapper = castedPorousSolid.createKernelUpdates(); - updatePorosityAndPermeabilityFromPressureAndAperture( porousWrapper, subRegion, pressure, oldHydraulicAperture, newHydraulicAperture ); + // if( m_isFixedStressPoromechanicsUpdate ) + // { + // arrayView1d< real64 const > const & pressure_n = subRegion.getField< flow::pressure_n >(); + // arrayView1d< real64 const > const & pressure_k = subRegion.getField< flow::pressure_k >(); + // arrayView1d< real64 const > const & temperature_n = subRegion.getField< flow::temperature_n >(); + // arrayView1d< real64 const > const & temperature_k = subRegion.getField< flow::temperature_k >(); + // updatePorosityAndPermeabilityFixedStress( porousWrapper, subRegion, pressure, pressure_k, pressure_n, temperature, temperature_k, temperature_n ); + // } + // else + if constexpr (std::is_same_v< typename TYPEOFREF( castedPorousSolid )::PermType, constitutive::ParallelPlatesPermeability >) { + updatePorosityAndPermeabilityFromPressureAndAperture( porousWrapper, subRegion, pressure, oldHydraulicAperture, newHydraulicAperture ); + } + else if constexpr ( std::is_same_v< typename TYPEOFREF( castedPorousSolid )::PermType, constitutive::SlipDependentPermeability > || + std::is_same_v< typename TYPEOFREF( castedPorousSolid )::PermType, constitutive::WillisRichardsPermeability > ) + { + updatePorosityAndPermeabilityFromPressurApertureJumpAndTraction(porousWrapper, + subRegion, pressure, oldHydraulicAperture, newHydraulicAperture, oldHydraulicAperture/*dHydraulicAperture_dNormalJump dummy*/, + dispJump, fractureTraction); + } } ); } From 8ea435f1727cab9b7a551daa5eaf465e86fc82a0 Mon Sep 17 00:00:00 2001 From: jacques franc Date: Tue, 17 Mar 2026 15:35:47 +0100 Subject: [PATCH 2/8] throw error if not one of the model --- .../fluidFlow/FlowSolverBase.cpp | 20 +++++++++---------- 1 file changed, 10 insertions(+), 10 deletions(-) diff --git a/src/coreComponents/physicsSolvers/fluidFlow/FlowSolverBase.cpp b/src/coreComponents/physicsSolvers/fluidFlow/FlowSolverBase.cpp index 6dea8e1506f..e738dfee774 100644 --- a/src/coreComponents/physicsSolvers/fluidFlow/FlowSolverBase.cpp +++ b/src/coreComponents/physicsSolvers/fluidFlow/FlowSolverBase.cpp @@ -676,15 +676,7 @@ void FlowSolverBase::updatePorosityAndPermeability( SurfaceElementSubRegion & su constitutive::ConstitutivePassThru< CompressibleSolidBase >::execute( porousSolid, [=, &subRegion] ( auto & castedPorousSolid ) { typename TYPEOFREF( castedPorousSolid ) ::KernelWrapper porousWrapper = castedPorousSolid.createKernelUpdates(); - // if( m_isFixedStressPoromechanicsUpdate ) - // { - // arrayView1d< real64 const > const & pressure_n = subRegion.getField< flow::pressure_n >(); - // arrayView1d< real64 const > const & pressure_k = subRegion.getField< flow::pressure_k >(); - // arrayView1d< real64 const > const & temperature_n = subRegion.getField< flow::temperature_n >(); - // arrayView1d< real64 const > const & temperature_k = subRegion.getField< flow::temperature_k >(); - // updatePorosityAndPermeabilityFixedStress( porousWrapper, subRegion, pressure, pressure_k, pressure_n, temperature, temperature_k, temperature_n ); - // } - // else + if constexpr (std::is_same_v< typename TYPEOFREF( castedPorousSolid )::PermType, constitutive::ParallelPlatesPermeability >) { updatePorosityAndPermeabilityFromPressureAndAperture( porousWrapper, subRegion, pressure, oldHydraulicAperture, newHydraulicAperture ); } @@ -692,9 +684,17 @@ void FlowSolverBase::updatePorosityAndPermeability( SurfaceElementSubRegion & su std::is_same_v< typename TYPEOFREF( castedPorousSolid )::PermType, constitutive::WillisRichardsPermeability > ) { updatePorosityAndPermeabilityFromPressurApertureJumpAndTraction(porousWrapper, - subRegion, pressure, oldHydraulicAperture, newHydraulicAperture, oldHydraulicAperture/*dHydraulicAperture_dNormalJump dummy*/, + subRegion, pressure, oldHydraulicAperture, newHydraulicAperture, oldHydraulicAperture/*dHydraulicAperture_dNormalJump dummy entry*/, dispJump, fractureTraction); } + else { + GEOS_ERROR( GEOS_FMT("{} permeability model is not yet supported. {}, {} and {} are supported.", + TYPEOFREF(castedPorousSolid)::PermType::catalogName(), + constitutive::ParallelPlatesPermeability::catalogName(), + constitutive::SlipDependentPermeability::catalogName(), + constitutive::WillisRichardsPermeability::catalogName() + ), getDataContext() ); + } } ); } From 9cd028f7f57d8047a9f62a24646efd26b45d1e01 Mon Sep 17 00:00:00 2001 From: Jian HUANG Date: Tue, 17 Mar 2026 15:30:29 -0500 Subject: [PATCH 3/8] update slip calculation for ALM --- .../contact/SolidMechanicsAugmentedLagrangianContact.cpp | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/src/coreComponents/physicsSolvers/solidMechanics/contact/SolidMechanicsAugmentedLagrangianContact.cpp b/src/coreComponents/physicsSolvers/solidMechanics/contact/SolidMechanicsAugmentedLagrangianContact.cpp index 17b720bbd64..c436004298c 100644 --- a/src/coreComponents/physicsSolvers/solidMechanics/contact/SolidMechanicsAugmentedLagrangianContact.cpp +++ b/src/coreComponents/physicsSolvers/solidMechanics/contact/SolidMechanicsAugmentedLagrangianContact.cpp @@ -762,9 +762,9 @@ void SolidMechanicsAugmentedLagrangianContact::implicitStepComplete( real64 cons GEOS_HOST_DEVICE ( localIndex const kfe ) { // Compute the slip - real64 const deltaDisp[2] = { deltaDispJump[kfe][1], - deltaDispJump[kfe][2] }; - slip[kfe] = LvArray::tensorOps::l2Norm< 2 >( deltaDisp ); + real64 const shearDisp[2] = { dispJump[kfe][1], + dispJump[kfe][2] }; + slip[kfe] = LvArray::tensorOps::l2Norm< 2 >( shearDisp ); // Compute current Tau and limit Tau real64 const tau[2] = { traction[kfe][1], From 00817750a789074742c9d30bbad24aca7b8b5bf4 Mon Sep 17 00:00:00 2001 From: Jian HUANG Date: Tue, 17 Mar 2026 15:40:03 -0500 Subject: [PATCH 4/8] reset bubble displacement after initialization --- .../solidMechanics/SolidMechanicsStateReset.cpp | 4 ++++ 1 file changed, 4 insertions(+) diff --git a/src/coreComponents/physicsSolvers/solidMechanics/SolidMechanicsStateReset.cpp b/src/coreComponents/physicsSolvers/solidMechanics/SolidMechanicsStateReset.cpp index 51086bdc70c..e15ff3a98f1 100644 --- a/src/coreComponents/physicsSolvers/solidMechanics/SolidMechanicsStateReset.cpp +++ b/src/coreComponents/physicsSolvers/solidMechanics/SolidMechanicsStateReset.cpp @@ -107,6 +107,10 @@ bool SolidMechanicsStateReset::execute( real64 const time_n, subRegion.getField< solidMechanics::averagePlasticStrain >().zero(); } ); + FaceManager & faceManager = mesh.getFaceManager(); + faceManager.getField< contact::totalBubbleDisplacement >().zero(); + faceManager.getField< contact::incrementalBubbleDisplacement >().zero(); + } // Option 2: enable / disable inelastic behavior From df3843c78be54d28be0d9f19e27fcf8c098140c5 Mon Sep 17 00:00:00 2001 From: Jian HUANG Date: Tue, 17 Mar 2026 15:48:40 -0500 Subject: [PATCH 5/8] add ExponentialDecayPermeability model --- .../physicsSolvers/fluidFlow/FlowSolverBase.cpp | 9 ++++++--- 1 file changed, 6 insertions(+), 3 deletions(-) diff --git a/src/coreComponents/physicsSolvers/fluidFlow/FlowSolverBase.cpp b/src/coreComponents/physicsSolvers/fluidFlow/FlowSolverBase.cpp index e738dfee774..5967721040f 100644 --- a/src/coreComponents/physicsSolvers/fluidFlow/FlowSolverBase.cpp +++ b/src/coreComponents/physicsSolvers/fluidFlow/FlowSolverBase.cpp @@ -35,9 +35,10 @@ #include "physicsSolvers/fluidFlow/FlowSolverBaseFields.hpp" #include "physicsSolvers/solidMechanics/contact/ContactFields.hpp" //should not be here -- testing only design flaw +#include "constitutive/permeability/ExponentialDecayPermeability.hpp" +#include "constitutive/permeability/ParallelPlatesPermeability.hpp" #include "constitutive/permeability/SlipDependentPermeability.hpp" #include "constitutive/permeability/WillisRichardsPermeability.hpp" -#include "constitutive/permeability/ParallelPlatesPermeability.hpp" #include "physicsSolvers/fluidFlow/kernels/MinPoreVolumeMaxPorosityKernel.hpp" #include "physicsSolvers/fluidFlow/kernels/StencilWeightsUpdateKernel.hpp" @@ -681,7 +682,8 @@ void FlowSolverBase::updatePorosityAndPermeability( SurfaceElementSubRegion & su updatePorosityAndPermeabilityFromPressureAndAperture( porousWrapper, subRegion, pressure, oldHydraulicAperture, newHydraulicAperture ); } else if constexpr ( std::is_same_v< typename TYPEOFREF( castedPorousSolid )::PermType, constitutive::SlipDependentPermeability > || - std::is_same_v< typename TYPEOFREF( castedPorousSolid )::PermType, constitutive::WillisRichardsPermeability > ) + std::is_same_v< typename TYPEOFREF( castedPorousSolid )::PermType, constitutive::WillisRichardsPermeability > || + std::is_same_v< typename TYPEOFREF( castedPorousSolid )::PermType, constitutive::ExponentialDecayPermeability > ) { updatePorosityAndPermeabilityFromPressurApertureJumpAndTraction(porousWrapper, subRegion, pressure, oldHydraulicAperture, newHydraulicAperture, oldHydraulicAperture/*dHydraulicAperture_dNormalJump dummy entry*/, @@ -692,7 +694,8 @@ void FlowSolverBase::updatePorosityAndPermeability( SurfaceElementSubRegion & su TYPEOFREF(castedPorousSolid)::PermType::catalogName(), constitutive::ParallelPlatesPermeability::catalogName(), constitutive::SlipDependentPermeability::catalogName(), - constitutive::WillisRichardsPermeability::catalogName() + constitutive::WillisRichardsPermeability::catalogName(), + constitutive::ExponentialDecayPermeability::catalogName() ), getDataContext() ); } From 581dc462d09d7e48760691ac650cee9a7c3e197e Mon Sep 17 00:00:00 2001 From: Jian HUANG Date: Tue, 17 Mar 2026 15:59:54 -0500 Subject: [PATCH 6/8] add a check for totalBubbleDisplacement --- .../solidMechanics/SolidMechanicsStateReset.cpp | 7 +++++-- 1 file changed, 5 insertions(+), 2 deletions(-) diff --git a/src/coreComponents/physicsSolvers/solidMechanics/SolidMechanicsStateReset.cpp b/src/coreComponents/physicsSolvers/solidMechanics/SolidMechanicsStateReset.cpp index e15ff3a98f1..2db864d5f02 100644 --- a/src/coreComponents/physicsSolvers/solidMechanics/SolidMechanicsStateReset.cpp +++ b/src/coreComponents/physicsSolvers/solidMechanics/SolidMechanicsStateReset.cpp @@ -108,8 +108,11 @@ bool SolidMechanicsStateReset::execute( real64 const time_n, } ); FaceManager & faceManager = mesh.getFaceManager(); - faceManager.getField< contact::totalBubbleDisplacement >().zero(); - faceManager.getField< contact::incrementalBubbleDisplacement >().zero(); + if( faceManager.hasField< contact::totalBubbleDisplacement >() ) + { + faceManager.getField< contact::totalBubbleDisplacement >().zero(); + faceManager.getField< contact::incrementalBubbleDisplacement >().zero(); + } } From 3cbd051557998d1016e80105f54867243cacceef Mon Sep 17 00:00:00 2001 From: jacques franc Date: Wed, 18 Mar 2026 11:56:21 +0100 Subject: [PATCH 7/8] missing include --- .../physicsSolvers/solidMechanics/SolidMechanicsStateReset.cpp | 1 + 1 file changed, 1 insertion(+) diff --git a/src/coreComponents/physicsSolvers/solidMechanics/SolidMechanicsStateReset.cpp b/src/coreComponents/physicsSolvers/solidMechanics/SolidMechanicsStateReset.cpp index 2db864d5f02..73e082f34ce 100644 --- a/src/coreComponents/physicsSolvers/solidMechanics/SolidMechanicsStateReset.cpp +++ b/src/coreComponents/physicsSolvers/solidMechanics/SolidMechanicsStateReset.cpp @@ -20,6 +20,7 @@ #include "SolidMechanicsStateReset.hpp" #include "physicsSolvers/PhysicsSolverManager.hpp" +#include "physicsSolvers/solidMechanics/contact/ContactFields.hpp" //should not be here -- testing only design flaw #include "physicsSolvers/solidMechanics/SolidMechanicsLagrangianFEM.hpp" #include "physicsSolvers/LogLevelsInfo.hpp" #include "mesh/DomainPartition.hpp" From f0d2b21843e8b3e4a25a751761dbfb3b6ac61468 Mon Sep 17 00:00:00 2001 From: jacques franc Date: Wed, 18 Mar 2026 18:53:47 +0100 Subject: [PATCH 8/8] more reset --- .../solidMechanics/SolidMechanicsStateReset.cpp | 10 ++++++++++ 1 file changed, 10 insertions(+) diff --git a/src/coreComponents/physicsSolvers/solidMechanics/SolidMechanicsStateReset.cpp b/src/coreComponents/physicsSolvers/solidMechanics/SolidMechanicsStateReset.cpp index 73e082f34ce..b3d68b7def7 100644 --- a/src/coreComponents/physicsSolvers/solidMechanics/SolidMechanicsStateReset.cpp +++ b/src/coreComponents/physicsSolvers/solidMechanics/SolidMechanicsStateReset.cpp @@ -108,6 +108,16 @@ bool SolidMechanicsStateReset::execute( real64 const time_n, subRegion.getField< solidMechanics::averagePlasticStrain >().zero(); } ); + + elemManager.forElementSubRegions< SurfaceElementSubRegion >( regionNames, + [&]( localIndex const, + ElementSubRegionBase & subRegion ) + { + subRegion.getField< contact::dispJump >().zero(); //MAYBE NOT IF Fault Failing from zero ! + subRegion.getField< contact::slip >().zero(); + } ); + + FaceManager & faceManager = mesh.getFaceManager(); if( faceManager.hasField< contact::totalBubbleDisplacement >() ) {