fix: Fault Perm Update for Contact solvers

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 ); }
 

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()  );
+      }
 
   } );
 }

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

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],