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..5967721040f 100644 --- a/src/coreComponents/physicsSolvers/fluidFlow/FlowSolverBase.cpp +++ b/src/coreComponents/physicsSolvers/fluidFlow/FlowSolverBase.cpp @@ -33,6 +33,13 @@ #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/ExponentialDecayPermeability.hpp" +#include "constitutive/permeability/ParallelPlatesPermeability.hpp" +#include "constitutive/permeability/SlipDependentPermeability.hpp" +#include "constitutive/permeability/WillisRichardsPermeability.hpp" + #include "physicsSolvers/fluidFlow/kernels/MinPoreVolumeMaxPorosityKernel.hpp" #include "physicsSolvers/fluidFlow/kernels/StencilWeightsUpdateKernel.hpp" @@ -104,6 +111,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 +668,36 @@ 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 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 > || + std::is_same_v< typename TYPEOFREF( castedPorousSolid )::PermType, constitutive::ExponentialDecayPermeability > ) + { + updatePorosityAndPermeabilityFromPressurApertureJumpAndTraction(porousWrapper, + 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(), + constitutive::ExponentialDecayPermeability::catalogName() + ), getDataContext() ); + } } ); } diff --git a/src/coreComponents/physicsSolvers/solidMechanics/SolidMechanicsStateReset.cpp b/src/coreComponents/physicsSolvers/solidMechanics/SolidMechanicsStateReset.cpp index 51086bdc70c..b3d68b7def7 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" @@ -107,6 +108,23 @@ 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 >() ) + { + faceManager.getField< contact::totalBubbleDisplacement >().zero(); + faceManager.getField< contact::incrementalBubbleDisplacement >().zero(); + } + } // Option 2: enable / disable inelastic behavior 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],