[cig-commits] r18467 - in short/3D/PyLith/trunk: . libsrc/pylith/faults libsrc/pylith/friction
brad at geodynamics.org
brad at geodynamics.org
Thu May 26 11:27:00 PDT 2011
Author: brad
Date: 2011-05-26 11:27:00 -0700 (Thu, 26 May 2011)
New Revision: 18467
Modified:
short/3D/PyLith/trunk/TODO
short/3D/PyLith/trunk/libsrc/pylith/faults/FaultCohesiveDyn.cc
short/3D/PyLith/trunk/libsrc/pylith/friction/RateStateAgeing.cc
Log:
Fixed bug in fault friction implementation. Need to account for possibility of overshooting slip estimate.
Modified: short/3D/PyLith/trunk/TODO
===================================================================
--- short/3D/PyLith/trunk/TODO 2011-05-26 16:13:29 UTC (rev 18466)
+++ short/3D/PyLith/trunk/TODO 2011-05-26 18:27:00 UTC (rev 18467)
@@ -4,8 +4,7 @@
* BUILDBOT
- (1) cygwin
- + packaging
+ cygwin packaging
----------------------------------------
BUGS
@@ -15,9 +14,6 @@
rate and state friction - sliding test
predictor/corrector scheme?
-* Manual
- top of pg 86 , k -> p
-
* Unform global refinement
Check refinement in parallel thoroughly.
@@ -36,11 +32,11 @@
MISSING FEATURES
----------------------------------------
-* Output to HDF5 files. [BRAD and MATT]
+* Output to HDF5 files. [BRAD]
- (2) Add time dataset to vertex_fields and cell_fields. [BRAD]
+ (1) Add time dataset to vertex_fields and cell_fields.
- (3) Add creation of .xmf metadata file for ParaView/Visit. [BRAD]
+ (2) Add creation of .xmf metadata file for ParaView/Visit.
* Field split.
@@ -50,7 +46,7 @@
Need to check performance of custom fault preconditioner.
-* Cleanup of materials [stable]
+* Cleanup of materials
+ Account for initial strain in viscoelastic models
@@ -65,11 +61,11 @@
FaultCohesiveLagrange::calcPreconditioner() [need unit test]
-* Cleanup
- + memory model
- + full-scale testing
- + Code cleanup
+----------------------------------------
+MISCELLANEOUS
+----------------------------------------
+
* Memory model [MATT and BRAD]
(1) The memory logging for distribution logs Creation and Stratification
Modified: short/3D/PyLith/trunk/libsrc/pylith/faults/FaultCohesiveDyn.cc
===================================================================
--- short/3D/PyLith/trunk/libsrc/pylith/faults/FaultCohesiveDyn.cc 2011-05-26 16:13:29 UTC (rev 18466)
+++ short/3D/PyLith/trunk/libsrc/pylith/faults/FaultCohesiveDyn.cc 2011-05-26 18:27:00 UTC (rev 18467)
@@ -435,6 +435,7 @@
#if 0 // DEBUGGING
slipSection->view("SLIP");
+ slipRateSection->view("SLIP RATE");
areaSection->view("AREA");
dispTSection->view("DISP (t)");
dispTIncrSection->view("DISP INCR (t->t+dt)");
@@ -584,6 +585,7 @@
dLagrangeTpdtSection->view("AFTER dLagrange");
dispTIncrSection->view("AFTER DISP INCR (t->t+dt)");
slipSection->view("AFTER SLIP");
+ slipRateSection->view("AFTER SLIP RATE");
#endif
} // constrainSolnSpace
@@ -1943,17 +1945,20 @@
// if in compression and no opening
const double frictionStress = _friction->calcFriction(slipMag, slipRateMag,
tractionNormal);
- if (tractionShearMag > frictionStress) {
- // traction is limited by friction, so have sliding
+ if (tractionShearMag > frictionStress || slipRateMag > 0.0) {
+ // traction is limited by friction, so have sliding OR
+ // friction exceeds traction due to overshoot in slip
- // Update slip based on value required to stick versus friction
+ // Update traction increment based on value required to stick
+ // versus friction
const double dlp = -(tractionShearMag - frictionStress) * area *
tractionTpdt[0] / tractionShearMag;
(*dLagrangeTpdt)[0] = dlp;
(*dLagrangeTpdt)[1] = 0.0;
} else {
- // else friction exceeds value necessary, so stick
+ // friction exceeds value necessary to stick
// no changes to solution
+ assert(0.0 == slipRateMag);
} // if/else
} else {
// if in tension, then traction is zero.
@@ -1989,13 +1994,16 @@
// if in compression and no opening
const double frictionStress =
_friction->calcFriction(slipShearMag, slipRateMag, tractionNormal);
- if (tractionShearMag > frictionStress) {
- // traction is limited by friction, so have sliding
- // Update slip based on value required to stick versus friction
+ if (tractionShearMag > frictionStress || slipRateMag > 0.0) {
+ // traction is limited by friction, so have sliding OR
+ // friction exceeds traction due to overshoot in slip
+
+ // Update traction increment based on value required to stick
+ // versus friction
const double dlp = -(tractionShearMag - frictionStress) * area *
- tractionTpdt[0] / tractionShearMag;
+ tractionTpdt[0] / tractionShearMag;
const double dlq = -(tractionShearMag - frictionStress) * area *
- tractionTpdt[1] / tractionShearMag;
+ tractionTpdt[1] / tractionShearMag;
(*dLagrangeTpdt)[0] = dlp;
(*dLagrangeTpdt)[1] = dlq;
@@ -2004,6 +2012,7 @@
} else {
// else friction exceeds value necessary, so stick
// no changes to solution
+ assert(0.0 == slipRateMag);
} // if/else
} else {
// if in tension, then traction is zero.
Modified: short/3D/PyLith/trunk/libsrc/pylith/friction/RateStateAgeing.cc
===================================================================
--- short/3D/PyLith/trunk/libsrc/pylith/friction/RateStateAgeing.cc 2011-05-26 16:13:29 UTC (rev 18466)
+++ short/3D/PyLith/trunk/libsrc/pylith/friction/RateStateAgeing.cc 2011-05-26 18:27:00 UTC (rev 18467)
@@ -356,18 +356,22 @@
// thetaTpdt = thetaT * exp(-slipRate/L * dt)
// + dt - 0.5*(sliprate/L)*dt**2 + 1.0/6.0*(slipRate/L)*dt**3;
+ // Since regulatized friction -> 0 as slipRate -> 0, limit slip
+ // rate to some minimum value
+ const double slipRateEff = std::max(1.0e-12, slipRate);
+
const double dt = _dt;
const double thetaTVertex = stateVars[s_state];
const double L = properties[p_L];
- const double vDtL = slipRate * dt / L;
+ const double vDtL = slipRateEff * dt / L;
const double expTerm = exp(-vDtL);
double thetaTpdtVertex = 0.0;
if (vDtL > 1.0e-20) {
- thetaTpdtVertex = thetaTVertex * expTerm + L / slipRate * (1 - expTerm);
+ thetaTpdtVertex = thetaTVertex * expTerm + L / slipRateEff * (1 - expTerm);
PetscLogFlops(7);
} else {
- thetaTpdtVertex = thetaTVertex * expTerm + dt - 0.5 * slipRate/L * dt*dt;
+ thetaTpdtVertex = thetaTVertex * expTerm + dt - 0.5 * slipRateEff/L * dt*dt;
PetscLogFlops(9);
} // if/else
More information about the CIG-COMMITS
mailing list