[cig-commits] r18348 - in short/3D/PyLith/trunk: . examples libsrc/faults libsrc/friction libsrc/problems tests/2d
brad at geodynamics.org
brad at geodynamics.org
Wed May 11 17:23:11 PDT 2011
Author: brad
Date: 2011-05-11 17:23:10 -0700 (Wed, 11 May 2011)
New Revision: 18348
Added:
short/3D/PyLith/trunk/examples/2d/
short/3D/PyLith/trunk/tests/2d/frictionslide/
Modified:
short/3D/PyLith/trunk/configure.ac
short/3D/PyLith/trunk/examples/Makefile.am
short/3D/PyLith/trunk/libsrc/faults/FaultCohesiveDyn.cc
short/3D/PyLith/trunk/libsrc/friction/RateStateAgeing.cc
short/3D/PyLith/trunk/libsrc/problems/SolverNonlinear.cc
Log:
Merge from stable.
Modified: short/3D/PyLith/trunk/configure.ac
===================================================================
--- short/3D/PyLith/trunk/configure.ac 2011-05-11 23:48:12 UTC (rev 18347)
+++ short/3D/PyLith/trunk/configure.ac 2011-05-12 00:23:10 UTC (rev 18348)
@@ -351,11 +351,14 @@
doc/presentations/Makefile
examples/Makefile
examples/3d/Makefile
- examples/bar_shearwave/Makefile
- examples/twocells/Makefile
examples/3d/tet4/Makefile
examples/3d/hex8/Makefile
examples/3d/hex8/output/Makefile
+ examples/2d/Makefile
+ examples/2d/subduction/Makefile
+ examples/2d/subduction/output/Makefile
+ examples/bar_shearwave/Makefile
+ examples/twocells/Makefile
examples/bar_shearwave/hex8/Makefile
examples/bar_shearwave/hex8/output/Makefile
examples/bar_shearwave/quad4/Makefile
Copied: short/3D/PyLith/trunk/examples/2d (from rev 18347, short/3D/PyLith/branches/v1.5-stable/examples/2d)
Modified: short/3D/PyLith/trunk/examples/Makefile.am
===================================================================
--- short/3D/PyLith/trunk/examples/Makefile.am 2011-05-11 23:48:12 UTC (rev 18347)
+++ short/3D/PyLith/trunk/examples/Makefile.am 2011-05-12 00:23:10 UTC (rev 18348)
@@ -17,6 +17,7 @@
#
SUBDIRS = \
+ 2d \
3d \
bar_shearwave \
greensfns \
Modified: short/3D/PyLith/trunk/libsrc/faults/FaultCohesiveDyn.cc
===================================================================
--- short/3D/PyLith/trunk/libsrc/faults/FaultCohesiveDyn.cc 2011-05-11 23:48:12 UTC (rev 18347)
+++ short/3D/PyLith/trunk/libsrc/faults/FaultCohesiveDyn.cc 2011-05-12 00:23:10 UTC (rev 18348)
@@ -1900,7 +1900,7 @@
} // _sensitivitySolveLumped3D
// ----------------------------------------------------------------------
-// Constrain solution space with lumped Jacobian in 1-D.
+// Constrain solution space in 1-D.
void
pylith::faults::FaultCohesiveDyn::_constrainSolnSpace1D(double_array* dLagrangeTpdt,
const double_array& slip,
@@ -1923,7 +1923,7 @@
} // _constrainSolnSpace1D
// ----------------------------------------------------------------------
-// Constrain solution space with lumped Jacobian in 2-D.
+// Constrain solution space in 2-D.
void
pylith::faults::FaultCohesiveDyn::_constrainSolnSpace2D(double_array* dLagrangeTpdt,
const double_array& slip,
@@ -1948,7 +1948,7 @@
// Update slip based on value required to stick versus friction
const double dlp = -(tractionShearMag - frictionStress) * area *
- tractionTpdt[0] / tractionShearMag;
+ tractionTpdt[0] / tractionShearMag;
(*dLagrangeTpdt)[0] = dlp;
(*dLagrangeTpdt)[1] = 0.0;
} else {
@@ -1965,7 +1965,7 @@
} // _constrainSolnSpace2D
// ----------------------------------------------------------------------
-// Constrain solution space with lumped Jacobian in 3-D.
+// Constrain solution space in 3-D.
void
pylith::faults::FaultCohesiveDyn::_constrainSolnSpace3D(double_array* dLagrangeTpdt,
const double_array& slip,
Modified: short/3D/PyLith/trunk/libsrc/friction/RateStateAgeing.cc
===================================================================
--- short/3D/PyLith/trunk/libsrc/friction/RateStateAgeing.cc 2011-05-11 23:48:12 UTC (rev 18347)
+++ short/3D/PyLith/trunk/libsrc/friction/RateStateAgeing.cc 2011-05-12 00:23:10 UTC (rev 18348)
@@ -294,9 +294,14 @@
double mu_f = 0.0;
if (normalTraction <= 0.0) {
// if fault is in compression
- // Using Regularized Rate and State equation
+
+ // regularized rate and state equation
const double f0 = properties[p_coef];
+ // 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 slipRate0 = properties[p_slipRate0];
const double a = properties[p_a];
@@ -305,7 +310,7 @@
const double b = properties[p_b];
const double bLnTerm = b * log(slipRate0 * theta / L);
const double expTerm = exp((f0 + bLnTerm)/a);
- const double sinhArg = 0.5 * slipRate / slipRate0 * expTerm;
+ const double sinhArg = 0.5 * slipRateEff / slipRate0 * expTerm;
mu_f = a * asinh(sinhArg);
friction = -mu_f * normalTraction + properties[p_cohesion];
@@ -332,31 +337,42 @@
assert(numStateVars);
assert(_RateStateAgeing::numStateVars == numStateVars);
+ // 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
Modified: short/3D/PyLith/trunk/libsrc/problems/SolverNonlinear.cc
===================================================================
--- short/3D/PyLith/trunk/libsrc/problems/SolverNonlinear.cc 2011-05-11 23:48:12 UTC (rev 18347)
+++ short/3D/PyLith/trunk/libsrc/problems/SolverNonlinear.cc 2011-05-12 00:23:10 UTC (rev 18348)
@@ -210,6 +210,10 @@
Formulation* formulation = (Formulation*) context;
assert(0 != formulation);
+
+ // TEMPORARY - ASK MATT
+ formulation->constrainSolnSpace(&tmpSolutionVec);
+
// Reform residual
formulation->reformResidual(&tmpResidualVec, &tmpSolutionVec);
Copied: short/3D/PyLith/trunk/tests/2d/frictionslide (from rev 18347, short/3D/PyLith/branches/v1.5-stable/tests/2d/frictionslide)
More information about the CIG-COMMITS
mailing list