feat: Rate-and-state friction with explicit time integration (#3450)

* Add rate-and-state kernel for explicit time stepping and quasi-dynamic solver using Runge-Kutta 3(2) * Update strings for added rate-and-state fields * Change to using a single multi-d array for the Runge-Kutta stage rates * Add PID controller to adaptive time step update * Bracket slip rate after Newton update * Add kernel for EmbeddedRungeKutta time stepping of slip and state and refactor QuasiDynamicEQRK32 solver * Add some comments to QuasiDynamicEQRK32 solver * Add Bogacki-Shampine 3(2) Runge-Kutta method * Add PIDController class for time step control * Add xml input for explicit solver, update subrepo, add docs * Rename constructor arguments to prevent shadowing * rebaseline --------- Co-authored-by: Matteo Cusini <[email protected]>
GEOS-DEV · Dec 10, 2024 · 7e2c33b · 7e2c33b
1 parent eef8de4
commit 7e2c33b
Show file tree

Hide file tree

Showing 16 changed files with 1,539 additions and 108 deletions.
diff --git a/.integrated_tests.yaml b/.integrated_tests.yaml
@@ -1,6 +1,6 @@
 baselines:
   bucket: geosx
-  baseline: integratedTests/baseline_integratedTests-pr3480-9217-caaecb8
+  baseline: integratedTests/baseline_integratedTests-pr3450-9221-37d940c
 allow_fail:
   all: ''
   streak: ''
diff --git a/BASELINE_NOTES.md b/BASELINE_NOTES.md
@@ -6,6 +6,10 @@ This file is designed to track changes to the integrated test baselines.
 Any developer who updates the baseline ID in the .integrated_tests.yaml file is expected to create an entry in this file with the pull request number, date, and their justification for rebaselining.
 These notes should be in reverse-chronological order, and use the following time format: (YYYY-MM-DD).
 
+PR #3450 (2024-12-08)
+=====================
+Added test for explicit runge kutta sprinslider.
+
 PR #3480 (2024-12-06)
 =====================
 Add "logLevel" parameter under /Problem/Outputs in baseline files

diff --git a/inputFiles/inducedSeismicity/SpringSliderExplicit_base.xml b/inputFiles/inducedSeismicity/SpringSliderExplicit_base.xml
@@ -0,0 +1,167 @@
+<?xml version="1.0" ?>
+<Problem>
+  <Solvers>
+    <QuasiDynamicEQRK32
+      name="SpringSlider"
+      targetRegions="{ Fault }"
+      shearImpedance="4.41"
+      initialDt="1e-5"
+      logLevel="1"
+      discretization="FE1">
+    </QuasiDynamicEQRK32>
+
+    <SurfaceGenerator
+      name="SurfaceGen"
+      targetRegions="{ Domain }"
+      rockToughness="1.0"
+      mpiCommOrder="1"
+      fractureRegion="Fault"/>  
+  </Solvers>
+
+  <NumericalMethods>
+    <FiniteElements>
+      <FiniteElementSpace
+        name="FE1"
+        order="1"/>
+    </FiniteElements>
+  </NumericalMethods>
+
+  <Mesh>
+    <InternalMesh
+      name="mesh"
+      elementTypes="{ C3D8 }"
+      xCoords="{ 0, 1 }"
+      yCoords="{ 0, 2 }"
+      zCoords="{ 0, 1 }"
+      nx="{ 1 }"
+      ny="{ 2 }"
+      nz="{ 1 }"
+      cellBlockNames="{ cb1 }"/>
+  </Mesh>
+
+  <Geometry>    
+    <ThickPlane
+      name="faultPlane"
+      normal="{ 0, 1, 0 }"
+      origin="{ 0, 1, 0 }"
+      thickness="0.1"/>
+  </Geometry>
+
+  <ElementRegions>
+    <CellElementRegion
+      name="Domain"
+      cellBlocks="{ cb1 }"
+      materialList="{}"/>
+
+    <SurfaceElementRegion
+      name="Fault"
+      materialList="{frictionLaw}"
+      defaultAperture="1e-3"/>  
+    </ElementRegions>
+
+  <Constitutive>
+    <RateAndStateFriction
+      name="frictionLaw"
+      defaultA="0.01"
+      defaultB="0.015"
+      defaultDc="1.0e-5"  
+      defaultReferenceVelocity="1.0e-6"
+      defaultReferenceFrictionCoefficient="0.6"/>  
+  </Constitutive>
+
+  <FieldSpecifications>
+    <FieldSpecification
+      name="fault"
+      fieldName="ruptureState"
+      initialCondition="1"
+      objectPath="faceManager"
+      scale="1"
+      setNames="{ faultPlane }"/>
+
+    <FieldSpecification
+      name="normalTraction"
+      fieldName="traction"
+      initialCondition="1"
+      objectPath="ElementRegions/Fault/FractureSubRegion"
+      component="0"
+      scale="50"
+      setNames="{all}"/>
+
+    <FieldSpecification
+      name="shearTraction1"
+      fieldName="traction"
+      initialCondition="1"
+      objectPath="ElementRegions/Fault/FractureSubRegion"
+      component="1"
+      scale="21.2132034356"
+      setNames="{all}"/>
+    <FieldSpecification
+      name="shearTraction2"
+      fieldName="traction"
+      initialCondition="1"
+      objectPath="ElementRegions/Fault/FractureSubRegion"
+      component="2"
+      scale="21.2132034356"
+      setNames="{all}"/>  
+    <FieldSpecification
+      name="stateVariable"
+      fieldName="stateVariable"
+      initialCondition="1"
+      objectPath="ElementRegions/Fault/FractureSubRegion"
+      scale="0.6"
+      setNames="{all}"/>
+    <FieldSpecification
+      name="slipRate"
+      fieldName="slipRate"
+      initialCondition="1"
+      objectPath="ElementRegions/Fault/FractureSubRegion"
+      scale="1.0e-6"
+      setNames="{all}"/>
+    <FieldSpecification
+      name="slipVelocity1"
+      fieldName="slipVelocity"
+      initialCondition="1"
+      objectPath="ElementRegions/Fault/FractureSubRegion"
+      component="0"
+      scale="0.70710678118e-6"
+      setNames="{all}"/>               
+    <FieldSpecification
+      name="slipVelocity2"
+      fieldName="slipVelocity"
+      initialCondition="1"
+      objectPath="ElementRegions/Fault/FractureSubRegion"
+      component="1"
+      scale="0.70710678118e-6"
+      setNames="{all}"/>               
+  </FieldSpecifications>
+
+  <Outputs>
+    <VTK
+      name="vtkOutput"
+      plotFileRoot="springSliderExplicit"/>
+    <Restart
+      name="restart"/>  
+
+    <TimeHistory
+      name="timeHistoryOutput"
+      sources="{/Tasks/slipCollection,/Tasks/slipRateCollection,/Tasks/stateVariableCollection}"
+      filename="springSliderExplicit"/>
+  </Outputs>
+
+  <Tasks>
+    <PackCollection
+      name="slipCollection"
+      objectPath="ElementRegions/Fault/FractureSubRegion"
+      fieldName="displacementJump"/>
+
+    <PackCollection
+      name="slipRateCollection"
+      objectPath="ElementRegions/Fault/FractureSubRegion"
+      fieldName="slipRate"/>
+
+    <PackCollection
+      name="stateVariableCollection"
+      objectPath="ElementRegions/Fault/FractureSubRegion"
+      fieldName="stateVariable"/>  
+  </Tasks>
+</Problem>
diff --git a/inputFiles/inducedSeismicity/SpringSliderExplicit_smoke.xml b/inputFiles/inducedSeismicity/SpringSliderExplicit_smoke.xml
@@ -0,0 +1,32 @@
+<?xml version="1.0" ?>
+<Problem>
+  <Included>
+    <File
+      name="./SpringSliderExplicit_base.xml"/>
+  </Included>
+
+  <Events
+    maxTime="4.5e4">
+
+    <SoloEvent 
+      name="generateFault"
+      target="/Solvers/SurfaceGen"/>
+
+    <PeriodicEvent
+      name="vtkOutput"
+      cycleFrequency="1"
+      targetExactTimestep="0"
+      target="/Outputs/vtkOutput"/>  
+
+    <PeriodicEvent
+      name="solverApplications"
+      maxEventDt="1e4"
+      target="/Solvers/SpringSlider"/>
+
+    <PeriodicEvent
+      name="resarts"
+      timeFrequency="2e4"
+      targetExactTimestep="0"
+      target="/Outputs/restart"/>  
+  </Events>
+</Problem>
diff --git a/inputFiles/inducedSeismicity/SpringSlider_base.xml b/inputFiles/inducedSeismicity/SpringSlider_base.xml
@@ -123,15 +123,15 @@
       initialCondition="1"
       objectPath="ElementRegions/Fault/FractureSubRegion"
       component="0"
-      scale="0."
+      scale="0.70710678118e-6"
       setNames="{all}"/>               
     <FieldSpecification
       name="slipVelocity2"
       fieldName="slipVelocity"
       initialCondition="1"
       objectPath="ElementRegions/Fault/FractureSubRegion"
       component="1"
-      scale="0."
+      scale="0.70710678118e-6"
       setNames="{all}"/>               
   </FieldSpecifications>
 

diff --git a/inputFiles/inducedSeismicity/inducedSeismicity.ats b/inputFiles/inducedSeismicity/inducedSeismicity.ats
@@ -28,6 +28,13 @@ decks = [
         partitions=((1, 1, 1), ),
         restart_step=0,
         check_step=3262,
-        restartcheck_params=RestartcheckParameters(atol=1e-4, rtol=1e-3)) 
+        restartcheck_params=RestartcheckParameters(atol=1e-4, rtol=1e-3)),
+    TestDeck(
+        name="SpringSliderExplicit_smoke",
+        description="Spring slider 0D system",
+        partitions=((1, 1, 1), ),
+        restart_step=0,
+        check_step=532,
+        restartcheck_params=RestartcheckParameters(atol=1e-4, rtol=1e-3))     
     ]
 generate_geos_tests(decks)
diff --git a/src/coreComponents/physicsSolvers/contact/ContactFields.hpp b/src/coreComponents/physicsSolvers/contact/ContactFields.hpp
@@ -74,6 +74,14 @@ DECLARE_FIELD( dispJump,
                WRITE_AND_READ,
                "Displacement jump vector in the local reference system" );
 
+DECLARE_FIELD( dispJump_n,
+               "displacementJump",
+               array2d< real64 >,
+               0,
+               NOPLOT,
+               WRITE_AND_READ,
+               "Displacement jump vector in the local reference system at the current time-step" );
+
 DECLARE_FIELD( slip,
                "slip",
                array1d< real64 >,
@@ -114,6 +122,14 @@ DECLARE_FIELD( traction,
                WRITE_AND_READ,
                "Fracture traction vector in the local reference system." );
 
+DECLARE_FIELD( traction_n,
+               "traction_n",
+               array2d< real64 >,
+               0,
+               NOPLOT,
+               WRITE_AND_READ,
+               "Initial fracture traction vector in the local reference system at this time-step." );
+
 DECLARE_FIELD( deltaTraction,
                "deltaTraction",
                array2d< real64 >,

diff --git a/src/coreComponents/physicsSolvers/inducedSeismicity/CMakeLists.txt b/src/coreComponents/physicsSolvers/inducedSeismicity/CMakeLists.txt
@@ -4,6 +4,7 @@ set( physicsSolvers_headers
      inducedSeismicity/inducedSeismicityFields.hpp
      inducedSeismicity/rateAndStateFields.hpp
      inducedSeismicity/QuasiDynamicEQ.hpp
+     inducedSeismicity/QuasiDynamicEQRK32.hpp
      inducedSeismicity/SeismicityRate.hpp
      inducedSeismicity/kernels/RateAndStateKernels.hpp
      inducedSeismicity/kernels/SeismicityRateKernels.hpp 
@@ -12,6 +13,7 @@ set( physicsSolvers_headers
 # Specify solver sources
 set( physicsSolvers_sources
      ${physicsSolvers_sources}
-     inducedSeismicity/QuasiDynamicEQ.cpp 
+     inducedSeismicity/QuasiDynamicEQ.cpp
+     inducedSeismicity/QuasiDynamicEQRK32.cpp 
      inducedSeismicity/SeismicityRate.cpp 
      PARENT_SCOPE )
diff --git a/src/coreComponents/physicsSolvers/inducedSeismicity/QuasiDynamicEQ.cpp b/src/coreComponents/physicsSolvers/inducedSeismicity/QuasiDynamicEQ.cpp
@@ -32,6 +32,7 @@ namespace geos
 using namespace dataRepository;
 using namespace fields;
 using namespace constitutive;
+using namespace rateAndStateKernels;
 
 QuasiDynamicEQ::QuasiDynamicEQ( const string & name,
                                 Group * const parent ):
@@ -149,8 +150,8 @@ real64 QuasiDynamicEQ::solverStep( real64 const & time_n,
 
   /// 2. Solve for slip rate and state variable and, compute slip
   GEOS_LOG_LEVEL_RANK_0( 1, "Rate and State solver" );
-
-  integer const maxNewtonIter = m_nonlinearSolverParameters.m_maxIterNewton;
+  integer const maxIterNewton = m_nonlinearSolverParameters.m_maxIterNewton;
+  real64 const newtonTol = m_nonlinearSolverParameters.m_newtonTol;
   forDiscretizationOnMeshTargets( domain.getMeshBodies(), [&]( string const &,
                                                                MeshLevel & mesh,
                                                                arrayView1d< string const > const & regionNames )
@@ -161,7 +162,8 @@ real64 QuasiDynamicEQ::solverStep( real64 const & time_n,
                                                                                 SurfaceElementSubRegion & subRegion )
     {
       // solve rate and state equations.
-      rateAndStateKernels::createAndLaunch< parallelDevicePolicy<> >( subRegion, viewKeyStruct::frictionLawNameString(), m_shearImpedance, maxNewtonIter, time_n, dtStress );
+      createAndLaunch< ImplicitFixedStressRateAndStateKernel, parallelDevicePolicy<> >( subRegion, viewKeyStruct::frictionLawNameString(), m_shearImpedance, maxIterNewton, newtonTol, time_n,
+                                                                                        dtStress );
       // save old state
       saveOldStateAndUpdateSlip( subRegion, dtStress );
     } );

diff --git a/src/coreComponents/physicsSolvers/inducedSeismicity/QuasiDynamicEQ.hpp b/src/coreComponents/physicsSolvers/inducedSeismicity/QuasiDynamicEQ.hpp
@@ -47,7 +47,7 @@ class QuasiDynamicEQ : public PhysicsSolverBase
   struct viewKeyStruct : public PhysicsSolverBase::viewKeyStruct
   {
     /// stress solver name
-    static constexpr char const * stressSolverNameString() { return "stressSolverName"; }
+    constexpr static char const * stressSolverNameString() { return "stressSolverName"; }
     /// Friction law name string
     constexpr static char const * frictionLawNameString() { return "frictionLawName"; }
     /// Friction law name string