[cig-commits] r18346 - short/3D/PyLith/branches/v1.5-stable/libsrc/friction

[brad at geodynamics.org]

Author: brad
Date: 2011-05-11 16:45:49 -0700 (Wed, 11 May 2011)
New Revision: 18346

Modified:
   short/3D/PyLith/branches/v1.5-stable/libsrc/friction/RateStateAgeing.cc
Log:
Clean up of state variable update for rate and state friction.

Modified: short/3D/PyLith/branches/v1.5-stable/libsrc/friction/RateStateAgeing.cc
===================================================================
--- short/3D/PyLith/branches/v1.5-stable/libsrc/friction/RateStateAgeing.cc	2011-05-11 23:45:17 UTC (rev 18345)
+++ short/3D/PyLith/branches/v1.5-stable/libsrc/friction/RateStateAgeing.cc	2011-05-11 23:45:49 UTC (rev 18346)
@@ -335,31 +335,42 @@
   assert(0 != numStateVars);
   assert(0 != numProperties);
 
+  // d(theta)/dt = (1 - slipRate * theta / L)
+  //
+  // Use separation of variables to integrate above ODE from t->t+dt,
+  // keeping slip rate constant.
+  //
+  // thetaTpdt = thetaT * exp(-slipRate/L * dt)
+  //             + L/slipRate * (1 -  exp(-slipRate/L * dt))
+  //
+  // As slipRate --> 0, L/sliprate --> infinity and
+  // exp(-sliprate/L*dt) --> 1.  To determine, d(theta)/dt near
+  // sliprate == 0, we expand the exponential term in a Taylor series:
+  //
+  // exp(-x) = 1 - x +1/2*x**2 + 1/6*x**3
+  //
+  // This leads to (in the vicinity of slipRate == 0):
+  //
+  // thetaTpdt = thetaT * exp(-slipRate/L * dt)
+  //             + dt - 0.5*(sliprate/L)*dt**2 + 1.0/6.0*(slipRate/L)*dt**3;
+
   const double dt = _dt;
   const double thetaTVertex = stateVars[s_state];
   const double L = properties[p_L];
-  const double expTerm = exp(-slipRate * dt / L);
   const double vDtL = slipRate * dt / L;
+  const double expTerm = exp(-vDtL);
 
-  // Ageing law
-  // d(theta)/dt = (1 - slipRate * theta / L)
-  // Above ODE is integrated from t->t+dt keeping slipRate constant gives
-  // thetaTpdt = thetaT * exp(- slipRate * theta / L)
-  //             + L / slipRate * (1 -  exp(- slipRate * theta / L))
   double thetaTpdtVertex = 0.0;
-  if (vDtL < 0.00001)
-    // As (slipRate * dt / L) --> 0, exp(-slipRate * dt / L) --> 1
-    // So using first three terms in the Taylor series expansion of 
-    // exp(- slipRate * theta / L) i.e., exp(-x) = 1 - x + (x^2)/2;
-    thetaTpdtVertex = thetaTVertex * expTerm + 
-                        dt - 0.5 * slipRate * pow(dt,2) / L;
-  else
-    thetaTpdtVertex = thetaTVertex * expTerm +
-                       L / slipRate * (1 - expTerm);
-
+  if (vDtL > 1.0e-20) {
+    thetaTpdtVertex = thetaTVertex * expTerm + L / slipRate * (1 - expTerm);
+    PetscLogFlops(7);
+  } else {
+    thetaTpdtVertex = thetaTVertex * expTerm + dt - 0.5 * slipRate/L * dt*dt;
+    PetscLogFlops(9);
+  } // if/else
+  
   stateVars[s_state] = thetaTpdtVertex;
 
-  PetscLogFlops(6);
 } // _updateStateVars