diff --git a/.integrated_tests.yaml b/.integrated_tests.yaml
index 3aa31a7e4c..31ce96119a 100644
--- 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
index acbeb39e38..3b09722325 100644
--- 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
new file mode 100644
index 0000000000..a0c0381fb2
--- /dev/null
+++ b/inputFiles/inducedSeismicity/SpringSliderExplicit_base.xml
@@ -0,0 +1,167 @@
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
\ No newline at end of file
diff --git a/inputFiles/inducedSeismicity/SpringSliderExplicit_smoke.xml b/inputFiles/inducedSeismicity/SpringSliderExplicit_smoke.xml
new file mode 100644
index 0000000000..2de09b79a0
--- /dev/null
+++ b/inputFiles/inducedSeismicity/SpringSliderExplicit_smoke.xml
@@ -0,0 +1,32 @@
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
\ No newline at end of file
diff --git a/inputFiles/inducedSeismicity/SpringSlider_base.xml b/inputFiles/inducedSeismicity/SpringSlider_base.xml
index 9f14103160..5f6aefc94a 100644
--- a/inputFiles/inducedSeismicity/SpringSlider_base.xml
+++ b/inputFiles/inducedSeismicity/SpringSlider_base.xml
@@ -123,7 +123,7 @@
initialCondition="1"
objectPath="ElementRegions/Fault/FractureSubRegion"
component="0"
- scale="0."
+ scale="0.70710678118e-6"
setNames="{all}"/>
diff --git a/inputFiles/inducedSeismicity/inducedSeismicity.ats b/inputFiles/inducedSeismicity/inducedSeismicity.ats
index 08b031177f..813bdf08b8 100644
--- 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
index 9b712c0e70..df37934c0c 100644
--- 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
index 7edf040454..c8f23c5532 100644
--- 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
index 50ca6bb3cc..3ddf5fc149 100644
--- 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
index e46d608d8c..6d85175d3c 100644
--- 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
diff --git a/src/coreComponents/physicsSolvers/inducedSeismicity/QuasiDynamicEQRK32.cpp b/src/coreComponents/physicsSolvers/inducedSeismicity/QuasiDynamicEQRK32.cpp
new file mode 100644
index 0000000000..7b56911893
--- /dev/null
+++ b/src/coreComponents/physicsSolvers/inducedSeismicity/QuasiDynamicEQRK32.cpp
@@ -0,0 +1,478 @@
+/*
+ * ------------------------------------------------------------------------------------------------------------
+ * SPDX-License-Identifier: LGPL-2.1-only
+ *
+ * Copyright (c) 2016-2024 Lawrence Livermore National Security LLC
+ * Copyright (c) 2018-2024 Total, S.A
+ * 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 QuasiDynamicEQRK32.cpp
+ */
+
+#include "QuasiDynamicEQRK32.hpp"
+
+#include "dataRepository/InputFlags.hpp"
+#include "mesh/DomainPartition.hpp"
+#include "kernels/RateAndStateKernels.hpp"
+#include "rateAndStateFields.hpp"
+#include "physicsSolvers/contact/ContactFields.hpp"
+#include "fieldSpecification/FieldSpecificationManager.hpp"
+
+namespace geos
+{
+
+using namespace dataRepository;
+using namespace fields;
+using namespace constitutive;
+using namespace rateAndStateKernels;
+
+QuasiDynamicEQRK32::QuasiDynamicEQRK32( const string & name,
+ Group * const parent ):
+ PhysicsSolverBase( name, parent ),
+ m_stressSolver( nullptr ),
+ m_stressSolverName( "SpringSlider" ),
+ m_shearImpedance( 0.0 ),
+ m_butcherTable( BogackiShampine32Table()), // TODO: The butcher table should be specified in the XML input.
+ m_successfulStep( false ),
+ m_controller( PIDController( { 1.0/18.0, 1.0/9.0, 1.0/18.0 },
+ 1.0e-6, 1.0e-6, 0.81 )) // TODO: The control parameters should be specified in the XML input
+{
+ this->registerWrapper( viewKeyStruct::shearImpedanceString(), &m_shearImpedance ).
+ setInputFlag( InputFlags::REQUIRED ).
+ setDescription( "Shear impedance." );
+
+ this->registerWrapper( viewKeyStruct::stressSolverNameString(), &m_stressSolverName ).
+ setInputFlag( InputFlags::OPTIONAL ).
+ setDescription( "Name of solver for computing stress. If empty, the spring-slider model is run." );
+}
+
+void QuasiDynamicEQRK32::postInputInitialization()
+{
+
+ // Initialize member stress solver as specified in XML input
+ if( !m_stressSolverName.empty() )
+ {
+ m_stressSolver = &this->getParent().getGroup< PhysicsSolverBase >( m_stressSolverName );
+ }
+
+ PhysicsSolverBase::postInputInitialization();
+}
+
+QuasiDynamicEQRK32::~QuasiDynamicEQRK32()
+{
+ // TODO Auto-generated destructor stub
+}
+
+
+void QuasiDynamicEQRK32::registerDataOnMesh( Group & meshBodies )
+{
+ PhysicsSolverBase::registerDataOnMesh( meshBodies );
+
+ forDiscretizationOnMeshTargets( meshBodies, [&] ( string const &,
+ MeshLevel & mesh,
+ arrayView1d< string const > const & regionNames )
+ {
+ ElementRegionManager & elemManager = mesh.getElemManager();
+
+ elemManager.forElementSubRegions< SurfaceElementSubRegion >( regionNames,
+ [&]( localIndex const,
+ SurfaceElementSubRegion & subRegion )
+ {
+ // Scalar functions on fault
+ subRegion.registerField< rateAndState::stateVariable >( getName() );
+ subRegion.registerField< rateAndState::stateVariable_n >( getName() );
+ subRegion.registerField< rateAndState::slipRate >( getName() );
+
+ // Tangent (2-component) functions on fault
+ string const labels2Comp[2] = {"tangent1", "tangent2" };
+ subRegion.registerField< rateAndState::slipVelocity >( getName() ).
+ setDimLabels( 1, labels2Comp ).reference().resizeDimension< 1 >( 2 );
+ subRegion.registerField< rateAndState::slipVelocity_n >( getName() ).
+ setDimLabels( 1, labels2Comp ).reference().resizeDimension< 1 >( 2 );
+ subRegion.registerField< rateAndState::deltaSlip >( getName() ).
+ setDimLabels( 1, labels2Comp ).reference().resizeDimension< 1 >( 2 );
+ subRegion.registerField< rateAndState::deltaSlip_n >( getName() ).
+ setDimLabels( 1, labels2Comp ).reference().resizeDimension< 1 >( 2 );
+
+ // Runge-Kutta stage rates and error
+ integer const numRKComponents = 3;
+ subRegion.registerField< rateAndState::rungeKuttaStageRates >( getName() ).reference().resizeDimension< 1, 2 >( m_butcherTable.numStages, numRKComponents );
+ subRegion.registerField< rateAndState::error >( getName() ).reference().resizeDimension< 1 >( numRKComponents );
+
+
+ if( !subRegion.hasWrapper( contact::dispJump::key() ))
+ {
+ // 3-component functions on fault
+ string const labels3Comp[3] = { "normal", "tangent1", "tangent2" };
+ subRegion.registerField< contact::dispJump >( getName() ).
+ setDimLabels( 1, labels3Comp ).
+ reference().resizeDimension< 1 >( 3 );
+ subRegion.registerField< contact::dispJump_n >( getName() ).
+ setDimLabels( 1, labels3Comp ).
+ reference().resizeDimension< 1 >( 3 );
+ subRegion.registerField< contact::traction >( getName() ).
+ setDimLabels( 1, labels3Comp ).
+ reference().resizeDimension< 1 >( 3 );
+ subRegion.registerField< contact::traction_n >( getName() ).
+ setDimLabels( 1, labels3Comp ).
+ reference().resizeDimension< 1 >( 3 );
+
+ subRegion.registerWrapper< string >( viewKeyStruct::frictionLawNameString() ).
+ setPlotLevel( PlotLevel::NOPLOT ).
+ setRestartFlags( RestartFlags::NO_WRITE ).
+ setSizedFromParent( 0 );
+
+ string & frictionLawName = subRegion.getReference< string >( viewKeyStruct::frictionLawNameString() );
+ frictionLawName = PhysicsSolverBase::getConstitutiveName< FrictionBase >( subRegion );
+ GEOS_ERROR_IF( frictionLawName.empty(), GEOS_FMT( "{}: FrictionBase model not found on subregion {}",
+ getDataContext(), subRegion.getDataContext() ) );
+ }
+ } );
+ } );
+}
+
+real64 QuasiDynamicEQRK32::solverStep( real64 const & time_n,
+ real64 const & dt,
+ int const cycleNumber,
+ DomainPartition & domain )
+{
+ if( cycleNumber == 0 )
+ {
+ /// Apply initial conditions to the Fault
+ FieldSpecificationManager & fieldSpecificationManager = FieldSpecificationManager::getInstance();
+
+ forDiscretizationOnMeshTargets( domain.getMeshBodies(), [&]( string const &,
+ MeshLevel & mesh,
+ arrayView1d< string const > const & )
+
+ {
+ fieldSpecificationManager.applyInitialConditions( mesh );
+
+ } );
+ saveState( domain );
+ }
+
+ real64 dtAdaptive = dt;
+
+ GEOS_LOG_LEVEL_RANK_0( 1, "Begin adaptive time step" );
+ while( true ) // Adaptive time step loop. Performs a Runge-Kutta time stepping with error control on state and slip
+ {
+ real64 dtStress; GEOS_UNUSED_VAR( dtStress );
+
+ // Initial Runge-Kutta stage
+ stepRateStateODEInitialSubstage( dtAdaptive, domain );
+ real64 dtStage = m_butcherTable.c[1]*dtAdaptive;
+ dtStress = updateStresses( time_n, dtStage, cycleNumber, domain );
+ updateSlipVelocity( time_n, dtStage, domain );
+
+ // Remaining stages
+ for( integer stageIndex = 1; stageIndex < m_butcherTable.numStages-1; stageIndex++ )
+ {
+ stepRateStateODESubstage( stageIndex, dtAdaptive, domain );
+ dtStage = m_butcherTable.c[stageIndex+1]*dtAdaptive;
+ dtStress = updateStresses( time_n, dtStage, cycleNumber, domain );
+ updateSlipVelocity( time_n, dtStage, domain );
+ }
+
+ stepRateStateODEAndComputeError( dtAdaptive, domain );
+ // Update timestep based on the time step error
+ real64 const dtNext = setNextDt( dtAdaptive, domain );
+ if( m_successfulStep ) // set in setNextDt
+ {
+ // Compute stresses, and slip velocity and save results at updated time,
+ if( !m_butcherTable.FSAL )
+ {
+ dtStress = updateStresses( time_n, dtAdaptive, cycleNumber, domain );
+ updateSlipVelocity( time_n, dtAdaptive, domain );
+ }
+ saveState( domain );
+ // update the time step and exit the adaptive time step loop
+ dtAdaptive = dtNext;
+ break;
+ }
+ else
+ {
+ // Retry with updated time step
+ dtAdaptive = dtNext;
+ }
+ }
+ // return time step size achieved by stress solver
+ return dtAdaptive;
+}
+
+void QuasiDynamicEQRK32::stepRateStateODEInitialSubstage( real64 const dt, DomainPartition & domain ) const
+{
+
+ forDiscretizationOnMeshTargets( domain.getMeshBodies(), [&]( string const &,
+ MeshLevel & mesh,
+ arrayView1d< string const > const & regionNames )
+
+ {
+ mesh.getElemManager().forElementSubRegions< SurfaceElementSubRegion >( regionNames,
+ [&]( localIndex const,
+ SurfaceElementSubRegion & subRegion )
+ {
+
+ string const & fricitonLawName = subRegion.template getReference< string >( viewKeyStruct::frictionLawNameString() );
+ RateAndStateFriction const & frictionLaw = getConstitutiveModel< RateAndStateFriction >( subRegion, fricitonLawName );
+ rateAndStateKernels::EmbeddedRungeKuttaKernel rkKernel( subRegion, frictionLaw, m_butcherTable );
+ arrayView3d< real64 > const rkStageRates = subRegion.getField< rateAndState::rungeKuttaStageRates >();
+
+ if( m_butcherTable.FSAL && m_successfulStep )
+ {
+ forAll< parallelDevicePolicy<> >( subRegion.size(), [=] GEOS_HOST_DEVICE ( localIndex const k )
+ {
+ rkKernel.updateStageRatesFSAL( k );
+ rkKernel.updateStageValues( k, 1, dt );
+ } );
+ }
+ else
+ {
+ forAll< parallelDevicePolicy<> >( subRegion.size(), [=] GEOS_HOST_DEVICE ( localIndex const k )
+ {
+ rkKernel.initialize( k );
+ rkKernel.updateStageRates( k, 0 );
+ rkKernel.updateStageValues( k, 1, dt );
+ } );
+ }
+ } );
+ } );
+}
+
+void QuasiDynamicEQRK32::stepRateStateODESubstage( integer const stageIndex,
+ real64 const dt,
+ DomainPartition & domain ) const
+{
+
+ forDiscretizationOnMeshTargets( domain.getMeshBodies(), [&]( string const &,
+ MeshLevel & mesh,
+ arrayView1d< string const > const & regionNames )
+
+ {
+ mesh.getElemManager().forElementSubRegions< SurfaceElementSubRegion >( regionNames,
+ [&]( localIndex const,
+ SurfaceElementSubRegion & subRegion )
+ {
+
+ string const & fricitonLawName = subRegion.template getReference< string >( viewKeyStruct::frictionLawNameString() );
+ RateAndStateFriction const & frictionLaw = getConstitutiveModel< RateAndStateFriction >( subRegion, fricitonLawName );
+ rateAndStateKernels::EmbeddedRungeKuttaKernel rkKernel( subRegion, frictionLaw, m_butcherTable );
+ arrayView3d< real64 > const rkStageRates = subRegion.getField< rateAndState::rungeKuttaStageRates >();
+
+ forAll< parallelDevicePolicy<> >( subRegion.size(), [=] GEOS_HOST_DEVICE ( localIndex const k )
+ {
+ rkKernel.updateStageRates( k, stageIndex );
+ rkKernel.updateStageValues( k, stageIndex+1, dt );
+ } );
+ } );
+ } );
+}
+
+void QuasiDynamicEQRK32::stepRateStateODEAndComputeError( real64 const dt, DomainPartition & domain ) const
+{
+ forDiscretizationOnMeshTargets( domain.getMeshBodies(), [&]( string const &,
+ MeshLevel & mesh,
+ arrayView1d< string const > const & regionNames )
+
+ {
+ mesh.getElemManager().forElementSubRegions< SurfaceElementSubRegion >( regionNames,
+ [&]( localIndex const,
+ SurfaceElementSubRegion & subRegion )
+ {
+
+ string const & fricitonLawName = subRegion.template getReference< string >( viewKeyStruct::frictionLawNameString() );
+ RateAndStateFriction const & frictionLaw = getConstitutiveModel< RateAndStateFriction >( subRegion, fricitonLawName );
+ rateAndStateKernels::EmbeddedRungeKuttaKernel rkKernel( subRegion, frictionLaw, m_butcherTable );
+ arrayView3d< real64 > const rkStageRates = subRegion.getField< rateAndState::rungeKuttaStageRates >();
+ if( m_butcherTable.FSAL )
+ {
+ forAll< parallelDevicePolicy<> >( subRegion.size(), [=] GEOS_HOST_DEVICE ( localIndex const k )
+ {
+ // Perform last stage rate update
+ rkKernel.updateStageRates( k, m_butcherTable.numStages-1 );
+ // Update solution to final time and compute errors
+ rkKernel.updateSolutionAndLocalErrorFSAL( k, dt, m_controller.absTol, m_controller.relTol );
+ } );
+ }
+ else
+ {
+ forAll< parallelDevicePolicy<> >( subRegion.size(), [=] GEOS_HOST_DEVICE ( localIndex const k )
+ {
+ // Perform last stage rate update
+ rkKernel.updateStageRates( k, m_butcherTable.numStages-1 );
+ // Update solution to final time and compute errors
+ rkKernel.updateSolutionAndLocalError( k, dt, m_controller.absTol, m_controller.relTol );
+ } );
+ }
+ } );
+ } );
+}
+
+real64 QuasiDynamicEQRK32::updateStresses( real64 const & time_n,
+ real64 const & dt,
+ const int cycleNumber,
+ DomainPartition & domain ) const
+{
+ GEOS_LOG_LEVEL_RANK_0( 1, "Stress solver" );
+ // Call member variable stress solver to update the stress state
+ if( m_stressSolver )
+ {
+ // 1. Solve the momentum balance
+ real64 const dtStress = m_stressSolver->solverStep( time_n, dt, cycleNumber, domain );
+
+ return dtStress;
+ }
+ else
+ {
+ // Spring-slider shear traction computation
+ forDiscretizationOnMeshTargets( domain.getMeshBodies(), [&]( string const &,
+ MeshLevel & mesh,
+ arrayView1d< string const > const & regionNames )
+
+ {
+ mesh.getElemManager().forElementSubRegions< SurfaceElementSubRegion >( regionNames,
+ [&]( localIndex const,
+ SurfaceElementSubRegion & subRegion )
+ {
+
+ arrayView2d< real64 const > const deltaSlip = subRegion.getField< rateAndState::deltaSlip >();
+ arrayView2d< real64 > const traction = subRegion.getField< fields::contact::traction >();
+ arrayView2d< real64 const > const traction_n = subRegion.getField< fields::contact::traction_n >();
+
+ string const & fricitonLawName = subRegion.template getReference< string >( viewKeyStruct::frictionLawNameString() );
+ RateAndStateFriction const & frictionLaw = getConstitutiveModel< RateAndStateFriction >( subRegion, fricitonLawName );
+
+ RateAndStateFriction::KernelWrapper frictionKernelWrapper = frictionLaw.createKernelUpdates();
+
+ forAll< parallelDevicePolicy<> >( subRegion.size(), [=] GEOS_HOST_DEVICE ( localIndex const k )
+ {
+ SpringSliderParameters springSliderParameters = SpringSliderParameters( traction[k][0],
+ frictionKernelWrapper.getACoefficient( k ),
+ frictionKernelWrapper.getBCoefficient( k ),
+ frictionKernelWrapper.getDcCoefficient( k ) );
+
+
+ traction[k][1] = traction_n[k][1] + springSliderParameters.tauRate * dt
+ - springSliderParameters.springStiffness * deltaSlip[k][0];
+ traction[k][2] = traction_n[k][2] + springSliderParameters.tauRate * dt
+ - springSliderParameters.springStiffness * deltaSlip[k][1];
+ } );
+ } );
+ } );
+ return dt;
+ }
+}
+
+void QuasiDynamicEQRK32::updateSlipVelocity( real64 const & time_n,
+ real64 const & dt,
+ DomainPartition & domain ) const
+{
+ GEOS_LOG_LEVEL_RANK_0( 1, "Rate and State solver" );
+ 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 )
+
+ {
+ mesh.getElemManager().forElementSubRegions< SurfaceElementSubRegion >( regionNames,
+ [&]( localIndex const,
+ SurfaceElementSubRegion & subRegion )
+ {
+ // solve rate and state equations.
+ rateAndStateKernels::createAndLaunch< rateAndStateKernels::ExplicitRateAndStateKernel, parallelDevicePolicy<> >( subRegion, viewKeyStruct::frictionLawNameString(), m_shearImpedance,
+ maxIterNewton, newtonTol, time_n, dt );
+ } );
+ } );
+}
+
+void QuasiDynamicEQRK32::saveState( DomainPartition & domain ) const
+{
+ forDiscretizationOnMeshTargets( domain.getMeshBodies(), [&]( string const &,
+ MeshLevel & mesh,
+ arrayView1d< string const > const & regionNames )
+
+ {
+ mesh.getElemManager().forElementSubRegions< SurfaceElementSubRegion >( regionNames,
+ [&]( localIndex const,
+ SurfaceElementSubRegion & subRegion )
+ {
+ arrayView1d< real64 const > const stateVariable = subRegion.getField< rateAndState::stateVariable >();
+ arrayView2d< real64 const > const slipVelocity = subRegion.getField< rateAndState::slipVelocity >();
+ arrayView2d< real64 const > const deltaSlip = subRegion.getField< rateAndState::deltaSlip >();
+ arrayView2d< real64 const > const dispJump = subRegion.getField< contact::dispJump >();
+ arrayView2d< real64 const > const traction = subRegion.getField< contact::traction >();
+
+ arrayView1d< real64 > const stateVariable_n = subRegion.getField< rateAndState::stateVariable_n >();
+ arrayView2d< real64 > const slipVelocity_n = subRegion.getField< rateAndState::slipVelocity_n >();
+ arrayView2d< real64 > const deltaSlip_n = subRegion.getField< rateAndState::deltaSlip >();
+ arrayView2d< real64 > const dispJump_n = subRegion.getField< contact::dispJump_n >();
+ arrayView2d< real64 > const traction_n = subRegion.getField< contact::traction_n >();
+
+ forAll< parallelDevicePolicy<> >( subRegion.size(), [=] GEOS_HOST_DEVICE ( localIndex const k )
+ {
+ stateVariable_n[k] = stateVariable[k];
+ LvArray::tensorOps::copy< 2 >( deltaSlip_n[k], deltaSlip[k] );
+ LvArray::tensorOps::copy< 2 >( slipVelocity_n[k], slipVelocity[k] );
+ LvArray::tensorOps::copy< 3 >( dispJump_n[k], dispJump[k] );
+ LvArray::tensorOps::copy< 3 >( traction_n[k], traction[k] );
+ } );
+ } );
+ } );
+}
+
+real64 QuasiDynamicEQRK32::setNextDt( real64 const & currentDt, DomainPartition & domain )
+{
+
+ // Spring-slider shear traction computation
+ forDiscretizationOnMeshTargets( domain.getMeshBodies(), [&]( string const &,
+ MeshLevel const & mesh,
+ arrayView1d< string const > const & regionNames )
+
+ {
+ mesh.getElemManager().forElementSubRegions< SurfaceElementSubRegion >( regionNames,
+ [&]( localIndex const,
+ SurfaceElementSubRegion const & subRegion )
+ {
+ arrayView2d< real64 const > const error = subRegion.getField< rateAndState::error >();
+
+ RAJA::ReduceSum< parallelDeviceReduce, real64 > scaledl2ErrorSquared( 0.0 );
+ integer const N = subRegion.size();
+ forAll< parallelDevicePolicy<> >( subRegion.size(), [=] GEOS_HOST_DEVICE ( localIndex const k )
+ {
+ scaledl2ErrorSquared += LvArray::tensorOps::l2NormSquared< 3 >( error[k] );
+ } );
+ m_controller.errors[0] = LvArray::math::sqrt( MpiWrapper::sum( scaledl2ErrorSquared.get() / (3.0*N) ));
+ } );
+ } );
+
+ // Compute update factor to currentDt using PID error controller + limiter
+ real64 const dtFactor = m_controller.computeUpdateFactor( m_butcherTable.algHighOrder, m_butcherTable.algLowOrder );
+ real64 const nextDt = dtFactor*currentDt;
+ // Check if step was acceptable
+ m_successfulStep = (dtFactor >= m_controller.acceptSafety) ? true : false;
+ if( m_successfulStep )
+ {
+ m_controller.errors[2] = m_controller.errors[1];
+ m_controller.errors[1] = m_controller.errors[0];
+ GEOS_LOG_LEVEL_RANK_0( 1, GEOS_FMT( "Adaptive time step successful. The next dt will be {:.2e} s", nextDt ));
+ }
+ else
+ {
+ GEOS_LOG_LEVEL_RANK_0( 1, GEOS_FMT( "Adaptive time step failed. The next dt will be {:.2e} s", nextDt ));
+ }
+
+ return nextDt;
+}
+
+REGISTER_CATALOG_ENTRY( PhysicsSolverBase, QuasiDynamicEQRK32, string const &, dataRepository::Group * const )
+
+} // namespace geos
diff --git a/src/coreComponents/physicsSolvers/inducedSeismicity/QuasiDynamicEQRK32.hpp b/src/coreComponents/physicsSolvers/inducedSeismicity/QuasiDynamicEQRK32.hpp
new file mode 100644
index 0000000000..0979fd22dd
--- /dev/null
+++ b/src/coreComponents/physicsSolvers/inducedSeismicity/QuasiDynamicEQRK32.hpp
@@ -0,0 +1,235 @@
+/*
+ * ------------------------------------------------------------------------------------------------------------
+ * SPDX-License-Identifier: LGPL-2.1-only
+ *
+ * Copyright (c) 2016-2024 Lawrence Livermore National Security LLC
+ * Copyright (c) 2018-2024 Total, S.A
+ * 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.
+ * ------------------------------------------------------------------------------------------------------------
+ */
+
+#ifndef GEOS_PHYSICSSOLVERS_INDUCED_QUASIDYNAMICEQRK32_HPP
+#define GEOS_PHYSICSSOLVERS_INDUCED_QUASIDYNAMICEQRK32_HPP
+
+#include "physicsSolvers/PhysicsSolverBase.hpp"
+#include "kernels/RateAndStateKernels.hpp"
+
+namespace geos
+{
+
+class QuasiDynamicEQRK32 : public PhysicsSolverBase
+{
+public:
+ /// The default nullary constructor is disabled to avoid compiler auto-generation:
+ QuasiDynamicEQRK32() = delete;
+
+ /// The constructor needs a user-defined "name" and a parent Group (to place this instance in the tree structure of classes)
+ QuasiDynamicEQRK32( const string & name,
+ Group * const parent );
+
+ /// Destructor
+ virtual ~QuasiDynamicEQRK32() override;
+
+ static string catalogName() { return "QuasiDynamicEQRK32"; }
+
+ /**
+ * @return Get the final class Catalog name
+ */
+ virtual string getCatalogName() const override { return catalogName(); }
+
+ /// This method ties properties with their supporting mesh
+ virtual void registerDataOnMesh( Group & meshBodies ) override;
+
+ struct viewKeyStruct : public PhysicsSolverBase::viewKeyStruct
+ {
+ /// stress solver name
+ constexpr static char const * stressSolverNameString() { return "stressSolverName"; }
+ /// Friction law name string
+ constexpr static char const * frictionLawNameString() { return "frictionLawName"; }
+ /// Friction law name string
+ constexpr static char const * shearImpedanceString() { return "shearImpedance"; }
+ /// target slip increment
+ constexpr static char const * timeStepTol() { return "timeStepTol"; }
+ };
+
+ virtual real64 solverStep( real64 const & time_n,
+ real64 const & dt,
+ integer const cycleNumber,
+ DomainPartition & domain ) override final;
+
+
+ virtual real64 setNextDt( real64 const & currentDt,
+ DomainPartition & domain ) override final;
+
+ /**
+ * @brief Computes stage rates for the initial Runge-Kutta substage and updates slip and state
+ * @param dt
+ * @param domain
+ */
+ void stepRateStateODEInitialSubstage( real64 const dt, DomainPartition & domain ) const;
+
+ /**
+ * @brief Computes stage rates at the Runge-Kutta substage specified by stageIndex and updates slip and state
+ * @param stageIndex
+ * @param dt
+ * @param domain
+ */
+ void stepRateStateODESubstage( integer const stageIndex,
+ real64 const dt,
+ DomainPartition & domain ) const;
+
+ /**
+ * @brief Updates slip and state to t + dt and approximates the error
+ * @param dt
+ * @param domain
+ */
+ void stepRateStateODEAndComputeError( real64 const dt, DomainPartition & domain ) const;
+
+ real64 updateStresses( real64 const & time_n,
+ real64 const & dt,
+ const int cycleNumber,
+ DomainPartition & domain ) const;
+
+ /**
+ * @brief Updates rate-and-state slip velocity
+ * @param domain
+ */
+ void updateSlipVelocity( real64 const & time_n,
+ real64 const & dt,
+ DomainPartition & domain ) const;
+
+ /**
+ * @brief save the current state
+ * @param domain
+ */
+ void saveState( DomainPartition & domain ) const;
+
+private:
+
+ virtual void postInputInitialization() override;
+
+
+ /// pointer to stress solver
+ PhysicsSolverBase * m_stressSolver;
+
+ /// stress solver name
+ string m_stressSolverName;
+
+ /// shear impedance
+ real64 m_shearImpedance;
+
+ /// Runge-Kutta Butcher table (specifies the embedded RK method)
+ // TODO: The specific type should not be hardcoded!
+ // Should be possible to change RK-method based on the table.
+ rateAndStateKernels::BogackiShampine32Table m_butcherTable;
+
+ bool m_successfulStep; // Flag indicating if the adative time step was accepted
+
+ /**
+ * @brief Proportional-integral-derivative controller used for updating time step
+ * based error estimate in the current and previous time steps.
+ */
+ class PIDController
+ {
+public:
+
+ GEOS_HOST_DEVICE
+ PIDController( std::array< const real64, 3 > const & cparams,
+ const real64 atol,
+ const real64 rtol,
+ const real64 safety ):
+ controlParameters{ cparams },
+ absTol( atol ),
+ relTol( rtol ),
+ acceptSafety( safety ),
+ errors{ {0.0, 0.0, 0.0} }
+ {}
+
+ /// Default copy constructor
+ PIDController( PIDController const & ) = default;
+
+ /// Default move constructor
+ PIDController( PIDController && ) = default;
+
+ /// Deleted default constructor
+ PIDController() = delete;
+
+ /// Deleted copy assignment operator
+ PIDController & operator=( PIDController const & ) = delete;
+
+ /// Deleted move assignment operator
+ PIDController & operator=( PIDController && ) = delete;
+
+ /// Parameters for the PID error controller
+ const std::array< const real64, 3 > controlParameters; // Controller parameters
+
+ real64 const absTol; // absolut tolerence
+
+ real64 const relTol; // relative tolerence
+
+ real64 const acceptSafety; // Acceptance safety
+
+ std::array< real64, 3 > errors; // Errors for current and two previous updates
+ // stored as [n+1, n, n-1]
+
+ real64 computeUpdateFactor( integer const algHighOrder, integer const algLowOrder )
+ {
+ // PID error controller + limiter
+ real64 const k = LvArray::math::min( algHighOrder, algLowOrder ) + 1.0;
+ real64 const eps0 = 1.0/(errors[0] + std::numeric_limits< real64 >::epsilon()); // n + 1
+ real64 const eps1 = 1.0/(errors[1] + std::numeric_limits< real64 >::epsilon()); // n
+ real64 const eps2 = 1.0/(errors[2] + std::numeric_limits< real64 >::epsilon()); // n-1
+ // Compute update factor eps0^(beta0/k)*eps1^(beta1/k)*eps2^(beta2/k) where
+ // beta0 - beta2 are the control parameters. Also apply limiter to smoothen changes.
+ // Limiter is 1.0 + atan(x - 1.0). Here use atan(x) = atan2(x, 1.0).
+ return 1.0 + LvArray::math::atan2( pow( eps0, controlParameters[0] / k ) *
+ pow( eps1, controlParameters[1] / k ) *
+ pow( eps2, controlParameters[2] / k ) - 1.0, 1.0 );
+ }
+ };
+
+ PIDController m_controller;
+
+
+ class SpringSliderParameters
+ {
+public:
+
+ GEOS_HOST_DEVICE
+ SpringSliderParameters( real64 const normalTraction, real64 const a, real64 const b, real64 const Dc ):
+ tauRate( 1e-4 ),
+ springStiffness( 0.0 )
+ {
+ real64 const criticalStiffness = normalTraction * (b - a) / Dc;
+ springStiffness = 0.9 * criticalStiffness;
+ }
+
+ /// Default copy constructor
+ SpringSliderParameters( SpringSliderParameters const & ) = default;
+
+ /// Default move constructor
+ SpringSliderParameters( SpringSliderParameters && ) = default;
+
+ /// Deleted default constructor
+ SpringSliderParameters() = delete;
+
+ /// Deleted copy assignment operator
+ SpringSliderParameters & operator=( SpringSliderParameters const & ) = delete;
+
+ /// Deleted move assignment operator
+ SpringSliderParameters & operator=( SpringSliderParameters && ) = delete;
+
+ real64 tauRate;
+
+ real64 springStiffness;
+ };
+};
+
+} /* namespace geos */
+
+#endif /* GEOS_PHYSICSSOLVERS_INDUCED_QUASIDYNAMICEQRK32_HPP */
diff --git a/src/coreComponents/physicsSolvers/inducedSeismicity/kernels/RateAndStateKernels.hpp b/src/coreComponents/physicsSolvers/inducedSeismicity/kernels/RateAndStateKernels.hpp
index 54e50af7ad..d91dff75ef 100644
--- a/src/coreComponents/physicsSolvers/inducedSeismicity/kernels/RateAndStateKernels.hpp
+++ b/src/coreComponents/physicsSolvers/inducedSeismicity/kernels/RateAndStateKernels.hpp
@@ -27,20 +27,38 @@ namespace geos
namespace rateAndStateKernels
{
+
+// TBD: Pass the kernel and add getters for relevant fields to make this function general purpose and avoid
+// wrappers?
+GEOS_HOST_DEVICE
+static void projectSlipRateBase( localIndex const k,
+ real64 const frictionCoefficient,
+ real64 const shearImpedance,
+ arrayView2d< real64 const > const traction,
+ arrayView1d< real64 const > const slipRate,
+ arrayView2d< real64 > const slipVelocity )
+{
+ // Project slip rate onto shear traction to get slip velocity components
+ real64 const frictionForce = traction[k][0] * frictionCoefficient;
+ real64 const projectionScaling = 1.0 / ( shearImpedance + frictionForce / slipRate[k] );
+ slipVelocity[k][0] = projectionScaling * traction[k][1];
+ slipVelocity[k][1] = projectionScaling * traction[k][2];
+}
+
/**
- * @class RateAndStateKernel
+ * @class ImplicitFixedStressRateAndStateKernel
*
* @brief
*
* @details
*/
-class RateAndStateKernel
+class ImplicitFixedStressRateAndStateKernel
{
public:
- RateAndStateKernel( SurfaceElementSubRegion & subRegion,
- constitutive::RateAndStateFriction const & frictionLaw,
- real64 const shearImpedance ):
+ ImplicitFixedStressRateAndStateKernel( SurfaceElementSubRegion & subRegion,
+ constitutive::RateAndStateFriction const & frictionLaw,
+ real64 const shearImpedance ):
m_slipRate( subRegion.getField< fields::rateAndState::slipRate >() ),
m_stateVariable( subRegion.getField< fields::rateAndState::stateVariable >() ),
m_stateVariable_n( subRegion.getField< fields::rateAndState::stateVariable_n >() ),
@@ -90,7 +108,7 @@ class RateAndStateKernel
stack.jacobian[0][0] = dFriction[0]; // derivative of Eq 1 w.r.t. stateVariable
stack.jacobian[0][1] = dFriction[1]; // derivative of Eq 1 w.r.t. slipRate
stack.jacobian[1][0] = dStateEvolutionLaw[0]; // derivative of Eq 2 w.r.t. stateVariable
- stack.jacobian[1][1] = dStateEvolutionLaw[1]; // derivative of Eq 2 w.r.t. m_slipRate
+ stack.jacobian[1][1] = dStateEvolutionLaw[1]; // derivative of Eq 2 w.r.t. slipRate
}
GEOS_HOST_DEVICE
@@ -109,11 +127,8 @@ class RateAndStateKernel
GEOS_HOST_DEVICE
void projectSlipRate( localIndex const k ) const
{
- // Project slip rate onto shear traction to get slip velocity components
- real64 const frictionForce = m_traction[k][0] * m_frictionLaw.frictionCoefficient( k, m_slipRate[k], m_stateVariable[k] );
- real64 const projectionScaling = 1.0 / ( m_shearImpedance + frictionForce / m_slipRate[k] );
- m_slipVelocity[k][0] = projectionScaling * m_traction[k][1];
- m_slipVelocity[k][1] = projectionScaling * m_traction[k][2];
+ real64 const frictionCoefficient = m_frictionLaw.frictionCoefficient( k, m_slipRate[k], m_stateVariable[k] );
+ projectSlipRateBase( k, frictionCoefficient, m_shearImpedance, m_traction, m_slipRate, m_slipVelocity );
}
GEOS_HOST_DEVICE
@@ -144,18 +159,124 @@ class RateAndStateKernel
};
+/**
+ * @class ExplicitRateAndStateKernel
+ *
+ * @brief
+ *
+ * @details
+ */
+class ExplicitRateAndStateKernel
+{
+public:
+
+ ExplicitRateAndStateKernel( SurfaceElementSubRegion & subRegion,
+ constitutive::RateAndStateFriction const & frictionLaw,
+ real64 const shearImpedance ):
+ m_slipRate( subRegion.getField< fields::rateAndState::slipRate >() ),
+ m_stateVariable( subRegion.getField< fields::rateAndState::stateVariable >() ),
+ m_traction( subRegion.getField< fields::contact::traction >() ),
+ m_slipVelocity( subRegion.getField< fields::rateAndState::slipVelocity >() ),
+ m_shearImpedance( shearImpedance ),
+ m_frictionLaw( frictionLaw.createKernelUpdates() )
+ {}
+
+ /**
+ * @struct StackVariables
+ * @brief Kernel variables located on the stack
+ */
+ struct StackVariables
+ {
+public:
+
+ GEOS_HOST_DEVICE
+ StackVariables( )
+ {}
+
+ real64 jacobian;
+ real64 rhs;
+
+ };
+
+ GEOS_HOST_DEVICE
+ void setup( localIndex const k,
+ real64 const dt,
+ StackVariables & stack ) const
+ {
+ GEOS_UNUSED_VAR( dt );
+ real64 const normalTraction = m_traction[k][0];
+ real64 const shearTractionMagnitude = LvArray::math::sqrt( m_traction[k][1] * m_traction[k][1] + m_traction[k][2] * m_traction[k][2] );
+
+ // Slip rate is bracketed between [0, shear traction magnitude / shear impedance]
+ // If slip rate is outside the bracket, re-initialize to the middle value
+ real64 const upperBound = shearTractionMagnitude/m_shearImpedance;
+ real64 const bracketedSlipRate = m_slipRate[k] > upperBound ? 0.5*upperBound : m_slipRate[k];
+
+ stack.rhs = shearTractionMagnitude - m_shearImpedance *bracketedSlipRate - normalTraction * m_frictionLaw.frictionCoefficient( k, bracketedSlipRate, m_stateVariable[k] );
+ stack.jacobian = -m_shearImpedance - normalTraction * m_frictionLaw.dFrictionCoefficient_dSlipRate( k, bracketedSlipRate, m_stateVariable[k] );
+ }
+
+ GEOS_HOST_DEVICE
+ void solve( localIndex const k,
+ StackVariables & stack ) const
+ {
+ m_slipRate[k] -= stack.rhs/stack.jacobian;
+
+ // Slip rate is bracketed between [0, shear traction magnitude / shear impedance]
+ // Check that the update did not end outside of the bracket.
+ real64 const shearTractionMagnitude = LvArray::math::sqrt( m_traction[k][1] * m_traction[k][1] + m_traction[k][2] * m_traction[k][2] );
+ real64 const upperBound = shearTractionMagnitude/m_shearImpedance;
+ if( m_slipRate[k] > upperBound ) m_slipRate[k] = 0.5*upperBound;
+
+ }
+
+
+ GEOS_HOST_DEVICE
+ camp::tuple< int, real64 > checkConvergence( StackVariables const & stack,
+ real64 const tol ) const
+ {
+ real64 const residualNorm = LvArray::math::abs( stack.rhs );
+ int const converged = residualNorm < tol ? 1 : 0;
+ camp::tuple< int, real64 > result { converged, residualNorm };
+ return result;
+ }
+
+ GEOS_HOST_DEVICE
+ void projectSlipRate( localIndex const k ) const
+ {
+ real64 const frictionCoefficient = m_frictionLaw.frictionCoefficient( k, m_slipRate[k], m_stateVariable[k] );
+ projectSlipRateBase( k, frictionCoefficient, m_shearImpedance, m_traction, m_slipRate, m_slipVelocity );
+ }
+
+private:
+
+ arrayView1d< real64 > const m_slipRate;
+
+ arrayView1d< real64 > const m_stateVariable;
+
+ arrayView2d< real64 const > const m_traction;
+
+ arrayView2d< real64 > const m_slipVelocity;
+
+ real64 const m_shearImpedance;
+
+ constitutive::RateAndStateFriction::KernelWrapper m_frictionLaw;
+
+};
/**
* @brief Performs the kernel launch
+ * @tparam KernelType The Rate-and-state kernel to launch
* @tparam POLICY the policy used in the RAJA kernels
*/
-template< typename POLICY >
+template< typename KernelType, typename POLICY >
static void
createAndLaunch( SurfaceElementSubRegion & subRegion,
string const & frictionLawNameKey,
real64 const shearImpedance,
- integer const maxNewtonIter,
+ integer const maxIterNewton,
+ real64 const newtonTol,
real64 const time_n,
real64 const dt )
{
@@ -165,20 +286,20 @@ createAndLaunch( SurfaceElementSubRegion & subRegion,
string const & frictionaLawName = subRegion.getReference< string >( frictionLawNameKey );
constitutive::RateAndStateFriction const & frictionLaw = subRegion.getConstitutiveModel< constitutive::RateAndStateFriction >( frictionaLawName );
- RateAndStateKernel kernel( subRegion, frictionLaw, shearImpedance );
+ KernelType kernel( subRegion, frictionLaw, shearImpedance );
// Newton loop (outside of the kernel launch)
bool allConverged = false;
- for( integer iter = 0; iter < maxNewtonIter; iter++ )
+ for( integer iter = 0; iter < maxIterNewton; iter++ )
{
RAJA::ReduceMin< parallelDeviceReduce, int > converged( 1 );
RAJA::ReduceMax< parallelDeviceReduce, real64 > residualNorm( 0.0 );
forAll< POLICY >( subRegion.size(), [=] GEOS_HOST_DEVICE ( localIndex const k )
{
- RateAndStateKernel::StackVariables stack;
+ typename KernelType::StackVariables stack;
kernel.setup( k, dt, stack );
kernel.solve( k, stack );
- auto const [elementConverged, elementResidualNorm] = kernel.checkConvergence( stack, 1.0e-6 );
+ auto const [elementConverged, elementResidualNorm] = kernel.checkConvergence( stack, newtonTol );
converged.min( elementConverged );
residualNorm.max( elementResidualNorm );
} );
@@ -202,8 +323,260 @@ createAndLaunch( SurfaceElementSubRegion & subRegion,
} );
}
+/**
+ * @brief Butcher table for embedded RK3(2) method using Kuttas third order
+ * method for the high-order update, and an explicit trapezoidal rule
+ * based on the first and third stage rates for the low-order update.
+ */
+struct Kutta32Table
+{
+ integer constexpr static algHighOrder = 3; // High-order update order
+ integer constexpr static algLowOrder = 2; // Low-order update order
+ integer constexpr static numStages = 3; // Number of stages
+ real64 const a[2][2] = { { 1.0/2.0, 0.0 }, // Coefficients for stage value updates
+ { -1.0, 2.0 } }; // (lower-triangular part of table).
+ real64 const c[3] = { 0.0, 1.0/2.0, 1.0 }; // Coefficients for time increments of substages
+ real64 const b[3] = { 1.0/6.0, 4.0/6.0, 1.0/6.0 }; // Quadrature weights used to step the solution to next time
+ real64 const bStar[3] = { 1.0/2.0, 0.0, 1.0/2.0 }; // Quadrature weights used for low-order comparision solution
+ real64 constexpr static FSAL = false; // Not first same as last
+};
+
+/**
+ * @brief Butcher table for the BogackiShampine 3(2) method.
+ */
+struct BogackiShampine32Table
+{
+ integer constexpr static algHighOrder = 3; // High-order update order
+ integer constexpr static algLowOrder = 2; // Low-order update order
+ integer constexpr static numStages = 4; // Number of stages
+ real64 const a[3][3] = { { 1.0/2.0, 0.0, 0.0 }, // Coefficients for stage value updates
+ { 0.0, 3.0/4.0, 0.0 }, // (lower-triangular part of table).
+ { 2.0/9.0, 1.0/3.0, 4.0/9.0 } };
+ real64 const c[4] = { 0.0, 1.0/2.0, 3.0/4.0, 1.0 }; // Coefficients for time increments of substages
+ real64 const b[4] = { 2.0/9.0, 1.0/3.0, 4.0/9.0, 0.0 }; // Quadrature weights used to step the solution to next time
+ real64 const bStar[4] = { 7.0/24.0, 1.0/4.0, 1.0/3.0, 1.0/8.0}; // Quadrature weights used for low-order comparision solution
+ bool constexpr static FSAL = true; // First same as last (can reuse the last stage rate in next
+ // update)
+};
+
+/**
+ * @brief Runge-Kutta method used to time integrate slip and state. Uses of a high order
+ * update used to integrate the solutions, and a lower order update to estimate the error
+ * in the time step.
+ *
+ * @tparam Butcher table defining the Runge-Kutta method.
+ */
+template< typename TABLE_TYPE > class EmbeddedRungeKuttaKernel
+{
+
+public:
+ EmbeddedRungeKuttaKernel( SurfaceElementSubRegion & subRegion,
+ constitutive::RateAndStateFriction const & frictionLaw,
+ TABLE_TYPE butcherTable ):
+ m_stateVariable( subRegion.getField< fields::rateAndState::stateVariable >() ),
+ m_stateVariable_n( subRegion.getField< fields::rateAndState::stateVariable_n >() ),
+ m_slipRate( subRegion.getField< fields::rateAndState::slipRate >() ),
+ m_slipVelocity( subRegion.getField< fields::rateAndState::slipVelocity >() ),
+ m_slipVelocity_n( subRegion.getField< fields::rateAndState::slipVelocity_n >() ),
+ m_deltaSlip( subRegion.getField< fields::rateAndState::deltaSlip >() ),
+ m_deltaSlip_n( subRegion.getField< fields::rateAndState::deltaSlip_n >() ),
+ m_dispJump( subRegion.getField< fields::contact::dispJump >() ),
+ m_dispJump_n( subRegion.getField< fields::contact::dispJump_n >() ),
+ m_error( subRegion.getField< fields::rateAndState::error >() ),
+ m_stageRates( subRegion.getField< fields::rateAndState::rungeKuttaStageRates >() ),
+ m_frictionLaw( frictionLaw.createKernelUpdates() ),
+ m_butcherTable( butcherTable )
+ {}
+
+ /**
+ * @brief Initialize slip and state buffers
+ */
+ GEOS_HOST_DEVICE
+ void initialize( localIndex const k ) const
+ {
+ LvArray::tensorOps::copy< 2 >( m_slipVelocity[k], m_slipVelocity_n[k] );
+ m_slipRate[k] = LvArray::tensorOps::l2Norm< 2 >( m_slipVelocity_n[k] );
+ m_stateVariable[k] = m_stateVariable_n[k];
+ }
+
+ /**
+ * @brief Re-uses the last stage rate from the previous time step as the first
+ * in the next update. Only valid for FSAL (first-same-as-last) Runge-Kutta methods.
+ */
+ GEOS_HOST_DEVICE
+ void updateStageRatesFSAL( localIndex const k ) const
+ {
+ LvArray::tensorOps::copy< 3 >( m_stageRates[k][0], m_stageRates[k][m_butcherTable.numStages-1] );
+ }
+
+ /**
+ * @brief Updates the stage rates rates (the right-hand-side of the ODEs for slip and state)
+ */
+ GEOS_HOST_DEVICE
+ void updateStageRates( localIndex const k, integer const stageIndex ) const
+ {
+ m_stageRates[k][stageIndex][0] = m_slipVelocity[k][0];
+ m_stageRates[k][stageIndex][1] = m_slipVelocity[k][1];
+ m_stageRates[k][stageIndex][2] = m_frictionLaw.stateEvolution( k, m_slipRate[k], m_stateVariable[k] );
+ }
+
+ /**
+ * @brief Update stage values (slip, state and displacement jump) to a Runge-Kutta substage.
+ */
+ GEOS_HOST_DEVICE
+ void updateStageValues( localIndex const k, integer const stageIndex, real64 const dt ) const
+ {
+
+ real64 stateVariableIncrement = 0.0;
+ real64 deltaSlipIncrement[2] = {0.0, 0.0};
+
+ for( integer i = 0; i < stageIndex; i++ )
+ {
+ deltaSlipIncrement[0] += m_butcherTable.a[stageIndex-1][i] * m_stageRates[k][i][0];
+ deltaSlipIncrement[1] += m_butcherTable.a[stageIndex-1][i] * m_stageRates[k][i][1];
+ stateVariableIncrement += m_butcherTable.a[stageIndex-1][i] * m_stageRates[k][i][2];
+ }
+ m_deltaSlip[k][0] = m_deltaSlip_n[k][0] + dt*deltaSlipIncrement[0];
+ m_deltaSlip[k][1] = m_deltaSlip_n[k][1] + dt*deltaSlipIncrement[1];
+ m_stateVariable[k] = m_stateVariable_n[k] + dt*stateVariableIncrement;
+
+ m_dispJump[k][1] = m_dispJump_n[k][1] + m_deltaSlip[k][0];
+ m_dispJump[k][2] = m_dispJump_n[k][2] + m_deltaSlip[k][1];
+ }
+
+ /**
+ * @brief Updates slip, state and displacement jump to the next time computes error the local error
+ * in the time step
+ */
+ GEOS_HOST_DEVICE
+ void updateSolutionAndLocalError( localIndex const k, real64 const dt, real64 const absTol, real64 const relTol ) const
+ {
+
+ real64 deltaSlipIncrement[2] = {0.0, 0.0};
+ real64 deltaSlipIncrementLowOrder[2] = {0.0, 0.0};
+
+ real64 stateVariableIncrement = 0.0;
+ real64 stateVariableIncrementLowOrder = 0.0;
+
+ for( localIndex i = 0; i < m_butcherTable.numStages; i++ )
+ {
+
+ // High order update of solution
+ deltaSlipIncrement[0] += m_butcherTable.b[i] * m_stageRates[k][i][0];
+ deltaSlipIncrement[1] += m_butcherTable.b[i] * m_stageRates[k][i][1];
+ stateVariableIncrement += m_butcherTable.b[i] * m_stageRates[k][i][2];
+
+ // Low order update for error
+ deltaSlipIncrementLowOrder[0] += m_butcherTable.bStar[i] * m_stageRates[k][i][0];
+ deltaSlipIncrementLowOrder[1] += m_butcherTable.bStar[i] * m_stageRates[k][i][1];
+ stateVariableIncrementLowOrder += m_butcherTable.bStar[i] * m_stageRates[k][i][2];
+ }
+
+ m_deltaSlip[k][0] = m_deltaSlip_n[k][0] + dt * deltaSlipIncrement[0];
+ m_deltaSlip[k][1] = m_deltaSlip_n[k][1] + dt * deltaSlipIncrement[1];
+ m_stateVariable[k] = m_stateVariable_n[k] + dt * stateVariableIncrement;
+
+ real64 const deltaSlipLowOrder[2] = {m_deltaSlip_n[k][0] + dt * deltaSlipIncrementLowOrder[0],
+ m_deltaSlip_n[k][1] + dt * deltaSlipIncrementLowOrder[1]};
+ real64 const stateVariableLowOrder = m_stateVariable_n[k] + dt * stateVariableIncrementLowOrder;
+
+ m_dispJump[k][1] = m_dispJump_n[k][1] + m_deltaSlip[k][0];
+ m_dispJump[k][2] = m_dispJump_n[k][2] + m_deltaSlip[k][1];
+
+ // Compute error
+ m_error[k][0] = computeError( m_deltaSlip[k][0], deltaSlipLowOrder[0], absTol, relTol );
+ m_error[k][1] = computeError( m_deltaSlip[k][1], deltaSlipLowOrder[1], absTol, relTol );
+ m_error[k][2] = computeError( m_stateVariable[k], stateVariableLowOrder, absTol, relTol );
+ }
+
+ /**
+ * @brief Updates slip, state and displacement jump to the next time computes error the local error
+ * in the time step. Uses the FSAL (first-same-as-last) property.
+ */
+ GEOS_HOST_DEVICE
+ void updateSolutionAndLocalErrorFSAL( localIndex const k, real64 const dt, real64 const absTol, real64 const relTol ) const
+ {
+
+ real64 deltaSlipIncrementLowOrder[2] = {0.0, 0.0};
+ real64 stateVariableIncrementLowOrder = 0.0;
+
+ for( localIndex i = 0; i < m_butcherTable.numStages; i++ )
+ {
+ // In FSAL algorithms the last RK substage update coincides with the
+ // high-order update. Only need to compute increments for the the
+ // low-order updates for error computation.
+ deltaSlipIncrementLowOrder[0] += m_butcherTable.bStar[i] * m_stageRates[k][i][0];
+ deltaSlipIncrementLowOrder[1] += m_butcherTable.bStar[i] * m_stageRates[k][i][1];
+ stateVariableIncrementLowOrder += m_butcherTable.bStar[i] * m_stageRates[k][i][2];
+ }
+
+ real64 const deltaSlipLowOrder[2] = {m_deltaSlip_n[k][0] + dt * deltaSlipIncrementLowOrder[0],
+ m_deltaSlip_n[k][1] + dt * deltaSlipIncrementLowOrder[1]};
+ real64 const stateVariableLowOrder = m_stateVariable_n[k] + dt * stateVariableIncrementLowOrder;
+
+ m_dispJump[k][1] = m_dispJump_n[k][1] + m_deltaSlip[k][0];
+ m_dispJump[k][2] = m_dispJump_n[k][2] + m_deltaSlip[k][1];
+
+ // Compute error
+ m_error[k][0] = computeError( m_deltaSlip[k][0], deltaSlipLowOrder[0], absTol, relTol );
+ m_error[k][1] = computeError( m_deltaSlip[k][1], deltaSlipLowOrder[1], absTol, relTol );
+ m_error[k][2] = computeError( m_stateVariable[k], stateVariableLowOrder, absTol, relTol );
+ }
+
+ /**
+ * @brief Computes the relative error scaled by error tolerances
+ */
+ GEOS_HOST_DEVICE
+ real64 computeError( real64 const highOrderApprox, real64 const lowOrderApprox, real64 const absTol, real64 const relTol ) const
+ {
+ return (highOrderApprox - lowOrderApprox) /
+ ( absTol + relTol * LvArray::math::max( LvArray::math::abs( highOrderApprox ), LvArray::math::abs( lowOrderApprox ) ));
+ }
+
+private:
+
+ /// Current state variable
+ arrayView1d< real64 > const m_stateVariable;
+
+ /// State variable at t = t_n
+ arrayView1d< real64 > const m_stateVariable_n;
+
+ /// Current slip rate (magnitude of slip velocity)
+ arrayView1d< real64 > const m_slipRate;
+
+ /// Current slip velocity
+ arrayView2d< real64 > const m_slipVelocity;
+
+ /// Slip velocity at time t_n
+ arrayView2d< real64 > const m_slipVelocity_n;
+
+ /// Current slip change
+ arrayView2d< real64 > const m_deltaSlip;
+
+ /// Slip change at time t_n
+ arrayView2d< real64 > const m_deltaSlip_n;
+
+ /// Current displacment jump
+ arrayView2d< real64 > const m_dispJump;
+
+ /// Displacment jump at time t_n
+ arrayView2d< real64 > const m_dispJump_n;
+
+ /// Local error for each solution component stored as slip1, slip2, state
+ arrayView2d< real64 > const m_error;
+
+ /// Stage rates for each solution component stored as slip1, slip2, state
+ arrayView3d< real64 > const m_stageRates;
+
+ /// Friction law used for rate-and-state updates
+ constitutive::RateAndStateFriction::KernelWrapper m_frictionLaw;
+
+ /// Butcher table used for explicit time stepping of slip and state
+ TABLE_TYPE m_butcherTable;
+};
+
} /* namespace rateAndStateKernels */
-}/* namespace geos */
+} /* namespace geos */
#endif /* GEOS_PHYSICSSOLVERS_RATEANDSTATEKERNELS_HPP_ */
diff --git a/src/coreComponents/physicsSolvers/inducedSeismicity/rateAndStateFields.hpp b/src/coreComponents/physicsSolvers/inducedSeismicity/rateAndStateFields.hpp
index 0e7e51a805..e060dceb28 100644
--- a/src/coreComponents/physicsSolvers/inducedSeismicity/rateAndStateFields.hpp
+++ b/src/coreComponents/physicsSolvers/inducedSeismicity/rateAndStateFields.hpp
@@ -36,10 +36,26 @@ DECLARE_FIELD( slipRate,
"slipRate",
array1d< real64 >,
1.0e-6,
- LEVEL_0,
+ NOPLOT,
WRITE_AND_READ,
"Slip rate" );
+DECLARE_FIELD( slipVelocity,
+ "slipVelocity",
+ array2d< real64 >,
+ 0.70710678118e-6,
+ LEVEL_0,
+ WRITE_AND_READ,
+ "Slip velocity" );
+
+DECLARE_FIELD( slipVelocity_n,
+ "slipVelocity_n",
+ array2d< real64 >,
+ 0.70710678118e-6,
+ NOPLOT,
+ WRITE_AND_READ,
+ "Slip velocity at previous time step" );
+
DECLARE_FIELD( stateVariable,
"stateVariable",
array1d< real64 >,
@@ -48,21 +64,14 @@ DECLARE_FIELD( stateVariable,
WRITE_AND_READ,
"Rate- and state-dependent friction state variable" );
-DECLARE_FIELD( slipVelocity,
- "slipVelocity",
- array2d< real64 >,
- 1.0e-6,
- LEVEL_0,
- WRITE_AND_READ,
- "Slip velocity" );
-
DECLARE_FIELD( stateVariable_n,
"stateVariable_n",
array1d< real64 >,
0.6,
NOPLOT,
WRITE_AND_READ,
- "Rate- and state-dependent friction state variable at previous time step" );
+ "Initial rate- and state-dependent friction state variable at this time step" );
+
DECLARE_FIELD( deltaSlip,
"deltaSlip",
@@ -72,6 +81,31 @@ DECLARE_FIELD( deltaSlip,
WRITE_AND_READ,
"Slip increment" );
+DECLARE_FIELD( deltaSlip_n,
+ "deltaSlip_n",
+ array2d< real64 >,
+ 0.0,
+ NOPLOT,
+ WRITE_AND_READ,
+ "Initial slip increment at this time step" );
+
+
+DECLARE_FIELD( rungeKuttaStageRates,
+ "rungeKuttaStageRates",
+ array3d< real64 >,
+ 0.0,
+ NOPLOT,
+ WRITE_AND_READ,
+ "Runge-Kutta stage rates for rate-and-state variables" );
+
+
+DECLARE_FIELD( error,
+ "error",
+ array2d< real64 >,
+ 0.0,
+ LEVEL_0,
+ WRITE_AND_READ,
+ "Error for rate-and-state fields" );
}
diff --git a/src/coreComponents/schema/schema.xsd b/src/coreComponents/schema/schema.xsd
index 7a6e7967c0..e130b1e56d 100644
--- a/src/coreComponents/schema/schema.xsd
+++ b/src/coreComponents/schema/schema.xsd
@@ -399,6 +399,10 @@
+
+
+
+
@@ -2109,6 +2113,11 @@ the relative residual norm satisfies:
+
+
@@ -2125,6 +2134,11 @@ the relative residual norm satisfies:
+
+
@@ -2139,6 +2153,11 @@ the relative residual norm satisfies:
+
+
@@ -2147,6 +2166,11 @@ the relative residual norm satisfies:
+
+
@@ -2157,6 +2181,11 @@ the relative residual norm satisfies:
+
+
@@ -2195,7 +2224,7 @@ the relative residual norm satisfies:
-
+
@@ -2262,6 +2291,7 @@ the relative residual norm satisfies:
+
@@ -2647,6 +2677,8 @@ Level 0 outputs no specific information for this solver. Higher levels require m
+
+
@@ -2734,6 +2766,8 @@ Level 0 outputs no specific information for this solver. Higher levels require m
+
+
@@ -3211,7 +3245,7 @@ Level 0 outputs no specific information for this solver. Higher levels require m
-
+
@@ -3234,7 +3268,7 @@ Local- Add jump stabilization on interior of macro elements-->
-
+
@@ -3538,6 +3572,41 @@ Level 0 outputs no specific information for this solver. Higher levels require m
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
diff --git a/src/coreComponents/schema/schema.xsd.other b/src/coreComponents/schema/schema.xsd.other
index 33386f1995..3abf5b863f 100644
--- a/src/coreComponents/schema/schema.xsd.other
+++ b/src/coreComponents/schema/schema.xsd.other
@@ -486,7 +486,7 @@
-
+
@@ -538,6 +538,7 @@
+
@@ -571,7 +572,7 @@
-
+
@@ -608,7 +609,7 @@
-
+
@@ -659,7 +660,7 @@
-
+
@@ -700,7 +701,7 @@
-
+
@@ -733,7 +734,7 @@
-
+
@@ -744,7 +745,7 @@
-
+
@@ -757,7 +758,7 @@
-
+
@@ -770,7 +771,7 @@
-
+
@@ -786,7 +787,7 @@
-
+
@@ -820,7 +821,7 @@
-
+
@@ -883,7 +884,7 @@
-
+
@@ -914,7 +915,7 @@
-
+
@@ -927,7 +928,7 @@
-
+
@@ -940,7 +941,7 @@
-
+
@@ -953,7 +954,7 @@
-
+
@@ -966,7 +967,7 @@
-
+
@@ -981,7 +982,7 @@
-
+
@@ -992,7 +993,7 @@
-
+
@@ -1005,7 +1006,7 @@
-
+
@@ -1016,7 +1017,7 @@
-
+
@@ -1027,7 +1028,18 @@
-
+
+
+
+
+
+
+
+
+
+
+
+
@@ -1038,7 +1050,7 @@
-
+
@@ -1051,7 +1063,7 @@
-
+
@@ -1062,7 +1074,7 @@
-
+
@@ -1073,7 +1085,7 @@
-
+
@@ -1086,7 +1098,7 @@
-
+
@@ -1101,7 +1113,7 @@
-
+
@@ -1116,7 +1128,7 @@
-
+
@@ -1129,7 +1141,7 @@
-
+
@@ -1144,7 +1156,7 @@
-
+
@@ -1155,7 +1167,7 @@
-
+
@@ -1168,7 +1180,7 @@
-
+
@@ -1181,7 +1193,7 @@
-
+
@@ -1196,7 +1208,7 @@
-
+
@@ -1212,7 +1224,7 @@
-
+
@@ -1227,7 +1239,7 @@
-
+
@@ -1244,7 +1256,7 @@
-
+
@@ -1261,7 +1273,7 @@
-
+
@@ -1278,7 +1290,7 @@
-
+
@@ -1293,7 +1305,7 @@
-
+
@@ -1306,7 +1318,7 @@
-
+
@@ -1345,7 +1357,7 @@
-
+
@@ -1374,7 +1386,7 @@
-
+
@@ -1467,7 +1479,7 @@
-
+
@@ -3091,7 +3103,7 @@
-
+
@@ -3119,7 +3131,7 @@
-
+
@@ -3138,11 +3150,11 @@
-
+
-
+
@@ -3152,7 +3164,7 @@
-
+
@@ -3162,11 +3174,11 @@
-
+
-
+
@@ -3176,7 +3188,7 @@
-
+
@@ -3186,7 +3198,7 @@
-
+
@@ -3196,7 +3208,7 @@
-
+
@@ -3220,7 +3232,7 @@
-
+
@@ -3238,7 +3250,7 @@
-
+
@@ -3250,7 +3262,7 @@
-
+
@@ -3262,7 +3274,7 @@
-
+
@@ -3270,11 +3282,11 @@
-
+
-
+
@@ -3297,7 +3309,7 @@
-
+
@@ -3323,7 +3335,7 @@
-
+
@@ -3344,7 +3356,7 @@
-
+
@@ -3374,7 +3386,7 @@
-
+
@@ -3388,7 +3400,7 @@
-
+
@@ -3415,7 +3427,7 @@
-
+
@@ -3454,7 +3466,7 @@
-
+