[cig-commits] r22189 - short/3D/PyLith/branches/v1.7-trunk/libsrc/pylith/faults
brad at geodynamics.org
brad at geodynamics.org
Fri Jun 7 17:06:59 PDT 2013
Author: brad
Date: 2013-06-07 17:06:59 -0700 (Fri, 07 Jun 2013)
New Revision: 22189
Modified:
short/3D/PyLith/branches/v1.7-trunk/libsrc/pylith/faults/FaultCohesiveDyn.cc
short/3D/PyLith/branches/v1.7-trunk/libsrc/pylith/faults/FaultCohesiveDyn.hh
Log:
Added patch from trunk- Newton friction solve.
Modified: short/3D/PyLith/branches/v1.7-trunk/libsrc/pylith/faults/FaultCohesiveDyn.cc
===================================================================
--- short/3D/PyLith/branches/v1.7-trunk/libsrc/pylith/faults/FaultCohesiveDyn.cc 2013-06-07 23:46:33 UTC (rev 22188)
+++ short/3D/PyLith/branches/v1.7-trunk/libsrc/pylith/faults/FaultCohesiveDyn.cc 2013-06-08 00:06:59 UTC (rev 22189)
@@ -48,10 +48,6 @@
#include <sstream> // USES std::ostringstream
#include <stdexcept> // USES std::runtime_error
-// Precomputing geometry significantly increases storage but gives a
-// slight speed improvement.
-//#define PRECOMPUTE_GEOMETRY
-
// ----------------------------------------------------------------------
typedef pylith::topology::Mesh::SieveMesh SieveMesh;
typedef pylith::topology::Mesh::RealSection RealSection;
@@ -545,6 +541,7 @@
const scalar_array&,
const scalar_array&,
const scalar_array&,
+ const PylithScalar,
const bool);
assert(fields);
@@ -735,12 +732,13 @@
// Use fault constitutive model to compute traction associated with
// friction.
dTractionTpdtVertex = 0.0;
+ const PylithScalar jacobianShearVertex = 0.0;
const bool iterating = true; // Iterating to get friction
CALL_MEMBER_FN(*this,
constrainSolnSpaceFn)(&dTractionTpdtVertex,
t, slipTpdtVertex, slipRateVertex,
tractionTpdtVertex,
- iterating);
+ jacobianShearVertex, iterating);
// Rotate increment in traction back to global coordinate system.
dLagrangeTpdtVertex = 0.0;
@@ -1172,6 +1170,7 @@
const scalar_array&,
const scalar_array&,
const scalar_array&,
+ const PylithScalar,
const bool);
assert(fields);
@@ -1375,6 +1374,7 @@
slipVertex = 0.0;
slipRateVertex = 0.0;
tractionTpdtVertex = 0.0;
+ PylithScalar jacobianShearVertex = 0.0;
for (int iDim=0; iDim < spaceDim; ++iDim) {
for (int jDim=0; jDim < spaceDim; ++jDim) {
slipVertex[iDim] += orientationVertex[iDim*spaceDim+jDim] *
@@ -1383,6 +1383,7 @@
velRelVertex[jDim];
tractionTpdtVertex[iDim] += orientationVertex[iDim*spaceDim+jDim] *
(lagrangeTVertex[jDim] + lagrangeTIncrVertex[jDim]);
+ jacobianShearVertex += orientationVertex[iDim*spaceDim+jDim] * (-1.0 / (areaVertex * (1.0 / jacobianVertexN[iDim] + 1.0 / jacobianVertexP[iDim])));
} // for
} // for
@@ -1397,7 +1398,7 @@
constrainSolnSpaceFn)(&dTractionTpdtVertex,
t, slipVertex, slipRateVertex,
tractionTpdtVertex,
- iterating);
+ jacobianShearVertex, iterating);
// Rotate traction back to global coordinate system.
dLagrangeTpdtVertex = 0.0;
@@ -2247,6 +2248,7 @@
const scalar_array&,
const scalar_array&,
const scalar_array&,
+ const PylithScalar,
const bool);
// Update time step in friction (can vary).
@@ -2431,12 +2433,13 @@
// Use fault constitutive model to compute traction associated with
// friction.
tractionMisfitVertex = 0.0;
+ const PylithScalar jacobianShearVertex = 0.0;
const bool iterating = true; // Iterating to get friction
CALL_MEMBER_FN(*this,
constrainSolnSpaceFn)(&tractionMisfitVertex, t,
slipTpdtVertex, slipRateVertex,
tractionTpdtVertex,
- iterating);
+ jacobianShearVertex, iterating);
#if 0 // DEBUGGING
std::cout << "alpha: " << alpha
@@ -2491,6 +2494,7 @@
const scalar_array& slip,
const scalar_array& sliprate,
const scalar_array& tractionTpdt,
+ const PylithScalar jacobianShear,
const bool iterating)
{ // _constrainSolnSpace1D
assert(dTractionTpdt);
@@ -2515,11 +2519,12 @@
const scalar_array& slip,
const scalar_array& slipRate,
const scalar_array& tractionTpdt,
+ const PylithScalar jacobianShear,
const bool iterating)
{ // _constrainSolnSpace2D
assert(dTractionTpdt);
- const PylithScalar slipMag = fabs(slip[0]);
+ PylithScalar slipMag = fabs(slip[0]);
const PylithScalar slipRateMag = fabs(slipRate[0]);
const PylithScalar tractionNormal = tractionTpdt[1];
@@ -2527,17 +2532,44 @@
if (fabs(slip[1]) < _zeroTolerance && tractionNormal < -_zeroTolerance) {
// if in compression and no opening
- const PylithScalar frictionStress =
- _friction->calcFriction(t, slipMag, slipRateMag, tractionNormal);
+ PylithScalar frictionStress = _friction->calcFriction(t, slipMag, slipRateMag, tractionNormal);
+
if (tractionShearMag > frictionStress || (iterating && slipRateMag > 0.0)) {
// traction is limited by friction, so have sliding OR
// friction exceeds traction due to overshoot in slip
if (tractionShearMag > 0.0) {
+#if 1 // New Newton stuff
+ if (fabs(jacobianShear) > 0.0) {
+ // Use Newton to get better update
+ const int maxiter = 16;
+ PylithScalar slipMagCur = slipMag;
+ PylithScalar slipRateMagCur = slipRateMag;
+ PylithScalar tractionShearMagCur = tractionShearMag;
+ for (int iter=0; iter < maxiter; ++iter) {
+ const PylithScalar frictionDeriv = _friction->calcFrictionDeriv(t, slipMagCur, slipRateMagCur, tractionNormal);
+ slipMag = slipMagCur;
+ if (slipMag > 0.0) {
+ // Use Newton (in log slip space) to get better update in slip & traction.
+ // D_{i+1} = exp(ln(D_i) - (T-T_f)/(D_i * (jacobian - frictionDeriv))
+ slipMagCur = exp(log(slipMag) - (tractionShearMagCur - frictionStress) / (slipMag * (jacobianShear - frictionDeriv)));
+ } else {
+ // Use Newton (in linear slip space) to get better update in slip & traction.
+ // D_{i+1} = D_i - (T-T_f)/(jacobian - frictionDeriv)
+ slipMagCur = slipMag - (tractionShearMagCur - frictionStress) / (jacobianShear - frictionDeriv);
+ } // if
+ tractionShearMagCur += (slipMagCur - slipMag) * jacobianShear;
+ slipRateMagCur += (slipMagCur - slipMag) / _dt;
+ frictionStress = _friction->calcFriction(t, slipMagCur, slipRateMagCur, tractionNormal);
+ if (fabs(tractionShearMagCur - frictionStress) < _zeroTolerance) {
+ break;
+ } // if
+ } // for
+ } // if
+#endif
// Update traction increment based on value required to stick
// versus friction
- const PylithScalar dlp = -(tractionShearMag - frictionStress) *
- tractionTpdt[0] / tractionShearMag;
+ const PylithScalar dlp = -(tractionShearMag - frictionStress) * tractionTpdt[0] / tractionShearMag;
(*dTractionTpdt)[0] = dlp;
} else {
// No shear stress and no friction.
@@ -2566,36 +2598,58 @@
const scalar_array& slip,
const scalar_array& slipRate,
const scalar_array& tractionTpdt,
+ const PylithScalar jacobianShear,
const bool iterating)
{ // _constrainSolnSpace3D
assert(dTractionTpdt);
- const PylithScalar slipShearMag = sqrt(slip[0] * slip[0] +
- slip[1] * slip[1]);
- PylithScalar slipRateMag = sqrt(slipRate[0]*slipRate[0] +
- slipRate[1]*slipRate[1]);
+ PylithScalar slipMag = sqrt(slip[0] * slip[0] + slip[1] * slip[1]);
+ const PylithScalar slipRateMag = sqrt(slipRate[0]*slipRate[0] + slipRate[1]*slipRate[1]);
const PylithScalar tractionNormal = tractionTpdt[2];
- const PylithScalar tractionShearMag =
- sqrt(tractionTpdt[0] * tractionTpdt[0] +
- tractionTpdt[1] * tractionTpdt[1]);
+ const PylithScalar tractionShearMag = sqrt(tractionTpdt[0] * tractionTpdt[0] + tractionTpdt[1] * tractionTpdt[1]);
if (fabs(slip[2]) < _zeroTolerance && tractionNormal < -_zeroTolerance) {
// if in compression and no opening
- const PylithScalar frictionStress =
- _friction->calcFriction(t, slipShearMag, slipRateMag, tractionNormal);
+ PylithScalar frictionStress = _friction->calcFriction(t, slipMag, slipRateMag, tractionNormal);
if (tractionShearMag > frictionStress || (iterating && slipRateMag > 0.0)) {
// traction is limited by friction, so have sliding OR
// friction exceeds traction due to overshoot in slip
if (tractionShearMag > 0.0) {
+#if 1 // New Newton stuff
+ if (fabs(jacobianShear) > 0.0) {
+ // Use Newton to get better update
+ const int maxiter = 16;
+ PylithScalar slipMagCur = slipMag;
+ PylithScalar slipRateMagCur = slipRateMag;
+ PylithScalar tractionShearMagCur = tractionShearMag;
+ for (int iter=0; iter < maxiter; ++iter) {
+ const PylithScalar frictionDeriv = _friction->calcFrictionDeriv(t, slipMagCur, slipRateMagCur, tractionNormal);
+ slipMag = slipMagCur;
+ if (slipMag > 0.0) {
+ // Use Newton (in log slip space) to get better update in slip & traction.
+ // D_{i+1} = exp(ln(D_i) - (T-T_f)/(D_i * (jacobian - frictionDeriv))
+ slipMagCur = exp(log(slipMag) - (tractionShearMagCur - frictionStress) / (slipMag * (jacobianShear - frictionDeriv)));
+ } else {
+ // Use Newton (in linear slip space) to get better update in slip & traction.
+ // D_{i+1} = D_i - (T-T_f)/(jacobian - frictionDeriv)
+ slipMagCur = slipMag - (tractionShearMagCur - frictionStress) / (jacobianShear - frictionDeriv);
+ } // if
+ tractionShearMagCur += (slipMagCur - slipMag) * jacobianShear;
+ slipRateMagCur += (slipMagCur - slipMag) / _dt;
+ frictionStress = _friction->calcFriction(t, slipMagCur, slipRateMagCur, tractionNormal);
+ if (fabs(tractionShearMagCur - frictionStress) < _zeroTolerance) {
+ break;
+ } // if
+ } // for
+ } // if
+#endif
// Update traction increment based on value required to stick
// versus friction
- const PylithScalar dlp = -(tractionShearMag - frictionStress) *
- tractionTpdt[0] / tractionShearMag;
- const PylithScalar dlq = -(tractionShearMag - frictionStress) *
- tractionTpdt[1] / tractionShearMag;
+ const PylithScalar dlp = -(tractionShearMag - frictionStress) * tractionTpdt[0] / tractionShearMag;
+ const PylithScalar dlq = -(tractionShearMag - frictionStress) * tractionTpdt[1] / tractionShearMag;
(*dTractionTpdt)[0] = dlp;
(*dTractionTpdt)[1] = dlq;
Modified: short/3D/PyLith/branches/v1.7-trunk/libsrc/pylith/faults/FaultCohesiveDyn.hh
===================================================================
--- short/3D/PyLith/branches/v1.7-trunk/libsrc/pylith/faults/FaultCohesiveDyn.hh 2013-06-07 23:46:33 UTC (rev 22188)
+++ short/3D/PyLith/branches/v1.7-trunk/libsrc/pylith/faults/FaultCohesiveDyn.hh 2013-06-08 00:06:59 UTC (rev 22189)
@@ -230,12 +230,15 @@
* @param slip Slip assoc. w/Lagrange multiplier vertex.
* @param slipRate Slip rate assoc. w/Lagrange multiplier vertex.
* @param tractionTpdt Fault traction assoc. w/Lagrange multiplier vertex.
+ * @param jacobianShear Derivative of shear traction with respect to slip (elasticity).
+ * @param iterating True if iterating on solution.
*/
void _constrainSolnSpace1D(scalar_array* dLagrangeTpdt,
const PylithScalar t,
const scalar_array& slip,
const scalar_array& slipRate,
const scalar_array& tractionTpdt,
+ const PylithScalar jacobianShear,
const bool iterating =true);
/** Constrain solution space in 2-D.
@@ -245,12 +248,15 @@
* @param slip Slip assoc. w/Lagrange multiplier vertex.
* @param slipRate Slip rate assoc. w/Lagrange multiplier vertex.
* @param tractionTpdt Fault traction assoc. w/Lagrange multiplier vertex.
+ * @param jacobianShear Derivative of shear traction with respect to slip (elasticity).
+ * @param iterating True if iterating on solution.
*/
void _constrainSolnSpace2D(scalar_array* dLagrangeTpdt,
const PylithScalar t,
const scalar_array& slip,
const scalar_array& slipRate,
const scalar_array& tractionTpdt,
+ const PylithScalar jacobianShear,
const bool iterating =true);
/** Constrain solution space in 3-D.
@@ -260,12 +266,15 @@
* @param slip Slip assoc. w/Lagrange multiplier vertex.
* @param slipRate Slip rate assoc. w/Lagrange multiplier vertex.
* @param tractionTpdt Fault traction assoc. w/Lagrange multiplier vertex.
+ * @param jacobianShear Derivative of shear traction with respect to slip (elasticity).
+ * @param iterating True if iterating on solution.
*/
void _constrainSolnSpace3D(scalar_array* dLagrangeTpdt,
const PylithScalar t,
const scalar_array& slip,
const scalar_array& slipRate,
const scalar_array& tractionTpdt,
+ const PylithScalar jacobianShear,
const bool iterating =true);
// PRIVATE MEMBERS ////////////////////////////////////////////////////
More information about the CIG-COMMITS
mailing list