GEOS
RateAndStateKernels.hpp
1 /*
2  * ------------------------------------------------------------------------------------------------------------
3  * SPDX-License-Identifier: LGPL-2.1-only
4  *
5  * Copyright (c) 2016-2024 Lawrence Livermore National Security LLC
6  * Copyright (c) 2018-2024 TotalEnergies
7  * Copyright (c) 2018-2024 The Board of Trustees of the Leland Stanford Junior University
8  * Copyright (c) 2018-2024 Chevron
9  * Copyright (c) 2019- GEOS/GEOSX Contributors
10  * All rights reserved
11  *
12  * See top level LICENSE, COPYRIGHT, CONTRIBUTORS, NOTICE, and ACKNOWLEDGEMENTS files for details.
13  * ------------------------------------------------------------------------------------------------------------
14  */
15 
16 #ifndef GEOS_PHYSICSSOLVERS_RATEANDSTATEKERNELS_HPP_
17 #define GEOS_PHYSICSSOLVERS_RATEANDSTATEKERNELS_HPP_
18 
19 #include "common/DataTypes.hpp"
20 #include "common/GEOS_RAJA_Interface.hpp"
21 #include "constitutive/contact/RateAndStateFriction.hpp"
22 #include "physicsSolvers/inducedSeismicity/rateAndStateFields.hpp"
23 #include "denseLinearAlgebra/denseLASolvers.hpp"
24 
25 namespace geos
26 {
27 
28 namespace rateAndStateKernels
29 {
30 
31 // TBD: Pass the kernel and add getters for relevant fields to make this function general purpose and avoid
32 // wrappers?
33 GEOS_HOST_DEVICE
34 static void projectSlipRateBase( localIndex const k,
35  real64 const frictionCoefficient,
36  real64 const shearImpedance,
37  arrayView2d< real64 const > const traction,
38  arrayView1d< real64 const > const slipRate,
39  arrayView2d< real64 > const slipVelocity )
40 {
41  // Project slip rate onto shear traction to get slip velocity components
42  real64 const frictionForce = traction[k][0] * frictionCoefficient;
43  real64 const projectionScaling = 1.0 / ( shearImpedance + frictionForce / slipRate[k] );
44  slipVelocity[k][0] = projectionScaling * traction[k][1];
45  slipVelocity[k][1] = projectionScaling * traction[k][2];
46 }
47 
55 class ImplicitFixedStressRateAndStateKernel
56 {
57 public:
58 
59  ImplicitFixedStressRateAndStateKernel( SurfaceElementSubRegion & subRegion,
60  constitutive::RateAndStateFriction const & frictionLaw,
61  real64 const shearImpedance ):
62  m_slipRate( subRegion.getField< fields::rateAndState::slipRate >() ),
63  m_stateVariable( subRegion.getField< fields::rateAndState::stateVariable >() ),
64  m_stateVariable_n( subRegion.getField< fields::rateAndState::stateVariable_n >() ),
65  m_traction( subRegion.getField< fields::contact::traction >() ),
66  m_slipVelocity( subRegion.getField< fields::rateAndState::slipVelocity >() ),
67  m_shearImpedance( shearImpedance ),
68  m_frictionLaw( frictionLaw.createKernelUpdates() )
69  {}
70 
75  struct StackVariables
76  {
77 public:
78 
79  GEOS_HOST_DEVICE
80  StackVariables( )
81  {}
82 
83  real64 jacobian[2][2]{};
84 
85  real64 rhs[2]{};
86 
87  };
88 
89  GEOS_HOST_DEVICE
90  void setup( localIndex const k,
91  real64 const dt,
92  StackVariables & stack ) const
93  {
94  real64 const normalTraction = m_traction[k][0];
95  real64 const shearTractionMagnitude = LvArray::math::sqrt( m_traction[k][1] * m_traction[k][1] + m_traction[k][2] * m_traction[k][2] );
96  // Eq 1: Scalar force balance for slipRate and shear traction magnitude
97  stack.rhs[0] = shearTractionMagnitude - m_shearImpedance * m_slipRate[k]
98  - normalTraction * m_frictionLaw.frictionCoefficient( k, m_slipRate[k], m_stateVariable[k] );
99  real64 const dFriction[2] = { -normalTraction * m_frictionLaw.dFrictionCoefficient_dStateVariable( k, m_slipRate[k], m_stateVariable[k] ),
100  -m_shearImpedance - normalTraction * m_frictionLaw.dFrictionCoefficient_dSlipRate( k, m_slipRate[k], m_stateVariable[k] ) };
101 
102  // Eq 2: slip law
103  stack.rhs[1] = (m_stateVariable[k] - m_stateVariable_n[k]) / dt - m_frictionLaw.stateEvolution( k, m_slipRate[k], m_stateVariable[k] );
104  real64 const dStateEvolutionLaw[2] = { 1 / dt - m_frictionLaw.dStateEvolution_dStateVariable( k, m_slipRate[k], m_stateVariable[k] ),
105  -m_frictionLaw.dStateEvolution_dSlipRate( k, m_slipRate[k], m_stateVariable[k] ) };
106 
107  // Assemble Jacobian matrix
108  stack.jacobian[0][0] = dFriction[0]; // derivative of Eq 1 w.r.t. stateVariable
109  stack.jacobian[0][1] = dFriction[1]; // derivative of Eq 1 w.r.t. slipRate
110  stack.jacobian[1][0] = dStateEvolutionLaw[0]; // derivative of Eq 2 w.r.t. stateVariable
111  stack.jacobian[1][1] = dStateEvolutionLaw[1]; // derivative of Eq 2 w.r.t. slipRate
112  }
113 
114  GEOS_HOST_DEVICE
115  void solve( localIndex const k,
116  StackVariables & stack ) const
117  {
119  real64 solution[2] = {0.0, 0.0};
120  denseLinearAlgebra::solve< 2 >( stack.jacobian, stack.rhs, solution );
121 
122  // Update variables
123  m_stateVariable[k] -= solution[0];
124  m_slipRate[k] -= solution[1];
125  }
126 
127  GEOS_HOST_DEVICE
128  void projectSlipRate( localIndex const k ) const
129  {
130  real64 const frictionCoefficient = m_frictionLaw.frictionCoefficient( k, m_slipRate[k], m_stateVariable[k] );
131  projectSlipRateBase( k, frictionCoefficient, m_shearImpedance, m_traction, m_slipRate, m_slipVelocity );
132  }
133 
134  GEOS_HOST_DEVICE
135  camp::tuple< int, real64 > checkConvergence( StackVariables const & stack,
136  real64 const tol ) const
137  {
138  real64 const residualNorm = LvArray::tensorOps::l2Norm< 2 >( stack.rhs );
139  int const converged = residualNorm < tol ? 1 : 0;
140  camp::tuple< int, real64 > result { converged, residualNorm };
141  return result;
142  }
143 
144 private:
145 
146  arrayView1d< real64 > const m_slipRate;
147 
148  arrayView1d< real64 > const m_stateVariable;
149 
150  arrayView1d< real64 const > const m_stateVariable_n;
151 
152  arrayView2d< real64 const > const m_traction;
153 
154  arrayView2d< real64 > const m_slipVelocity;
155 
156  real64 const m_shearImpedance;
157 
158  constitutive::RateAndStateFriction::KernelWrapper m_frictionLaw;
159 
160 };
161 
169 class ExplicitRateAndStateKernel
170 {
171 public:
172 
173  ExplicitRateAndStateKernel( SurfaceElementSubRegion & subRegion,
174  constitutive::RateAndStateFriction const & frictionLaw,
175  real64 const shearImpedance ):
176  m_slipRate( subRegion.getField< fields::rateAndState::slipRate >() ),
177  m_stateVariable( subRegion.getField< fields::rateAndState::stateVariable >() ),
178  m_traction( subRegion.getField< fields::contact::traction >() ),
179  m_slipVelocity( subRegion.getField< fields::rateAndState::slipVelocity >() ),
180  m_shearImpedance( shearImpedance ),
181  m_frictionLaw( frictionLaw.createKernelUpdates() )
182  {}
183 
188  struct StackVariables
189  {
190 public:
191 
192  GEOS_HOST_DEVICE
193  StackVariables( )
194  {}
195 
196  real64 jacobian;
197  real64 rhs;
198 
199  };
200 
201  GEOS_HOST_DEVICE
202  void setup( localIndex const k,
203  real64 const dt,
204  StackVariables & stack ) const
205  {
206  GEOS_UNUSED_VAR( dt );
207  real64 const normalTraction = m_traction[k][0];
208  real64 const shearTractionMagnitude = LvArray::math::sqrt( m_traction[k][1] * m_traction[k][1] + m_traction[k][2] * m_traction[k][2] );
209 
210  // Slip rate is bracketed between [0, shear traction magnitude / shear impedance]
211  // If slip rate is outside the bracket, re-initialize to the middle value
212  real64 const upperBound = shearTractionMagnitude/m_shearImpedance;
213  real64 const bracketedSlipRate = m_slipRate[k] > upperBound ? 0.5*upperBound : m_slipRate[k];
214 
215  stack.rhs = shearTractionMagnitude - m_shearImpedance *bracketedSlipRate - normalTraction * m_frictionLaw.frictionCoefficient( k, bracketedSlipRate, m_stateVariable[k] );
216  stack.jacobian = -m_shearImpedance - normalTraction * m_frictionLaw.dFrictionCoefficient_dSlipRate( k, bracketedSlipRate, m_stateVariable[k] );
217  }
218 
219  GEOS_HOST_DEVICE
220  void solve( localIndex const k,
221  StackVariables & stack ) const
222  {
223  m_slipRate[k] -= stack.rhs/stack.jacobian;
224 
225  // Slip rate is bracketed between [0, shear traction magnitude / shear impedance]
226  // Check that the update did not end outside of the bracket.
227  real64 const shearTractionMagnitude = LvArray::math::sqrt( m_traction[k][1] * m_traction[k][1] + m_traction[k][2] * m_traction[k][2] );
228  real64 const upperBound = shearTractionMagnitude/m_shearImpedance;
229  if( m_slipRate[k] > upperBound ) m_slipRate[k] = 0.5*upperBound;
230 
231  }
232 
233 
234  GEOS_HOST_DEVICE
235  camp::tuple< int, real64 > checkConvergence( StackVariables const & stack,
236  real64 const tol ) const
237  {
238  real64 const residualNorm = LvArray::math::abs( stack.rhs );
239  int const converged = residualNorm < tol ? 1 : 0;
240  camp::tuple< int, real64 > result { converged, residualNorm };
241  return result;
242  }
243 
244  GEOS_HOST_DEVICE
245  void projectSlipRate( localIndex const k ) const
246  {
247  real64 const frictionCoefficient = m_frictionLaw.frictionCoefficient( k, m_slipRate[k], m_stateVariable[k] );
248  projectSlipRateBase( k, frictionCoefficient, m_shearImpedance, m_traction, m_slipRate, m_slipVelocity );
249  }
250 
251 private:
252 
253  arrayView1d< real64 > const m_slipRate;
254 
255  arrayView1d< real64 > const m_stateVariable;
256 
257  arrayView2d< real64 const > const m_traction;
258 
259  arrayView2d< real64 > const m_slipVelocity;
260 
261  real64 const m_shearImpedance;
262 
263  constitutive::RateAndStateFriction::KernelWrapper m_frictionLaw;
264 
265 };
266 
267 
273 template< typename KernelType, typename POLICY >
274 static void
275 createAndLaunch( SurfaceElementSubRegion & subRegion,
276  string const & frictionLawNameKey,
277  real64 const shearImpedance,
278  integer const maxIterNewton,
279  real64 const newtonTol,
280  real64 const time_n,
281  real64 const dt )
282 {
283  GEOS_MARK_FUNCTION;
284 
285  GEOS_UNUSED_VAR( time_n );
286 
287  string const & frictionaLawName = subRegion.getReference< string >( frictionLawNameKey );
288  constitutive::RateAndStateFriction const & frictionLaw = subRegion.getConstitutiveModel< constitutive::RateAndStateFriction >( frictionaLawName );
289  KernelType kernel( subRegion, frictionLaw, shearImpedance );
290 
291  // Newton loop (outside of the kernel launch)
292  bool allConverged = false;
293  for( integer iter = 0; iter < maxIterNewton; iter++ )
294  {
295  RAJA::ReduceMin< parallelDeviceReduce, int > converged( 1 );
296  RAJA::ReduceMax< parallelDeviceReduce, real64 > residualNorm( 0.0 );
297  forAll< POLICY >( subRegion.size(), [=] GEOS_HOST_DEVICE ( localIndex const k )
298  {
299  typename KernelType::StackVariables stack;
300  kernel.setup( k, dt, stack );
301  kernel.solve( k, stack );
302  auto const [elementConverged, elementResidualNorm] = kernel.checkConvergence( stack, newtonTol );
303  converged.min( elementConverged );
304  residualNorm.max( elementResidualNorm );
305  } );
306 
307  real64 const maxResidualNorm = MpiWrapper::max( residualNorm.get() );
308  GEOS_LOG_RANK_0( GEOS_FMT( "-----iter {} : residual = {:.10e} ", iter, maxResidualNorm ) );
309 
310  if( converged.get() )
311  {
312  allConverged = true;
313  break;
314  }
315  }
316  if( !allConverged )
317  {
318  GEOS_ERROR( " Failed to converge" );
319  }
320  forAll< POLICY >( subRegion.size(), [=] GEOS_HOST_DEVICE ( localIndex const k )
321  {
322  kernel.projectSlipRate( k );
323  } );
324 }
325 
331 struct Kutta32Table
332 {
333  integer constexpr static algHighOrder = 3; // High-order update order
334  integer constexpr static algLowOrder = 2; // Low-order update order
335  integer constexpr static numStages = 3; // Number of stages
336  real64 const a[2][2] = { { 1.0/2.0, 0.0 }, // Coefficients for stage value updates
337  { -1.0, 2.0 } }; // (lower-triangular part of table).
338  real64 const c[3] = { 0.0, 1.0/2.0, 1.0 }; // Coefficients for time increments of substages
339  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
340  real64 const bStar[3] = { 1.0/2.0, 0.0, 1.0/2.0 }; // Quadrature weights used for low-order comparision solution
341  real64 constexpr static FSAL = false; // Not first same as last
342 };
343 
347 struct BogackiShampine32Table
348 {
349  integer constexpr static algHighOrder = 3; // High-order update order
350  integer constexpr static algLowOrder = 2; // Low-order update order
351  integer constexpr static numStages = 4; // Number of stages
352  real64 const a[3][3] = { { 1.0/2.0, 0.0, 0.0 }, // Coefficients for stage value updates
353  { 0.0, 3.0/4.0, 0.0 }, // (lower-triangular part of table).
354  { 2.0/9.0, 1.0/3.0, 4.0/9.0 } };
355  real64 const c[4] = { 0.0, 1.0/2.0, 3.0/4.0, 1.0 }; // Coefficients for time increments of substages
356  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
357  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
358  bool constexpr static FSAL = true; // First same as last (can reuse the last stage rate in next
359  // update)
360 };
361 
369 template< typename TABLE_TYPE > class EmbeddedRungeKuttaKernel
370 {
371 
372 public:
373  EmbeddedRungeKuttaKernel( SurfaceElementSubRegion & subRegion,
374  constitutive::RateAndStateFriction const & frictionLaw,
375  TABLE_TYPE butcherTable ):
376  m_stateVariable( subRegion.getField< fields::rateAndState::stateVariable >() ),
377  m_stateVariable_n( subRegion.getField< fields::rateAndState::stateVariable_n >() ),
378  m_slipRate( subRegion.getField< fields::rateAndState::slipRate >() ),
379  m_slipVelocity( subRegion.getField< fields::rateAndState::slipVelocity >() ),
380  m_slipVelocity_n( subRegion.getField< fields::rateAndState::slipVelocity_n >() ),
381  m_deltaSlip( subRegion.getField< fields::rateAndState::deltaSlip >() ),
382  m_deltaSlip_n( subRegion.getField< fields::rateAndState::deltaSlip_n >() ),
383  m_dispJump( subRegion.getField< fields::contact::dispJump >() ),
384  m_dispJump_n( subRegion.getField< fields::contact::dispJump_n >() ),
385  m_error( subRegion.getField< fields::rateAndState::error >() ),
386  m_stageRates( subRegion.getField< fields::rateAndState::rungeKuttaStageRates >() ),
387  m_frictionLaw( frictionLaw.createKernelUpdates() ),
388  m_butcherTable( butcherTable )
389  {}
390 
394  GEOS_HOST_DEVICE
395  void initialize( localIndex const k ) const
396  {
397  LvArray::tensorOps::copy< 2 >( m_slipVelocity[k], m_slipVelocity_n[k] );
398  m_slipRate[k] = LvArray::tensorOps::l2Norm< 2 >( m_slipVelocity_n[k] );
399  m_stateVariable[k] = m_stateVariable_n[k];
400  }
401 
406  GEOS_HOST_DEVICE
407  void updateStageRatesFSAL( localIndex const k ) const
408  {
409  LvArray::tensorOps::copy< 3 >( m_stageRates[k][0], m_stageRates[k][m_butcherTable.numStages-1] );
410  }
411 
415  GEOS_HOST_DEVICE
416  void updateStageRates( localIndex const k, integer const stageIndex ) const
417  {
418  m_stageRates[k][stageIndex][0] = m_slipVelocity[k][0];
419  m_stageRates[k][stageIndex][1] = m_slipVelocity[k][1];
420  m_stageRates[k][stageIndex][2] = m_frictionLaw.stateEvolution( k, m_slipRate[k], m_stateVariable[k] );
421  }
422 
426  GEOS_HOST_DEVICE
427  void updateStageValues( localIndex const k, integer const stageIndex, real64 const dt ) const
428  {
429 
430  real64 stateVariableIncrement = 0.0;
431  real64 deltaSlipIncrement[2] = {0.0, 0.0};
432 
433  for( integer i = 0; i < stageIndex; i++ )
434  {
435  deltaSlipIncrement[0] += m_butcherTable.a[stageIndex-1][i] * m_stageRates[k][i][0];
436  deltaSlipIncrement[1] += m_butcherTable.a[stageIndex-1][i] * m_stageRates[k][i][1];
437  stateVariableIncrement += m_butcherTable.a[stageIndex-1][i] * m_stageRates[k][i][2];
438  }
439  m_deltaSlip[k][0] = m_deltaSlip_n[k][0] + dt*deltaSlipIncrement[0];
440  m_deltaSlip[k][1] = m_deltaSlip_n[k][1] + dt*deltaSlipIncrement[1];
441  m_stateVariable[k] = m_stateVariable_n[k] + dt*stateVariableIncrement;
442 
443  m_dispJump[k][1] = m_dispJump_n[k][1] + m_deltaSlip[k][0];
444  m_dispJump[k][2] = m_dispJump_n[k][2] + m_deltaSlip[k][1];
445  }
446 
451  GEOS_HOST_DEVICE
452  void updateSolutionAndLocalError( localIndex const k, real64 const dt, real64 const absTol, real64 const relTol ) const
453  {
454 
455  real64 deltaSlipIncrement[2] = {0.0, 0.0};
456  real64 deltaSlipIncrementLowOrder[2] = {0.0, 0.0};
457 
458  real64 stateVariableIncrement = 0.0;
459  real64 stateVariableIncrementLowOrder = 0.0;
460 
461  for( localIndex i = 0; i < m_butcherTable.numStages; i++ )
462  {
463 
464  // High order update of solution
465  deltaSlipIncrement[0] += m_butcherTable.b[i] * m_stageRates[k][i][0];
466  deltaSlipIncrement[1] += m_butcherTable.b[i] * m_stageRates[k][i][1];
467  stateVariableIncrement += m_butcherTable.b[i] * m_stageRates[k][i][2];
468 
469  // Low order update for error
470  deltaSlipIncrementLowOrder[0] += m_butcherTable.bStar[i] * m_stageRates[k][i][0];
471  deltaSlipIncrementLowOrder[1] += m_butcherTable.bStar[i] * m_stageRates[k][i][1];
472  stateVariableIncrementLowOrder += m_butcherTable.bStar[i] * m_stageRates[k][i][2];
473  }
474 
475  m_deltaSlip[k][0] = m_deltaSlip_n[k][0] + dt * deltaSlipIncrement[0];
476  m_deltaSlip[k][1] = m_deltaSlip_n[k][1] + dt * deltaSlipIncrement[1];
477  m_stateVariable[k] = m_stateVariable_n[k] + dt * stateVariableIncrement;
478 
479  real64 const deltaSlipLowOrder[2] = {m_deltaSlip_n[k][0] + dt * deltaSlipIncrementLowOrder[0],
480  m_deltaSlip_n[k][1] + dt * deltaSlipIncrementLowOrder[1]};
481  real64 const stateVariableLowOrder = m_stateVariable_n[k] + dt * stateVariableIncrementLowOrder;
482 
483  m_dispJump[k][1] = m_dispJump_n[k][1] + m_deltaSlip[k][0];
484  m_dispJump[k][2] = m_dispJump_n[k][2] + m_deltaSlip[k][1];
485 
486  // Compute error
487  m_error[k][0] = computeError( m_deltaSlip[k][0], deltaSlipLowOrder[0], absTol, relTol );
488  m_error[k][1] = computeError( m_deltaSlip[k][1], deltaSlipLowOrder[1], absTol, relTol );
489  m_error[k][2] = computeError( m_stateVariable[k], stateVariableLowOrder, absTol, relTol );
490  }
491 
496  GEOS_HOST_DEVICE
497  void updateSolutionAndLocalErrorFSAL( localIndex const k, real64 const dt, real64 const absTol, real64 const relTol ) const
498  {
499 
500  real64 deltaSlipIncrementLowOrder[2] = {0.0, 0.0};
501  real64 stateVariableIncrementLowOrder = 0.0;
502 
503  for( localIndex i = 0; i < m_butcherTable.numStages; i++ )
504  {
505  // In FSAL algorithms the last RK substage update coincides with the
506  // high-order update. Only need to compute increments for the the
507  // low-order updates for error computation.
508  deltaSlipIncrementLowOrder[0] += m_butcherTable.bStar[i] * m_stageRates[k][i][0];
509  deltaSlipIncrementLowOrder[1] += m_butcherTable.bStar[i] * m_stageRates[k][i][1];
510  stateVariableIncrementLowOrder += m_butcherTable.bStar[i] * m_stageRates[k][i][2];
511  }
512 
513  real64 const deltaSlipLowOrder[2] = {m_deltaSlip_n[k][0] + dt * deltaSlipIncrementLowOrder[0],
514  m_deltaSlip_n[k][1] + dt * deltaSlipIncrementLowOrder[1]};
515  real64 const stateVariableLowOrder = m_stateVariable_n[k] + dt * stateVariableIncrementLowOrder;
516 
517  m_dispJump[k][1] = m_dispJump_n[k][1] + m_deltaSlip[k][0];
518  m_dispJump[k][2] = m_dispJump_n[k][2] + m_deltaSlip[k][1];
519 
520  // Compute error
521  m_error[k][0] = computeError( m_deltaSlip[k][0], deltaSlipLowOrder[0], absTol, relTol );
522  m_error[k][1] = computeError( m_deltaSlip[k][1], deltaSlipLowOrder[1], absTol, relTol );
523  m_error[k][2] = computeError( m_stateVariable[k], stateVariableLowOrder, absTol, relTol );
524  }
525 
529  GEOS_HOST_DEVICE
530  real64 computeError( real64 const highOrderApprox, real64 const lowOrderApprox, real64 const absTol, real64 const relTol ) const
531  {
532  return (highOrderApprox - lowOrderApprox) /
533  ( absTol + relTol * LvArray::math::max( LvArray::math::abs( highOrderApprox ), LvArray::math::abs( lowOrderApprox ) ));
534  }
535 
536 private:
537 
539  arrayView1d< real64 > const m_stateVariable;
540 
542  arrayView1d< real64 > const m_stateVariable_n;
543 
545  arrayView1d< real64 > const m_slipRate;
546 
548  arrayView2d< real64 > const m_slipVelocity;
549 
551  arrayView2d< real64 > const m_slipVelocity_n;
552 
554  arrayView2d< real64 > const m_deltaSlip;
555 
557  arrayView2d< real64 > const m_deltaSlip_n;
558 
560  arrayView2d< real64 > const m_dispJump;
561 
563  arrayView2d< real64 > const m_dispJump_n;
564 
566  arrayView2d< real64 > const m_error;
567 
569  arrayView3d< real64 > const m_stageRates;
570 
572  constitutive::RateAndStateFriction::KernelWrapper m_frictionLaw;
573 
575  TABLE_TYPE m_butcherTable;
576 };
577 
578 } /* namespace rateAndStateKernels */
579 
580 } /* namespace geos */
581 
582 #endif /* GEOS_PHYSICSSOLVERS_RATEANDSTATEKERNELS_HPP_ */
GEOS_HOST_DEVICE
#define GEOS_HOST_DEVICE
Marks a host-device function.
Definition: GeosxMacros.hpp:49
GEOS_UNUSED_VAR
#define GEOS_UNUSED_VAR(...)
Mark an unused variable and silence compiler warnings.
Definition: GeosxMacros.hpp:84
GEOS_ERROR
#define GEOS_ERROR(msg)
Raise a hard error and terminate the program.
Definition: Logger.hpp:157
GEOS_LOG_RANK_0
#define GEOS_LOG_RANK_0(msg)
Log a message on screen on rank 0.
Definition: Logger.hpp:101
GEOS_MARK_FUNCTION
#define GEOS_MARK_FUNCTION
Mark function with both Caliper and NVTX if enabled.
Definition: TimingMacros.hpp:134
geos::ObjectManagerBase::getField
GEOS_DECLTYPE_AUTO_RETURN getField() const
Get a view to the field associated with a trait from this ObjectManagerBase.
Definition: ObjectManagerBase.hpp:627
geos::rateAndStateKernels::EmbeddedRungeKuttaKernel
Runge-Kutta method used to time integrate slip and state. Uses of a high order update used to integra...
Definition: RateAndStateKernels.hpp:370
geos::rateAndStateKernels::EmbeddedRungeKuttaKernel::computeError
GEOS_HOST_DEVICE real64 computeError(real64 const highOrderApprox, real64 const lowOrderApprox, real64 const absTol, real64 const relTol) const
Computes the relative error scaled by error tolerances.
Definition: RateAndStateKernels.hpp:530
geos::rateAndStateKernels::EmbeddedRungeKuttaKernel::updateSolutionAndLocalError
GEOS_HOST_DEVICE void updateSolutionAndLocalError(localIndex const k, real64 const dt, real64 const absTol, real64 const relTol) const
Updates slip, state and displacement jump to the next time computes error the local error in the time...
Definition: RateAndStateKernels.hpp:452
geos::rateAndStateKernels::EmbeddedRungeKuttaKernel::updateStageRates
GEOS_HOST_DEVICE void updateStageRates(localIndex const k, integer const stageIndex) const
Updates the stage rates rates (the right-hand-side of the ODEs for slip and state)
Definition: RateAndStateKernels.hpp:416
geos::rateAndStateKernels::EmbeddedRungeKuttaKernel::updateStageRatesFSAL
GEOS_HOST_DEVICE void updateStageRatesFSAL(localIndex const k) const
Re-uses the last stage rate from the previous time step as the first in the next update....
Definition: RateAndStateKernels.hpp:407
geos::rateAndStateKernels::EmbeddedRungeKuttaKernel::updateStageValues
GEOS_HOST_DEVICE void updateStageValues(localIndex const k, integer const stageIndex, real64 const dt) const
Update stage values (slip, state and displacement jump) to a Runge-Kutta substage.
Definition: RateAndStateKernels.hpp:427
geos::rateAndStateKernels::EmbeddedRungeKuttaKernel::initialize
GEOS_HOST_DEVICE void initialize(localIndex const k) const
Initialize slip and state buffers.
Definition: RateAndStateKernels.hpp:395
geos::rateAndStateKernels::EmbeddedRungeKuttaKernel::updateSolutionAndLocalErrorFSAL
GEOS_HOST_DEVICE void updateSolutionAndLocalErrorFSAL(localIndex const k, real64 const dt, real64 const absTol, real64 const relTol) const
Updates slip, state and displacement jump to the next time computes error the local error in the time...
Definition: RateAndStateKernels.hpp:497
geos::rateAndStateKernels::ImplicitFixedStressRateAndStateKernel::solve
GEOS_HOST_DEVICE void solve(localIndex const k, StackVariables &stack) const
Definition: RateAndStateKernels.hpp:115
geos::arrayView1d
ArrayView< T, 1 > arrayView1d
Alias for 1D array view.
Definition: DataTypes.hpp:180
geos::real64
double real64
64-bit floating point type.
Definition: DataTypes.hpp:99
geos::localIndex
GEOS_LOCALINDEX_TYPE localIndex
Local index type (for indexing objects within an MPI partition).
Definition: DataTypes.hpp:85
geos::integer
std::int32_t integer
Signed integer type.
Definition: DataTypes.hpp:82
geos::arrayView2d
ArrayView< T, 2, USD > arrayView2d
Alias for 2D array view.
Definition: DataTypes.hpp:196
geos::arrayView3d
ArrayView< T, 3, USD > arrayView3d
Alias for 3D array view.
Definition: DataTypes.hpp:212
geos::MpiWrapper::max
static T max(T const &value, MPI_Comm comm=MPI_COMM_GEOS)
Convenience function for a MPI_Allreduce using a MPI_MAX operation.
Definition: MpiWrapper.hpp:1101
geos::rateAndStateKernels::BogackiShampine32Table
Butcher table for the BogackiShampine 3(2) method.
Definition: RateAndStateKernels.hpp:348
geos::rateAndStateKernels::Kutta32Table
Butcher table for embedded RK3(2) method using Kuttas third order method for the high-order update,...
Definition: RateAndStateKernels.hpp:332
Generated by doxygen 1.9.1