[cig-commits] r19640 - in short/3D/PyLith/branches/v1.6-stable: examples/3d/hex8 libsrc/pylith/faults libsrc/pylith/friction libsrc/pylith/problems libsrc/pylith/utils modulesrc/friction tests/2d/frictionslide tests/3d/cyclicfriction
brad at geodynamics.org
brad at geodynamics.org
Wed Feb 15 17:36:54 PST 2012
Author: brad
Date: 2012-02-15 17:36:53 -0800 (Wed, 15 Feb 2012)
New Revision: 19640
Modified:
short/3D/PyLith/branches/v1.6-stable/examples/3d/hex8/pylithapp.cfg
short/3D/PyLith/branches/v1.6-stable/examples/3d/hex8/test.cfg
short/3D/PyLith/branches/v1.6-stable/libsrc/pylith/faults/FaultCohesiveDyn.cc
short/3D/PyLith/branches/v1.6-stable/libsrc/pylith/faults/FaultCohesiveDyn.hh
short/3D/PyLith/branches/v1.6-stable/libsrc/pylith/friction/FrictionModel.cc
short/3D/PyLith/branches/v1.6-stable/libsrc/pylith/friction/FrictionModel.hh
short/3D/PyLith/branches/v1.6-stable/libsrc/pylith/friction/RateStateAgeing.cc
short/3D/PyLith/branches/v1.6-stable/libsrc/pylith/friction/RateStateAgeing.hh
short/3D/PyLith/branches/v1.6-stable/libsrc/pylith/friction/SlipWeakening.cc
short/3D/PyLith/branches/v1.6-stable/libsrc/pylith/friction/SlipWeakening.hh
short/3D/PyLith/branches/v1.6-stable/libsrc/pylith/friction/StaticFriction.cc
short/3D/PyLith/branches/v1.6-stable/libsrc/pylith/friction/StaticFriction.hh
short/3D/PyLith/branches/v1.6-stable/libsrc/pylith/friction/TimeWeakening.cc
short/3D/PyLith/branches/v1.6-stable/libsrc/pylith/friction/TimeWeakening.hh
short/3D/PyLith/branches/v1.6-stable/libsrc/pylith/problems/SolverNonlinear.cc
short/3D/PyLith/branches/v1.6-stable/libsrc/pylith/utils/constdefs.h
short/3D/PyLith/branches/v1.6-stable/modulesrc/friction/FrictionModel.i
short/3D/PyLith/branches/v1.6-stable/modulesrc/friction/RateStateAgeing.i
short/3D/PyLith/branches/v1.6-stable/modulesrc/friction/SlipWeakening.i
short/3D/PyLith/branches/v1.6-stable/modulesrc/friction/StaticFriction.i
short/3D/PyLith/branches/v1.6-stable/modulesrc/friction/TimeWeakening.i
short/3D/PyLith/branches/v1.6-stable/tests/2d/frictionslide/plot_friction.py
short/3D/PyLith/branches/v1.6-stable/tests/2d/frictionslide/pylithapp.cfg
short/3D/PyLith/branches/v1.6-stable/tests/3d/cyclicfriction/pylithapp.cfg
Log:
More work on quasi-static friction. Reworked constrainSolnSpace() to use line search in log space. Still need to fix fault opening with line search.
Modified: short/3D/PyLith/branches/v1.6-stable/examples/3d/hex8/pylithapp.cfg
===================================================================
--- short/3D/PyLith/branches/v1.6-stable/examples/3d/hex8/pylithapp.cfg 2012-02-15 22:48:27 UTC (rev 19639)
+++ short/3D/PyLith/branches/v1.6-stable/examples/3d/hex8/pylithapp.cfg 2012-02-16 01:36:53 UTC (rev 19640)
@@ -70,7 +70,6 @@
# PETSc
# ----------------------------------------------------------------------
# Set the solver options.
-
[pylithapp.petsc]
# Preconditioner settings.
Modified: short/3D/PyLith/branches/v1.6-stable/examples/3d/hex8/test.cfg
===================================================================
--- short/3D/PyLith/branches/v1.6-stable/examples/3d/hex8/test.cfg 2012-02-15 22:48:27 UTC (rev 19639)
+++ short/3D/PyLith/branches/v1.6-stable/examples/3d/hex8/test.cfg 2012-02-16 01:36:53 UTC (rev 19640)
@@ -142,7 +142,7 @@
friction.db_properties = spatialdata.spatialdb.UniformDB
friction.db_properties.label = Rate Stete Ageing
friction.db_properties.values = [reference-friction-coefficient,reference-slip-rate,characteristic-slip-distance,constitutive-parameter-a,constitutive-parameter-b,cohesion]
-friction.db_properties.data = [0.4, 1.0e-7*m/s, 2.0*m, 0.01, 0.02, 0.0*Pa]
+friction.db_properties.data = [0.4, 1.0e-7*m/s, 0.1*m, 0.01, 0.02, 0.0*Pa]
# Set spatial database for the initial value of the state variable.
friction.db_initial_state = spatialdata.spatialdb.UniformDB
Modified: short/3D/PyLith/branches/v1.6-stable/libsrc/pylith/faults/FaultCohesiveDyn.cc
===================================================================
--- short/3D/PyLith/branches/v1.6-stable/libsrc/pylith/faults/FaultCohesiveDyn.cc 2012-02-15 22:48:27 UTC (rev 19639)
+++ short/3D/PyLith/branches/v1.6-stable/libsrc/pylith/faults/FaultCohesiveDyn.cc 2012-02-16 01:36:53 UTC (rev 19640)
@@ -338,6 +338,13 @@
residualVertexN = areaVertex * (dispTpdtVertexL - initialTractionsVertex);
} else { // opening, normal traction should be zero
+ if (fabs(tractionNormal) > _zeroTolerance) {
+ std::cerr << "WARNING! Fault opening with nonzero traction."
+ << ", v_fault: " << v_fault
+ << ", opening: " << slipNormal
+ << ", normal traction: " << tractionNormal
+ << std::endl;
+ } // if
assert(fabs(tractionNormal) < _zeroTolerance);
} // if/else
residualVertexP = -residualVertexN;
@@ -356,12 +363,6 @@
residualSection->getFiberDimension(v_positive));
residualSection->updateAddPoint(v_positive, &residualVertexP[0]);
-#if 0 // TEST
- assert(residualVertexL.size() ==
- residualSection->getFiberDimension(v_lagrange));
- residualSection->updateAddPoint(v_lagrange, &residualVertexL[0]);
-#endif
-
#if defined(DETAILED_EVENT_LOGGING)
_logger->eventEnd(updateEvent);
#endif
@@ -510,18 +511,11 @@
const double_array&,
const double_array&,
const bool);
- /// Member prototype for _constrainSolnSpaceImproveXD()
- typedef void (pylith::faults::FaultCohesiveDyn::*constrainSolnSpaceImprove_fn_type)
- (double_array*,
- double_array*,
- const double_array&,
- const double_array&,
- const double_array&);
- assert(0 != fields);
- assert(0 != _quadrature);
- assert(0 != _fields);
- assert(0 != _friction);
+ assert(fields);
+ assert(_quadrature);
+ assert(_fields);
+ assert(_friction);
_sensitivitySetup(jacobian);
@@ -567,25 +561,18 @@
assert(!dLagrangeTpdtSection.isNull());
constrainSolnSpace_fn_type constrainSolnSpaceFn;
- constrainSolnSpaceImprove_fn_type constrainSolnSpaceImproveFn;
switch (spaceDim) { // switch
case 1:
constrainSolnSpaceFn =
&pylith::faults::FaultCohesiveDyn::_constrainSolnSpace1D;
- constrainSolnSpaceImproveFn =
- &pylith::faults::FaultCohesiveDyn::_constrainSolnSpaceImprove1D;
break;
case 2:
constrainSolnSpaceFn =
&pylith::faults::FaultCohesiveDyn::_constrainSolnSpace2D;
- constrainSolnSpaceImproveFn =
- &pylith::faults::FaultCohesiveDyn::_constrainSolnSpaceImprove2D;
break;
case 3:
constrainSolnSpaceFn =
&pylith::faults::FaultCohesiveDyn::_constrainSolnSpace3D;
- constrainSolnSpaceImproveFn =
- &pylith::faults::FaultCohesiveDyn::_constrainSolnSpaceImprove3D;
break;
default :
assert(0);
@@ -745,18 +732,9 @@
std::cout << ", tractionTpdtVertex: ";
for (int iDim=0; iDim < spaceDim; ++iDim)
std::cout << " " << tractionTpdtVertex[iDim];
- std::cout << ", lagrangeTVertex: ";
- for (int iDim=0; iDim < spaceDim; ++iDim)
- std::cout << " " << lagrangeTVertex[iDim];
- std::cout << ", lagrangeTIncrVertex: ";
- for (int iDim=0; iDim < spaceDim; ++iDim)
- std::cout << " " << lagrangeTIncrVertex[iDim];
std::cout << ", dTractionTpdtVertex: ";
for (int iDim=0; iDim < spaceDim; ++iDim)
std::cout << " " << dTractionTpdtVertex[iDim];
- std::cout << ", dLagrangeTpdtVertex: ";
- for (int iDim=0; iDim < spaceDim; ++iDim)
- std::cout << " " << dLagrangeTpdtVertex[iDim];
std::cout << std::endl;
#endif
@@ -794,30 +772,109 @@
// change in Lagrange multipliers imposed by friction criterion.
// Solve sensitivity problem for negative side of the fault.
- bool negativeSide = true;
- _sensitivityUpdateJacobian(negativeSide, jacobian, *fields);
- _sensitivityReformResidual(negativeSide);
+ bool negativeSideFlag = true;
+ _sensitivityUpdateJacobian(negativeSideFlag, jacobian, *fields);
+ _sensitivityReformResidual(negativeSideFlag);
_sensitivitySolve();
- _sensitivityUpdateSoln(negativeSide);
+ _sensitivityUpdateSoln(negativeSideFlag);
// Solve sensitivity problem for positive side of the fault.
- negativeSide = false;
- _sensitivityUpdateJacobian(negativeSide, jacobian, *fields);
- _sensitivityReformResidual(negativeSide);
+ negativeSideFlag = false;
+ _sensitivityUpdateJacobian(negativeSideFlag, jacobian, *fields);
+ _sensitivityReformResidual(negativeSideFlag);
_sensitivitySolve();
- _sensitivityUpdateSoln(negativeSide);
+ _sensitivityUpdateSoln(negativeSideFlag);
// Step 4: Update Lagrange multipliers and displacement fields based
// on changes imposed by friction criterion in Step 2 (change in
// Lagrange multipliers) and Step 3 (slip associated with change in
// Lagrange multipliers).
+ //
+ // Use line search to find best update. This improves convergence
+ // because it accounts for feedback between the fault constitutive
+ // model and the deformation. We also search in log space because
+ // some fault constitutive models depend on the log of slip rate.
+ const double residualTol = _zeroTolerance; // L2 misfit in tractions
+ const int maxIter = 32;
+ double logAlphaL = log10(_zeroTolerance); // minimum step
+ double logAlphaR = log10(1.0); // maximum step
+ double logAlphaM = 0.5*(logAlphaL + logAlphaR);
+ double logAlphaML = 0.5*(logAlphaL + logAlphaM);
+ double logAlphaMR = 0.5*(logAlphaM + logAlphaR);
+ double residualL = _constrainSolnSpaceNorm(pow(10.0, logAlphaL), fields);
+ double residualML = _constrainSolnSpaceNorm(pow(10.0, logAlphaML), fields);
+ double residualM = _constrainSolnSpaceNorm(pow(10.0, logAlphaM), fields);
+ double residualMR = _constrainSolnSpaceNorm(pow(10.0, logAlphaMR), fields);
+ double residualR = _constrainSolnSpaceNorm(pow(10.0, logAlphaR), fields);
+ for (int iter=0; iter < maxIter; ++iter) {
+ if (residualM < residualTol || residualR < _zeroTolerance)
+ break;
+
+ if (residualL < residualML && residualML < residualM) {
+ logAlphaL = logAlphaL;
+ logAlphaR = logAlphaM;
+ residualL = residualL;
+ residualR = residualM;
+ residualM = residualML;
+ } else if (residualL >= residualML && residualML <= residualM) {
+ logAlphaL = logAlphaL;
+ logAlphaR = logAlphaM;
+ residualL = residualL;
+ residualR = residualM;
+ residualM = residualML;
+ } else if (residualML >= residualM && residualM <= residualMR) {
+ logAlphaL = logAlphaML;
+ logAlphaR = logAlphaMR;
+ residualL = residualML;
+ residualR = residualMR;
+ residualM = residualM;
+ } else if (residualM >= residualMR && residualMR <= residualR) {
+ logAlphaL = logAlphaM;
+ logAlphaR = logAlphaR;
+ residualL = residualM;
+ residualR = residualR;
+ residualM = residualMR;
+ } else if (residualM > residualMR && residualMR > residualR) {
+ logAlphaL = logAlphaM;
+ logAlphaR = logAlphaR;
+ residualL = residualM;
+ residualR = residualR;
+ residualM = residualMR;
+ } else {
+ assert(0);
+ throw std::logic_error("Unknown case in constrain solution space "
+ "line search.");
+ } // if/else
+ logAlphaM = (logAlphaL + logAlphaR) / 2.0;
+ logAlphaML = (logAlphaL + logAlphaM) / 2.0;
+ logAlphaMR = (logAlphaM + logAlphaR) / 2.0;
+
+ residualML = _constrainSolnSpaceNorm(pow(10.0, logAlphaML), fields);
+ residualMR = _constrainSolnSpaceNorm(pow(10.0, logAlphaMR), fields);
+
+ } // for
+ // Account for possibility that end points have lowest residual
+ if (residualR < residualM) {
+ logAlphaM = logAlphaR;
+ residualM = residualR;
+ } else if (residualL < residualM) {
+ logAlphaM = logAlphaL;
+ residualM = residualL;
+ } // if/else
+ const double alpha = pow(10.0, logAlphaM); // alphaM is our best guess
+#if 1
+ std::cout << "ALPHA: " << alpha
+ << ", residual: " << residualM
+ << std::endl;
+#endif
+
double_array slipTVertex(spaceDim);
double_array dSlipTpdtVertex(spaceDim);
double_array dispRelVertex(spaceDim);
+ const ALE::Obj<RealSection>& sensDispRelSection = _fields->get("sensitivity relative disp").section();
+ assert(!sensDispRelSection.isNull());
- const ALE::Obj<RealSection>& sensDispRelSection =
- _fields->get("sensitivity relative disp").section();
for (int iVertex=0; iVertex < numVertices; ++iVertex) {
const int v_fault = _cohesiveVertices[iVertex].fault;
const int v_lagrange = _cohesiveVertices[iVertex].lagrange;
@@ -828,14 +885,6 @@
dLagrangeTpdtSection->restrictPoint(v_fault, &dLagrangeTpdtVertex[0],
dLagrangeTpdtVertex.size());
- // If no change in the Lagrange multiplier computed from friction
- // criterion, there are no updates, so continue.
- double dLagrangeTpdtMag = 0.0;
- for (int i=0; i < spaceDim; ++i)
- dLagrangeTpdtMag += dLagrangeTpdtVertex[i]*dLagrangeTpdtVertex[i];
- if (0.0 == dLagrangeTpdtMag)
- continue;
-
// Get change in relative displacement from sensitivity solve.
assert(spaceDim == sensDispRelSection->getFiberDimension(v_fault));
const double* sensDispRelVertex =
@@ -883,11 +932,8 @@
dispIncrSection->restrictPoint(v_lagrange);
assert(lagrangeTIncrVertex);
- // Step 4a: Prevent nonphysical trial solutions. The product of the
- // normal traction and normal slip must be nonnegative (forbid
- // interpenetration with tension or opening with compression).
-
- // Compute slip, change in slip, and tractions in fault coordinates.
+ // Scale perturbation in relative displacements and change in
+ // Lagrange multipliers by alpha using only shear components.
slipTVertex = 0.0;
slipTpdtVertex = 0.0;
dSlipTpdtVertex = 0.0;
@@ -908,6 +954,12 @@
dLagrangeTpdtVertex[jDim];
} // for
} // for
+ // Scale shear components of updates by alpha; leave normal
+ // components unchanged to get full updates.
+ for (int iDim=0; iDim < indexN; ++iDim) {
+ dSlipTpdtVertex[iDim] *= alpha;
+ dTractionTpdtVertex[iDim] *= alpha;
+ } // for
if (fabs(slipTpdtVertex[indexN]) < _zeroTolerance) {
slipTpdtVertex[indexN] = 0.0;
} // if
@@ -915,11 +967,15 @@
dSlipTpdtVertex[indexN] = 0.0;
} // if
+ // Step 5a: Prevent nonphysical trial solutions. The product of the
+ // normal traction and normal slip must be nonnegative (forbid
+ // interpenetration with tension or opening with compression).
+
if ((slipTpdtVertex[indexN] + dSlipTpdtVertex[indexN]) *
(tractionTpdtVertex[indexN] + dTractionTpdtVertex[indexN])
< 0.0) {
#if 0 // DEBUGGING
- std::cout << "STEP 4a CORRECTING NONPHYSICAL SLIP/TRACTIONS"
+ std::cout << "STEP 5a CORRECTING NONPHYSICAL SLIP/TRACTIONS"
<< ", v_fault: " << v_fault
<< ", slipNormal: " << slipTpdtVertex[indexN] + dSlipTpdtVertex[indexN]
<< ", tractionNormal: " << tractionTpdtVertex[indexN] + dTractionTpdtVertex[indexN]
@@ -941,7 +997,7 @@
// Do not allow fault interpenetration.
if (slipTpdtVertex[indexN] + dSlipTpdtVertex[indexN] < 0.0) {
#if 0 // DEBUGGING
- std::cout << "STEP 4a CORRECTING INTERPENETATION"
+ std::cout << "STEP 5a CORRECTING INTERPENETATION"
<< ", v_fault: " << v_fault
<< ", slipNormal: " << slipTpdtVertex[indexN] + dSlipTpdtVertex[indexN]
<< std::endl;
@@ -949,21 +1005,10 @@
dSlipTpdtVertex[indexN] = -slipTpdtVertex[indexN];
} // if
-
- // Improve estimate of slip and change in traction using dFriction/dD.
- // Get friction properties and state variables.
- _friction->retrievePropsStateVars(v_fault);
-
- CALL_MEMBER_FN(*this,
- constrainSolnSpaceImproveFn)(&dTractionTpdtVertex, &dSlipTpdtVertex,
- slipTVertex, slipTpdtVertex,
- tractionTpdtVertex);
-
-#if 0 // OBSOLETE? Move to ImproveFn?
- // Prevent over-correction in slip resulting in backslip
+ // Prevent over-correction in shear slip resulting in backslip
double slipDot = 0.0;
double tractionSlipDot = 0.0;
- for (int iDim=0; iDim < spaceDim-1; ++iDim) { // :TODO: Optimize this
+ for (int iDim=0; iDim < indexN; ++iDim) {
// Compute dot product between slip and increment in slip (want positive)
slipDot +=
(slipTpdtVertex[iDim] - slipTVertex[iDim]) *
@@ -975,14 +1020,19 @@
if (slipDot < 0.0 &&
sqrt(fabs(slipDot)) > _zeroTolerance &&
tractionSlipDot < 0.0) {
- // Correct backslip
- dTractionTpdtVertex *= 0.5; // Use bisection as guess for traction
- for (int iDim=0; iDim < spaceDim-1; ++iDim) {
- // Use bisection for slip
+#if 1 // DEBUGGING
+ std::cout << "STEP 5b CORRECTING BACKSLIP"
+ << ", v_fault: " << v_fault
+ << ", slipDot: " << slipDot
+ << ", tractionSlipDot: " << tractionSlipDot
+ << std::endl;
+#endif
+ // Correct backslip, use bisection as guess
+ for (int iDim=0; iDim < indexN; ++iDim) {
+ dTractionTpdtVertex[iDim] *= 0.5;
dSlipTpdtVertex[iDim] = 0.5*(slipTVertex[iDim] - slipTpdtVertex[iDim]);
} // for
} // if/else
-#endif
// Update current estimate of slip from t to t+dt.
slipTpdtVertex += dSlipTpdtVertex;
@@ -1022,12 +1072,6 @@
std::cout << std::endl;
#endif
-#if 0 // TEST
- // Set change in relative displacement.
- assert(dispRelVertex.size() ==
- dispRelSection->getFiberDimension(v_fault));
- dispRelSection->updatePoint(v_fault, &dispRelVertex[0]);
-#endif
// Update Lagrange multiplier increment.
assert(dLagrangeTpdtVertex.size() ==
dispIncrSection->getFiberDimension(v_lagrange));
@@ -1345,10 +1389,6 @@
dispIncrVertexP[iDim] -=
areaVertex * dLagrangeTpdtVertexGlobal[iDim] / jacobianVertexP[iDim];
- // Set increment in relative displacement.
- dispRelVertex[iDim] = -areaVertex * 2.0*dLagrangeTpdtVertexGlobal[iDim] /
- (jacobianVertexN[iDim] + jacobianVertexP[iDim]);
-
// Update increment in Lagrange multiplier.
lagrangeTIncrVertex[iDim] += dLagrangeTpdtVertexGlobal[iDim];
} // for
@@ -1381,12 +1421,6 @@
dispIncrSection->getFiberDimension(v_lagrange));
dispIncrSection->updatePoint(v_lagrange, &lagrangeTIncrVertex[0]);
- // Update the relative displacement estimate based on adjustment
- // to the Lagrange multiplier values.
- assert(dispRelVertex.size() ==
- dispRelSection->getFiberDimension(v_fault));
- dispRelSection->updateAddPoint(v_fault, &dispRelVertex[0]);
-
#if defined(DETAILED_EVENT_LOGGING)
_logger->eventEnd(updateEvent);
#endif
@@ -1402,8 +1436,6 @@
#if 0 // DEBUGGING
//dLagrangeTpdtSection->view("AFTER dLagrange");
//dispIncrSection->view("AFTER DISP INCR (t->t+dt)");
- dispRelSection->view("AFTER RELATIVE DISPLACEMENT");
- //velRelSection->view("AFTER RELATIVE VELOCITY");
#endif
} // adjustSolnLumped
@@ -2191,8 +2223,13 @@
double_array dispVertex(spaceDim);
const ALE::Obj<RealSection>& solutionSection =
_fields->get("sensitivity solution").section();
+ assert(!solutionSection.isNull());
const ALE::Obj<RealSection>& dispRelSection =
_fields->get("sensitivity relative disp").section();
+ assert(!dispRelSection.isNull());
+ const ALE::Obj<RealSection>& dLagrangeTpdtSection =
+ _fields->get("sensitivity dLagrange").section();
+ assert(!dLagrangeTpdtSection.isNull());
const double sign = (negativeSide) ? -1.0 : 1.0;
@@ -2202,6 +2239,20 @@
solutionSection->restrictPoint(v_fault, &dispVertex[0], dispVertex.size());
+ // If no change in the Lagrange multiplier computed from friction
+ // criterion, there are no updates, so continue.
+ assert(spaceDim == dLagrangeTpdtSection->getFiberDimension(v_fault));
+ const double* dLagrangeTpdtVertex = dLagrangeTpdtSection->restrictPoint(v_fault);
+ assert(dLagrangeTpdtVertex);
+ double dLagrangeTpdtVertexMag = 0.0;
+ for (int iDim=0; iDim < spaceDim; ++iDim) {
+ dLagrangeTpdtVertexMag += dLagrangeTpdtVertex[iDim]*dLagrangeTpdtVertex[iDim];
+ } // for
+ if (0.0 == dLagrangeTpdtVertexMag)
+ continue;
+
+ // Update relative displacements associated with sensitivity solve
+ // solution
dispVertex *= sign;
assert(dispVertex.size() == dispRelSection->getFiberDimension(v_fault));
@@ -2209,7 +2260,244 @@
} // for
} // _sensitivityUpdateSoln
+
// ----------------------------------------------------------------------
+// Compute norm of residual associated with matching fault
+// constitutive model using update from sensitivity solve. We use
+// this in a line search to find a good update (required because
+// fault constitutive model may have a complex nonlinear feedback
+// with deformation).
+double
+pylith::faults::FaultCohesiveDyn::_constrainSolnSpaceNorm(const double alpha,
+ topology::SolutionFields* const fields)
+{ // _constrainSolnSpaceNorm
+ /// Member prototype for _constrainSolnSpaceXD()
+ typedef void (pylith::faults::FaultCohesiveDyn::*constrainSolnSpace_fn_type)
+ (double_array*,
+ const double_array&,
+ const double_array&,
+ const double_array&,
+ const bool);
+
+ // Update time step in friction (can vary).
+ _friction->timeStep(_dt);
+ const double dt = _dt;
+
+ const int spaceDim = _quadrature->spaceDim();
+ const int indexN = spaceDim - 1;
+
+ constrainSolnSpace_fn_type constrainSolnSpaceFn;
+ switch (spaceDim) { // switch
+ case 1:
+ constrainSolnSpaceFn =
+ &pylith::faults::FaultCohesiveDyn::_constrainSolnSpace1D;
+ break;
+ case 2:
+ constrainSolnSpaceFn =
+ &pylith::faults::FaultCohesiveDyn::_constrainSolnSpace2D;
+ break;
+ case 3:
+ constrainSolnSpaceFn =
+ &pylith::faults::FaultCohesiveDyn::_constrainSolnSpace3D;
+ break;
+ default :
+ assert(0);
+ throw std::logic_error("Unknown spatial dimension in "
+ "FaultCohesiveDyn::constrainSolnSpace().");
+ } // switch
+
+ // Get sections
+ double_array slipTpdtVertex(spaceDim); // fault coordinates
+ double_array dSlipTpdtVertex(spaceDim); // fault coordinates
+ double_array slipRateVertex(spaceDim); // fault coordinates
+ double_array dSlipRateVertex(spaceDim); // fault coordinates
+ double_array tractionTpdtVertex(spaceDim); // fault coordinates
+ double_array dTractionTpdtVertex(spaceDim); // fault coordinates
+ double_array tractionMisfitVertex(spaceDim); // fault coordinates
+
+ const ALE::Obj<RealSection>& orientationSection =
+ _fields->get("orientation").section();
+ assert(!orientationSection.isNull());
+
+ const ALE::Obj<RealSection>& dLagrangeTpdtSection =
+ _fields->get("sensitivity dLagrange").section();
+ assert(!dLagrangeTpdtSection.isNull());
+
+ const ALE::Obj<RealSection>& sensDispRelSection =
+ _fields->get("sensitivity relative disp").section();
+ assert(!sensDispRelSection.isNull());
+
+ const ALE::Obj<RealSection>& dispTSection = fields->get("disp(t)").section();
+ assert(!dispTSection.isNull());
+
+ const ALE::Obj<RealSection>& dispIncrSection =
+ fields->get("dispIncr(t->t+dt)").section();
+ assert(!dispIncrSection.isNull());
+
+ double norm2 = 0.0;
+ const int numVertices = _cohesiveVertices.size();
+ for (int iVertex=0; iVertex < numVertices; ++iVertex) {
+ const int v_lagrange = _cohesiveVertices[iVertex].lagrange;
+ const int v_fault = _cohesiveVertices[iVertex].fault;
+ const int v_negative = _cohesiveVertices[iVertex].negative;
+ const int v_positive = _cohesiveVertices[iVertex].positive;
+
+ // Get displacement values
+ assert(spaceDim == dispTSection->getFiberDimension(v_negative));
+ const double* dispTVertexN = dispTSection->restrictPoint(v_negative);
+ assert(dispTVertexN);
+
+ assert(spaceDim == dispTSection->getFiberDimension(v_positive));
+ const double* dispTVertexP = dispTSection->restrictPoint(v_positive);
+ assert(dispTVertexP);
+
+ // Get displacement increment values.
+ assert(spaceDim == dispIncrSection->getFiberDimension(v_negative));
+ const double* dispIncrVertexN =
+ dispIncrSection->restrictPoint(v_negative);
+ assert(dispIncrVertexN);
+
+ assert(spaceDim == dispIncrSection->getFiberDimension(v_positive));
+ const double* dispIncrVertexP =
+ dispIncrSection->restrictPoint(v_positive);
+ assert(dispIncrVertexP);
+
+ // Get orientation
+ assert(spaceDim*spaceDim == orientationSection->getFiberDimension(v_fault));
+ const double* orientationVertex =
+ orientationSection->restrictPoint(v_fault);
+
+ // Get change in relative displacement from sensitivity solve.
+ assert(spaceDim == sensDispRelSection->getFiberDimension(v_fault));
+ const double* dDispRelVertex =
+ sensDispRelSection->restrictPoint(v_fault);
+ assert(dDispRelVertex);
+
+ // Get Lagrange multiplier values from disp(t), and dispIncr(t->t+dt)
+ assert(spaceDim == dispTSection->getFiberDimension(v_lagrange));
+ const double* lagrangeTVertex = dispTSection->restrictPoint(v_lagrange);
+ assert(lagrangeTVertex);
+
+ assert(spaceDim == dispIncrSection->getFiberDimension(v_lagrange));
+ const double* lagrangeTIncrVertex =
+ dispIncrSection->restrictPoint(v_lagrange);
+ assert(lagrangeTIncrVertex);
+
+ assert(spaceDim == dLagrangeTpdtSection->getFiberDimension(v_fault));
+ const double* dLagrangeTpdtVertex =
+ dLagrangeTpdtSection->restrictPoint(v_fault);
+ assert(dLagrangeTpdtVertex);
+
+ // Step 1: Prevent nonphysical trial solutions. The product of the
+ // normal traction and normal slip must be nonnegative (forbid
+ // interpenetration with tension or opening with compression).
+
+ // Compute slip, slip rate, and Lagrange multiplier at time t+dt
+ // in fault coordinate system.
+ slipTpdtVertex = 0.0;
+ dSlipTpdtVertex = 0.0;
+ slipRateVertex = 0.0;
+ dSlipRateVertex = 0.0;
+ tractionTpdtVertex = 0.0;
+ dTractionTpdtVertex = 0.0;
+ for (int iDim=0; iDim < spaceDim; ++iDim) {
+ for (int jDim=0; jDim < spaceDim; ++jDim) {
+ slipTpdtVertex[iDim] += orientationVertex[iDim*spaceDim+jDim] *
+ (dispTVertexP[jDim] + dispIncrVertexP[jDim]
+ - dispTVertexN[jDim] - dispIncrVertexN[jDim]);
+ dSlipTpdtVertex[iDim] += orientationVertex[iDim*spaceDim+jDim] *
+ dDispRelVertex[jDim];
+ slipRateVertex[iDim] += orientationVertex[iDim*spaceDim+jDim] *
+ (dispIncrVertexP[jDim] - dispIncrVertexN[jDim]) / dt;
+ dSlipRateVertex[iDim] += orientationVertex[iDim*spaceDim+jDim] *
+ dDispRelVertex[jDim] / dt;
+ tractionTpdtVertex[iDim] += orientationVertex[iDim*spaceDim+jDim] *
+ (lagrangeTVertex[jDim] + lagrangeTIncrVertex[jDim]);
+ dTractionTpdtVertex[iDim] += orientationVertex[iDim*spaceDim+jDim] *
+ dLagrangeTpdtVertex[jDim];
+ } // for
+ if (fabs(slipRateVertex[iDim]) < _zeroTolerance) {
+ slipRateVertex[iDim] = 0.0;
+ } // if
+ if (fabs(dSlipRateVertex[iDim]) < _zeroTolerance) {
+ dSlipRateVertex[iDim] = 0.0;
+ } // if
+ } // for
+ if (fabs(slipTpdtVertex[indexN]) < _zeroTolerance) {
+ slipTpdtVertex[indexN] = 0.0;
+ } // if
+ if (fabs(dSlipTpdtVertex[indexN]) < _zeroTolerance) {
+ dSlipTpdtVertex[indexN] = 0.0;
+ } // if
+
+ // Project shear slip and traction for line search
+ for (int iDim=0; iDim < indexN; ++iDim) {
+ slipTpdtVertex[iDim] += alpha*dSlipTpdtVertex[iDim];
+ slipRateVertex[iDim] += alpha*dSlipRateVertex[iDim];
+ tractionTpdtVertex[iDim] += alpha*dTractionTpdtVertex[iDim];
+ } // for
+ // Normal slip and tractions get full update
+ slipTpdtVertex[indexN] += dSlipTpdtVertex[indexN];
+ slipRateVertex[indexN] += dSlipRateVertex[indexN];
+ tractionTpdtVertex[indexN] += dTractionTpdtVertex[indexN];
+
+ if (slipTpdtVertex[indexN]*tractionTpdtVertex[indexN] < 0.0) {
+#if 0 // DEBUGGING
+ std::cout << "STEP 10 CORRECTING NONPHYSICAL SLIP/TRACTIONS"
+ << ", v_fault: " << v_fault
+ << ", slipNormal: " << slipTpdtVertex[indexN]
+ << ", tractionNormal: " << tractionTpdtVertex[indexN]
+ << std::endl;
+#endif
+ // Don't know what behavior is appropriate so set smaller of
+ // traction and slip to zero (should be appropriate if problem
+ // is nondimensionalized correctly).
+ if (fabs(slipTpdtVertex[indexN]) > fabs(tractionTpdtVertex[indexN])) {
+ // slip is bigger, so force normal traction back to zero
+ tractionTpdtVertex[indexN] = 0.0;
+ } else {
+ // traction is bigger, so force slip back to zero
+ slipTpdtVertex[indexN] = 0.0;
+ } // if/else
+ } // if
+ if (slipTpdtVertex[indexN] < 0.0) {
+#if 0 // DEBUGGING
+ std::cout << "STEP 10 CORRECTING INTERPENETRATION"
+ << ", v_fault: " << v_fault
+ << ", slipNormal: " << slipTpdtVertex[indexN]
+ << ", tractionNormal: " << tractionTpdtVertex[indexN]
+ << std::endl;
+#endif
+ slipTpdtVertex[indexN] = 0.0;
+ } // if
+
+ // Step 2: Apply friction criterion to trial solution to get
+ // change in Lagrange multiplier (dLagrangeTpdtVertex) in fault
+ // coordinate system.
+
+ // Get friction properties and state variables.
+ _friction->retrievePropsStateVars(v_fault);
+
+ // Use fault constitutive model to compute traction associated with
+ // friction.
+ tractionMisfitVertex = 0.0;
+ const bool iterating = true; // Iterating to get friction
+ CALL_MEMBER_FN(*this,
+ constrainSolnSpaceFn)(&tractionMisfitVertex,
+ slipTpdtVertex, slipRateVertex,
+ tractionTpdtVertex,
+ iterating);
+
+ for (int iDim=0; iDim < spaceDim; ++iDim) {
+ norm2 += tractionMisfitVertex[iDim]*tractionMisfitVertex[iDim];
+ } // for
+ } // for
+
+ return sqrt(norm2) / numVertices;
+} // _constrainSolnSpaceNorm
+
+
+// ----------------------------------------------------------------------
// Constrain solution space in 1-D.
void
pylith::faults::FaultCohesiveDyn::_constrainSolnSpace1D(double_array* dTractionTpdt,
@@ -2344,320 +2632,4 @@
} // _constrainSolnSpace3D
-// ----------------------------------------------------------------------
-// Constrain solution space in 1-D.
-void
-pylith::faults::FaultCohesiveDyn::_constrainSolnSpaceImprove1D(
- double_array* dTractionTpdt,
- double_array* dSlipTpdt,
- const double_array& slipT,
- const double_array& slipTpdt,
- const double_array& tractionTpdt)
-{ // _constrainSolnSpaceImprove3D
- assert(dTractionTpdt);
- assert(dSlipTpdt);
-
- // Improving slip estimate only applies in shear. Do nothing.
-
- PetscLogFlops(0); // :TODO: Update this.
-} // _constrainSolnSpaceImprove1D
-
-
-// ----------------------------------------------------------------------
-// Constrain solution space in 3-D.
-void
-pylith::faults::FaultCohesiveDyn::_constrainSolnSpaceImprove2D(
- double_array* dTractionTpdt,
- double_array* dSlipTpdt,
- const double_array& slipT,
- const double_array& slipTpdt,
- const double_array& tractionTpdt)
-{ // _constrainSolnSpaceImprove2D
- assert(dTractionTpdt);
- assert(dSlipTpdt);
-
-#if 0 // DEBUGGING
- std::cout << "BEFORE improvement"
- << ", dTractionTpdt:";
- for (int i=0; i < 2; ++i)
- std::cout << " " << (*dTractionTpdt)[i];
- std::cout << ", dSlipTpdt:";
- for (int i=0; i < 2; ++i)
- std::cout << " " << (*dSlipTpdt)[i];
- std::cout << std::endl;
-#endif
-
- // Compute magnitude of slip and slip rate (with current increment).
- const double slipShearMag = fabs(slipTpdt[0]);
- const double slipShearNewMag = fabs(slipTpdt[0]+(*dSlipTpdt)[0]);
- const double slipRateMag = fabs(slipTpdt[0]-slipT[0]) / _dt;
- const double dSlipShearNewMag = fabs((*dSlipTpdt)[0]);
-
- const double tractionNormal = tractionTpdt[1];
-
- // Friction stress for old estimate of slip is tractionTpdt +
- // dTractionTpdt.
- const double frictionStress =
- fabs(tractionTpdt[0]+(*dTractionTpdt)[0]);
- const double tractionShearMag = fabs(tractionTpdt[0]);
-
- if (fabs(slipTpdt[1] + (*dSlipTpdt)[1]) < _zeroTolerance &&
- tractionNormal < -_zeroTolerance &&
- dSlipShearNewMag > 0.0) {
- // if in compression and no opening, and changing slip
-
- // Calculate slope (Jacobian) of friction at slip before adding
- // new increment.
- const double slopeF =
- _friction->calcFrictionSlope(slipShearMag, slipRateMag, tractionNormal);
-
-#if 0 // linear space
- const double slopeT =
- (frictionStress - tractionShearMag) / dSlipShearNewMag;
-
- // Set adjustments to increments to original increments as default.
- double dSlipShearNew2Mag = dSlipShearNewMag;
- double dTractionShearNew2Mag = frictionStress - tractionShearMag;
- if (slopeF > 0.0 && tractionShearMag > frictionStress) {
- // Strengthening, so reduce increment in slip
- dSlipShearNew2Mag =
- (tractionShearMag - frictionStress) / (slopeF-slopeT));
- dTractionShearNew2Mag += slopeF * dSlipShearNew2Mag;
- } // if
- // Ignore other cases, because slip will exceed estimate based on
- // elasticity
-
-#else // log space
-
- const double slopeT = (frictionStress - tractionShearMag) /
- (log(slipShearNewMag) - log(slipShearMag));
-
- // Set adjustments to increments to original increments as default.
- double dSlipShearNew2Mag = dSlipShearNewMag;
- double dTractionShearNew2Mag = frictionStress - tractionShearMag;
- if (slopeF > 0.0 && tractionShearMag > frictionStress) {
- // Strengthening, so reduce increment in slip
- const double slipShearMagEff = std::max(slipShearMag, _zeroTolerance);
- dSlipShearNew2Mag = slipShearMagEff *
- (-1.0 + exp((tractionShearMag - frictionStress) / (slopeF-slopeT)));
- dTractionShearNew2Mag +=
- slopeF * log((slipShearMagEff + dSlipShearNew2Mag)/slipShearMagEff);
-
- } // if
- // Ignore other cases, because slip will exceed estimate based on
- // elasticity
-#endif
-
- // Project slip and traction into vector components
- (*dSlipTpdt)[0] *= dSlipShearNew2Mag / dSlipShearNewMag;
-
- const double dTractionTpdtMag = fabs((*dTractionTpdt)[0]);
- assert(dTractionTpdtMag > 0.0);
- (*dTractionTpdt)[0] *= fabs(dTractionShearNew2Mag) / dTractionTpdtMag;
-
- // Prevent over-correction in slip resulting in backslip.
- // Expect slip direction to match tractions
-
- // Compute dot product between slip and increment in slip (want positive)
- const double slipDot =
- (slipTpdt[0] - slipT[0]) * (slipTpdt[0] + (*dSlipTpdt)[0] - slipT[0]);
- // Compute dot product of traction and slip
- const double tractionSlipDot =
- (tractionTpdt[0] + (*dTractionTpdt)[0]) * (slipTpdt[0] + (*dSlipTpdt)[0]);
- if (slipDot < 0.0 &&
- sqrt(fabs(slipDot)) > _zeroTolerance &&
- tractionSlipDot < 0.0) {
- // Correct backslip
- (*dTractionTpdt)[0] *= 0.5; // Use bisection as guess for traction
- // Use bisection for slip
-#if 0 // linear space
- (*dSlipTpdt)[0] = 0.5*(slipT[0] - slipTpdt[0]);
-#else // log space
- assert(slipT[0] * slipTpdt[0] > 0.0);
- (*dSlipTpdt)[0] *= sqrt(slipT[0] * slipTpdt[0]) / dSlipShearNew2Mag;
-#endif
- } // if/else
-
-#if 0 // DEBUGGING
- std::cout << "AFTER improvement"
- << ", dTractionTpdt:";
- for (int i=0; i < 2; ++i)
- std::cout << " " << (*dTractionTpdt)[i];
- std::cout << ", dSlipTpdt:";
- for (int i=0; i < 2; ++i)
- std::cout << " " << (*dSlipTpdt)[i];
- std::cout << ", frictionStress: " << frictionStress
- << ", tractionShearMag: " << tractionShearMag
- << ", slopeF: " << slopeF
- << ", slopeT: " << slopeT
- << std::endl;
-#endif
-
- } // if
-
- PetscLogFlops(0); // :TODO: Update this.
-
-
-} // _constrainSolnSpaceImprove2D
-
-
-// ----------------------------------------------------------------------
-// Constrain solution space in 3-D.
-void
-pylith::faults::FaultCohesiveDyn::_constrainSolnSpaceImprove3D(
- double_array* dTractionTpdt,
- double_array* dSlipTpdt,
- const double_array& slipT,
- const double_array& slipTpdt,
- const double_array& tractionTpdt)
-{ // _constrainSolnSpaceImprove3D
- assert(dTractionTpdt);
- assert(dSlipTpdt);
-
-#if 0 // DEBUGGING
- std::cout << "BEFORE improvement"
- << ", dTractionTpdt:";
- for (int i=0; i < 2; ++i)
- std::cout << " " << (*dTractionTpdt)[i];
- std::cout << ", dSlipTpdt:";
- for (int i=0; i < 2; ++i)
- std::cout << " " << (*dSlipTpdt)[i];
- std::cout << std::endl;
-#endif
-
- // Compute magnitude of slip and slip rate (with current increment).
- const double slipShearMag =
- sqrt(slipTpdt[0]*slipTpdt[0] +
- slipTpdt[1]*slipTpdt[1]);
- const double slipShearNewMag =
- sqrt(pow(slipTpdt[0]+(*dSlipTpdt)[0], 2) +
- pow(slipTpdt[1]+(*dSlipTpdt)[1], 2));
- const double slipRateMag =
- sqrt(pow(slipTpdt[0]-slipT[0], 2) +
- pow(slipTpdt[1]-slipT[1], 2)) / _dt;
- const double dSlipShearNewMag =
- sqrt(pow((*dSlipTpdt)[0], 2) +
- pow((*dSlipTpdt)[1], 2));
-
- // Friction stress for old estimate of slip is tractionTpdt +
- // dTractionTpdt.
- const double frictionStress =
- sqrt(pow(tractionTpdt[0]+(*dTractionTpdt)[0], 2) +
- pow(tractionTpdt[1]+(*dTractionTpdt)[1], 2));
- const double tractionShearMag =
- sqrt(tractionTpdt[0]*tractionTpdt[0] +
- tractionTpdt[1]*tractionTpdt[1]);
-
- const double tractionNormal = tractionTpdt[2] + (*dTractionTpdt)[2];
-
- if (fabs(slipTpdt[2] + (*dSlipTpdt)[2]) < _zeroTolerance &&
- tractionNormal < -_zeroTolerance &&
- dSlipShearNewMag > 0.0) {
- // if in compression and no opening, and changing slip
-
- // Calculate slope (Jacobian) of friction at slip before adding
- // new increment.
- const double slopeF =
- _friction->calcFrictionSlope(slipShearMag, slipRateMag, tractionNormal);
-
-#if 0 // linear space
- const double slopeT =
- (frictionStress - tractionShearMag) / dSlipShearNewMag;
-
- // Set adjustments to increments to original increments as default.
- double dSlipShearNew2Mag = dSlipShearNewMag;
- double dTractionShearNew2Mag = frictionStress - tractionShearMag;
- if (slopeF > 0.0 && tractionShearMag > frictionStress) {
- // Strengthening, so reduce increment in slip
- dSlipShearNew2Mag =
- (tractionShearMag - frictionStress) / (slopeF-slopeT));
- dTractionShearNew2Mag += slopeF * dSlipShearNew2Mag;
- } // if
- // Ignore other cases, because slip will exceed estimate based on
- // elasticity
-
-#else // log space
-
- const double slopeT = (frictionStress - tractionShearMag) /
- (log(slipShearNewMag) - log(slipShearMag));
-
- // Set adjustments to increments to original increments as default.
- double dSlipShearNew2Mag = dSlipShearNewMag;
- double dTractionShearNew2Mag = frictionStress - tractionShearMag;
- if (slopeF > 0.0 && tractionShearMag > frictionStress) {
- // Strengthening, so reduce increment in slip
- const double slipShearMagEff = std::max(slipShearMag, _zeroTolerance);
- dSlipShearNew2Mag = slipShearMagEff *
- (-1.0 + exp((tractionShearMag - frictionStress) / (slopeF-slopeT)));
- dTractionShearNew2Mag +=
- slopeF * log((slipShearMagEff + dSlipShearNew2Mag)/slipShearMagEff);
-
- } // if
- // Ignore other cases, because slip will exceed estimate based on
- // elasticity
-#endif
-
- // Project slip and traction into vector components
- // keeping same direction as original
- (*dSlipTpdt)[0] *= dSlipShearNew2Mag / dSlipShearNewMag;
- (*dSlipTpdt)[1] *= dSlipShearNew2Mag / dSlipShearNewMag;
-
- const double dTractionTpdtMag =
- sqrt(pow((*dTractionTpdt)[0], 2) +
- pow((*dTractionTpdt)[1], 2));
- assert(dTractionTpdtMag > 0.0);
- (*dTractionTpdt)[0] *= fabs(dTractionShearNew2Mag) / dTractionTpdtMag;
- (*dTractionTpdt)[1] *= fabs(dTractionShearNew2Mag) / dTractionTpdtMag;
-
- // Prevent over-correction in slip resulting in backslip.
- // Expect slip direction to match tractions
-
- // Compute dot product between slip and increment in slip (want positive)
- const double slipDot =
- (slipTpdt[0] - slipT[0]) * (slipTpdt[0] + (*dSlipTpdt)[0] - slipT[0]) +
- (slipTpdt[1] - slipT[1]) * (slipTpdt[1] + (*dSlipTpdt)[1] - slipT[1]);
- // Compute dot product of traction and slip
- const double tractionSlipDot =
- (tractionTpdt[0] + (*dTractionTpdt)[0]) * (slipTpdt[0] + (*dSlipTpdt)[0])+
- (tractionTpdt[1] + (*dTractionTpdt)[1]) * (slipTpdt[1] + (*dSlipTpdt)[1]);
- if (slipDot < 0.0 &&
- sqrt(fabs(slipDot)) > _zeroTolerance &&
- tractionSlipDot < 0.0) {
- // Correct backslip
- (*dTractionTpdt)[0] *= 0.5; // Use bisection as guess for traction
- (*dTractionTpdt)[1] *= 0.5;
- // Use bisection for slip
-#if 0 // linear space
- (*dSlipTpdt)[0] = 0.5*(slipT[0] - slipTpdt[0]);
- (*dSlipTpdt)[1] = 0.5*(slipT[1] - slipTpdt[1]);
-#else // log space
- assert(slipT[0] * slipTpdt[0] > 0.0);
- assert(slipT[1] * slipTpdt[1] > 0.0);
-
- (*dSlipTpdt)[0] *= sqrt(slipT[0] * slipTpdt[0]) / dSlipShearNew2Mag;
- (*dSlipTpdt)[1] *= sqrt(slipT[1] * slipTpdt[1]) / dSlipShearNew2Mag;
-#endif
- } // if
-
-#if 0 // DEBUGGING
- std::cout << "AFTER improvement"
- << ", dTractionTpdt:";
- for (int i=0; i < 2; ++i)
- std::cout << " " << (*dTractionTpdt)[i];
- std::cout << ", dSlipTpdt:";
- for (int i=0; i < 2; ++i)
- std::cout << " " << (*dSlipTpdt)[i];
- std::cout << ", frictionStress: " << frictionStress
- << ", tractionShearMag: " << tractionShearMag
- << ", slopeF: " << slopeF
- << ", slopeT: " << slopeT
- << std::endl;
-#endif
-
-} // if
-
- PetscLogFlops(0); // :TODO: Update this.
-} // _constrainSolnSpaceImprove3D
-
-
// End of file
Modified: short/3D/PyLith/branches/v1.6-stable/libsrc/pylith/faults/FaultCohesiveDyn.hh
===================================================================
--- short/3D/PyLith/branches/v1.6-stable/libsrc/pylith/faults/FaultCohesiveDyn.hh 2012-02-15 22:48:27 UTC (rev 19639)
+++ short/3D/PyLith/branches/v1.6-stable/libsrc/pylith/faults/FaultCohesiveDyn.hh 2012-02-16 01:36:53 UTC (rev 19640)
@@ -199,6 +199,20 @@
*/
void _sensitivityUpdateSoln(const bool negativeSide);
+ /** Compute norm of residual associated with matching fault
+ * constitutive model using update from sensitivity solve. We use
+ * this in a line search to find a good update (required because
+ * fault constitutive model may have a complex nonlinear feedback
+ * with deformation).
+ *
+ * @param alpha Current step in line search.
+ * @param fields Solution fields.
+ *
+ * @returns L2 norm of residual.
+ */
+ double _constrainSolnSpaceNorm(const double alpha,
+ topology::SolutionFields* const fields);
+
/** Constrain solution space in 1-D.
*
* @param dTractionTpdt Adjustment to fault traction vector.
@@ -241,48 +255,6 @@
const double_array& tractionTpdt,
const bool iterating =true);
- /** Improve slip estimate when constraining solution space in 1-D.
- *
- * @param dTractionTpdt Adjustment to fault traction vector.
- * @param dSlipTpdt Adjustment to fault slip vector.
- * @param slipT Fault slip vector at time t.
- * @param slipTpdt Fault slip vector at time t+dt (without adjustment).
- * @param tractionTpdt Fault traction vector (without adjustment).
- */
- void _constrainSolnSpaceImprove1D(double_array* dTractionTpdt,
- double_array* dSlipTpdt,
- const double_array& slipT,
- const double_array& slipTpdt,
- const double_array& tractionTpdt);
-
- /** Improve slip estimate when constraining solution space in 2-D.
- *
- * @param dTractionTpdt Adjustment to fault traction vector.
- * @param dSlipTpdt Adjustment to fault slip vector.
- * @param slipT Fault slip vector at time t.
- * @param slipTpdt Fault slip vector at time t+dt (without adjustment).
- * @param tractionTpdt Fault traction vector (without adjustment).
- */
- void _constrainSolnSpaceImprove2D(double_array* dTractionTpdt,
- double_array* dSlipTpdt,
- const double_array& slipT,
- const double_array& slipTpdt,
- const double_array& tractionTpdt);
-
- /** Improve slip estimate when constraining solution space in 3-D.
- *
- * @param dTractionTpdt Adjustment to fault traction vector.
- * @param dSlipTpdt Adjustment to fault slip vector.
- * @param slipT Fault slip vector at time t.
- * @param slipTpdt Fault slip vector at time t+dt (without adjustment).
- * @param tractionTpdt Fault traction vector (without adjustment).
- */
- void _constrainSolnSpaceImprove3D(double_array* dTractionTpdt,
- double_array* dSlipTpdt,
- const double_array& slipT,
- const double_array& slipTpdt,
- const double_array& tractionTpdt);
-
// PRIVATE MEMBERS ////////////////////////////////////////////////////
private :
Modified: short/3D/PyLith/branches/v1.6-stable/libsrc/pylith/friction/FrictionModel.cc
===================================================================
--- short/3D/PyLith/branches/v1.6-stable/libsrc/pylith/friction/FrictionModel.cc 2012-02-15 22:48:27 UTC (rev 19639)
+++ short/3D/PyLith/branches/v1.6-stable/libsrc/pylith/friction/FrictionModel.cc 2012-02-16 01:36:53 UTC (rev 19640)
@@ -324,28 +324,6 @@
} // calcFriction
// ----------------------------------------------------------------------
-// Compute change in friction with change in slip.
-double
-pylith::friction::FrictionModel::calcFrictionSlope(const double slip,
- const double slipRate,
- const double normalTraction)
-{ // calcFrictionSlope
- assert(_fieldsPropsStateVars);
-
- assert(_propsFiberDim+_varsFiberDim == _propsStateVarsVertex.size());
- const double* propertiesVertex = &_propsStateVarsVertex[0];
- const double* stateVarsVertex = (_varsFiberDim > 0) ?
- &_propsStateVarsVertex[_propsFiberDim] : 0;
-
- const double slope =
- _calcFrictionSlope(slip, slipRate, normalTraction,
- propertiesVertex, _propsFiberDim,
- stateVarsVertex, _varsFiberDim);
-
- return slope;
-} // calcFrictionSlope
-
-// ----------------------------------------------------------------------
// Update state variables (for next time step).
void
pylith::friction::FrictionModel::updateStateVars(const double slip,
Modified: short/3D/PyLith/branches/v1.6-stable/libsrc/pylith/friction/FrictionModel.hh
===================================================================
--- short/3D/PyLith/branches/v1.6-stable/libsrc/pylith/friction/FrictionModel.hh 2012-02-15 22:48:27 UTC (rev 19639)
+++ short/3D/PyLith/branches/v1.6-stable/libsrc/pylith/friction/FrictionModel.hh 2012-02-16 01:36:53 UTC (rev 19640)
@@ -166,21 +166,6 @@
const double slipRate,
const double normalTraction);
- /** Compute change in friction for a change in slip (Jacobian).
- *
- * @pre Must call retrievePropsAndVars for cell before calling
- * calcFriction().
- *
- * @param slip Current slip at location.
- * @param slipRate Current slip rate at location.
- * @param normalTraction Normal traction at location.
- *
- * @returns Change in friction for a chance in slip (dT/dD).
- */
- double calcFrictionSlope(const double slip,
- const double slipRate,
- const double normalTraction);
-
/** Compute friction at vertex.
*
* @pre Must call retrievePropsAndVars for cell before calling
@@ -276,27 +261,6 @@
const double* stateVars,
const int numStateVars) = 0;
- /** Compute change in friction for a change in slip (Jacobian).
- *
- * @param slip Current slip at location.
- * @param slipRate Current slip rate at location.
- * @param normalTraction Normal traction at location.
- * @param properties Properties at location.
- * @param numProperties Number of properties.
- * @param stateVars State variables at location.
- * @param numStateVars Number of state variables.
- *
- * @returns Change in friction for a chance in slip (dT/dD).
- */
- virtual
- double _calcFrictionSlope(const double slip,
- const double slipRate,
- const double normalTraction,
- const double* properties,
- const int numProperties,
- const double* stateVars,
- const int numStateVars) = 0;
-
/** Update state variables (for next time step).
*
* @param slip Current slip at location.
Modified: short/3D/PyLith/branches/v1.6-stable/libsrc/pylith/friction/RateStateAgeing.cc
===================================================================
--- short/3D/PyLith/branches/v1.6-stable/libsrc/pylith/friction/RateStateAgeing.cc 2012-02-15 22:48:27 UTC (rev 19639)
+++ short/3D/PyLith/branches/v1.6-stable/libsrc/pylith/friction/RateStateAgeing.cc 2012-02-16 01:36:53 UTC (rev 19640)
@@ -369,42 +369,6 @@
} // _calcFriction
// ----------------------------------------------------------------------
-// Compute change in friction for a change in slip (Jacobian).
-double
-pylith::friction::RateStateAgeing::_calcFrictionSlope(const double slip,
- const double slipRate,
- const double normalTraction,
- const double* properties,
- const int numProperties,
- const double* stateVars,
- const int numStateVars)
-{ // _calcFrictionSlope
- assert(properties);
- assert(_RateStateAgeing::numProperties == numProperties);
- assert(numStateVars);
- assert(_RateStateAgeing::numStateVars == numStateVars);
-
- double slope = 0.0;
- if (normalTraction <= 0.0) {
- // if fault is in compression
-
- const double a = properties[p_a];
- const double slipRate0 = properties[p_slipRate0];
-
- const double slipRateLinear = _minSlipRate;
- const double slipRateEff = std::max(slipRate, slipRateLinear);
-
- //slope = -slipRateEff / (slipRate0 * normalTraction * a);
- slope = -normalTraction * a; // log space
- } // if
-
- PetscLogFlops(5);
-
- return slope;
-} // _calcFrictionSlope
-
-
-// ----------------------------------------------------------------------
// Update state variables (for next time step).
void
pylith::friction::RateStateAgeing::_updateStateVars(const double slip,
Modified: short/3D/PyLith/branches/v1.6-stable/libsrc/pylith/friction/RateStateAgeing.hh
===================================================================
--- short/3D/PyLith/branches/v1.6-stable/libsrc/pylith/friction/RateStateAgeing.hh 2012-02-15 22:48:27 UTC (rev 19639)
+++ short/3D/PyLith/branches/v1.6-stable/libsrc/pylith/friction/RateStateAgeing.hh 2012-02-16 01:36:53 UTC (rev 19640)
@@ -141,26 +141,6 @@
const double* stateVars,
const int numStateVars);
- /** Compute change in friction for a change in slip (Jacobian).
- *
- * @param slip Current slip at location.
- * @param slipRate Current slip rate at location.
- * @param normalTraction Normal traction at location.
- * @param properties Properties at location.
- * @param numProperties Number of properties.
- * @param stateVars State variables at location.
- * @param numStateVars Number of state variables.
- *
- * @returns Change in friction for a chance in slip (dT/dD).
- */
- double _calcFrictionSlope(const double slip,
- const double slipRate,
- const double normalTraction,
- const double* properties,
- const int numProperties,
- const double* stateVars,
- const int numStateVars);
-
/** Update state variables (for next time step).
*
* @param slip Current slip at location.
Modified: short/3D/PyLith/branches/v1.6-stable/libsrc/pylith/friction/SlipWeakening.cc
===================================================================
--- short/3D/PyLith/branches/v1.6-stable/libsrc/pylith/friction/SlipWeakening.cc 2012-02-15 22:48:27 UTC (rev 19639)
+++ short/3D/PyLith/branches/v1.6-stable/libsrc/pylith/friction/SlipWeakening.cc 2012-02-16 01:36:53 UTC (rev 19640)
@@ -295,7 +295,9 @@
mu_f = properties[p_coefD];
} // if/else
friction = -mu_f * normalTraction + properties[p_cohesion];
- } // if
+ } else { // else
+ friction = properties[p_cohesion];
+ } // if/else
PetscLogFlops(10);
@@ -303,43 +305,6 @@
} // _calcFriction
// ----------------------------------------------------------------------
-// Compute change in friction for a change in slip (Jacobian).
-double
-pylith::friction::SlipWeakening::_calcFrictionSlope(const double slip,
- const double slipRate,
- const double normalTraction,
- const double* properties,
- const int numProperties,
- const double* stateVars,
- const int numStateVars)
-{ // _calcFrictionSlope
- assert(properties);
- assert(_SlipWeakening::numProperties == numProperties);
- assert(stateVars);
- assert(_SlipWeakening::numStateVars == numStateVars);
-
- double slope = 0.0;
- if (normalTraction <= 0.0) {
- // if fault is in compression
- const double slipPrev = stateVars[s_slipPrev];
- const double slipCum = stateVars[s_slipCum] + fabs(slip - slipPrev);
-
- if (slipCum < properties[p_d0]) {
- // if/else linear slip-weakening form of mu_f
- slope = -normalTraction * (properties[p_coefS] - properties[p_coefD])
- / properties[p_d0];
- } else {
- slope = 0.0;
- } // if/else
- } // if
-
- PetscLogFlops(7);
-
- return slope;
-} // _calcFrictionSlope
-
-
-// ----------------------------------------------------------------------
// Update state variables (for next time step).
void
pylith::friction::SlipWeakening::_updateStateVars(const double slip,
Modified: short/3D/PyLith/branches/v1.6-stable/libsrc/pylith/friction/SlipWeakening.hh
===================================================================
--- short/3D/PyLith/branches/v1.6-stable/libsrc/pylith/friction/SlipWeakening.hh 2012-02-15 22:48:27 UTC (rev 19639)
+++ short/3D/PyLith/branches/v1.6-stable/libsrc/pylith/friction/SlipWeakening.hh 2012-02-16 01:36:53 UTC (rev 19640)
@@ -129,26 +129,6 @@
const double* stateVars,
const int numStateVars);
- /** Compute change in friction for a change in slip (Jacobian).
- *
- * @param slip Current slip at location.
- * @param slipRate Current slip rate at location.
- * @param normalTraction Normal traction at location.
- * @param properties Properties at location.
- * @param numProperties Number of properties.
- * @param stateVars State variables at location.
- * @param numStateVars Number of state variables.
- *
- * @returns Change in friction for a chance in slip (dT/dD).
- */
- double _calcFrictionSlope(const double slip,
- const double slipRate,
- const double normalTraction,
- const double* properties,
- const int numProperties,
- const double* stateVars,
- const int numStateVars);
-
/** Update state variables (for next time step).
*
* @param slip Current slip at location.
Modified: short/3D/PyLith/branches/v1.6-stable/libsrc/pylith/friction/StaticFriction.cc
===================================================================
--- short/3D/PyLith/branches/v1.6-stable/libsrc/pylith/friction/StaticFriction.cc 2012-02-15 22:48:27 UTC (rev 19639)
+++ short/3D/PyLith/branches/v1.6-stable/libsrc/pylith/friction/StaticFriction.cc 2012-02-16 01:36:53 UTC (rev 19640)
@@ -166,21 +166,4 @@
} // _calcFriction
-// ----------------------------------------------------------------------
-// Compute change in friction for a change in slip (Jacobian).
-double
-pylith::friction::StaticFriction::_calcFrictionSlope(const double slip,
- const double slipRate,
- const double normalTraction,
- const double* properties,
- const int numProperties,
- const double* stateVars,
- const int numStateVars)
-{ // _calcFrictionSlope
- const double slope = 0;
-
- return slope;
-} // _calcFrictionSlope
-
-
// End of file
Modified: short/3D/PyLith/branches/v1.6-stable/libsrc/pylith/friction/StaticFriction.hh
===================================================================
--- short/3D/PyLith/branches/v1.6-stable/libsrc/pylith/friction/StaticFriction.hh 2012-02-15 22:48:27 UTC (rev 19639)
+++ short/3D/PyLith/branches/v1.6-stable/libsrc/pylith/friction/StaticFriction.hh 2012-02-16 01:36:53 UTC (rev 19640)
@@ -94,26 +94,6 @@
const double* stateVars,
const int numStateVars);
- /** Compute change in friction for a change in slip (Jacobian).
- *
- * @param slip Current slip at location.
- * @param slipRate Current slip rate at location.
- * @param normalTraction Normal traction at location.
- * @param properties Properties at location.
- * @param numProperties Number of properties.
- * @param stateVars State variables at location.
- * @param numStateVars Number of state variables.
- *
- * @returns Change in friction for a chance in slip (dT/dD).
- */
- double _calcFrictionSlope(const double slip,
- const double slipRate,
- const double normalTraction,
- const double* properties,
- const int numProperties,
- const double* stateVars,
- const int numStateVars);
-
// PRIVATE MEMBERS ////////////////////////////////////////////////////
private :
Modified: short/3D/PyLith/branches/v1.6-stable/libsrc/pylith/friction/TimeWeakening.cc
===================================================================
--- short/3D/PyLith/branches/v1.6-stable/libsrc/pylith/friction/TimeWeakening.cc 2012-02-15 22:48:27 UTC (rev 19639)
+++ short/3D/PyLith/branches/v1.6-stable/libsrc/pylith/friction/TimeWeakening.cc 2012-02-16 01:36:53 UTC (rev 19640)
@@ -308,21 +308,4 @@
} // _updateStateVars
-// ----------------------------------------------------------------------
-// Compute change in friction for a change in slip (Jacobian).
-double
-pylith::friction::TimeWeakening::_calcFrictionSlope(const double slip,
- const double slipRate,
- const double normalTraction,
- const double* properties,
- const int numProperties,
- const double* stateVars,
- const int numStateVars)
-{ // _calcFrictionSlope
- const double slope = 0;
-
- return slope;
-} // _calcFrictionSlope
-
-
// End of file
Modified: short/3D/PyLith/branches/v1.6-stable/libsrc/pylith/friction/TimeWeakening.hh
===================================================================
--- short/3D/PyLith/branches/v1.6-stable/libsrc/pylith/friction/TimeWeakening.hh 2012-02-15 22:48:27 UTC (rev 19639)
+++ short/3D/PyLith/branches/v1.6-stable/libsrc/pylith/friction/TimeWeakening.hh 2012-02-16 01:36:53 UTC (rev 19640)
@@ -123,26 +123,6 @@
const double* stateVars,
const int numStateVars);
- /** Compute change in friction for a change in slip (Jacobian).
- *
- * @param slip Current slip at location.
- * @param slipRate Current slip rate at location.
- * @param normalTraction Normal traction at location.
- * @param properties Properties at location.
- * @param numProperties Number of properties.
- * @param stateVars State variables at location.
- * @param numStateVars Number of state variables.
- *
- * @returns Change in friction for a chance in slip (dT/dD).
- */
- double _calcFrictionSlope(const double slip,
- const double slipRate,
- const double normalTraction,
- const double* properties,
- const int numProperties,
- const double* stateVars,
- const int numStateVars);
-
/** Update state variables (for next time step).
*
* @param slip Current slip at location.
Modified: short/3D/PyLith/branches/v1.6-stable/libsrc/pylith/problems/SolverNonlinear.cc
===================================================================
--- short/3D/PyLith/branches/v1.6-stable/libsrc/pylith/problems/SolverNonlinear.cc 2012-02-15 22:48:27 UTC (rev 19639)
+++ short/3D/PyLith/branches/v1.6-stable/libsrc/pylith/problems/SolverNonlinear.cc 2012-02-16 01:36:53 UTC (rev 19640)
@@ -23,6 +23,7 @@
#include "Formulation.hh" // USES Formulation
#include "pylith/topology/SolutionFields.hh" // USES SolutionFields
#include "pylith/topology/Jacobian.hh" // USES Jacobian
+#include "pylith/utils/constdefs.h" // USES PYLITH_MAXDOUBLE
#include "pylith/utils/EventLogger.hh" // USES EventLogger
#include <petscsnes.h> // USES PetscSNES
@@ -243,7 +244,10 @@
*ynorm = snes->maxstep;
}
ierr = VecMaxPointwiseDivide(y,x,&rellength);CHKERRQ(ierr);
- minlambda = snes->steptol/rellength;
+
+ // Place reasonable absolute limit on minimum lambda
+ minlambda = std::max(snes->steptol/rellength, 1.0/PYLITH_MAXDOUBLE);
+
ierr = MatMult(snes->jacobian,y,w);CHKERRQ(ierr);
#if defined(PETSC_USE_COMPLEX)
ierr = VecDot(f,w,&cinitslope);CHKERRQ(ierr);
Modified: short/3D/PyLith/branches/v1.6-stable/libsrc/pylith/utils/constdefs.h
===================================================================
--- short/3D/PyLith/branches/v1.6-stable/libsrc/pylith/utils/constdefs.h 2012-02-15 22:48:27 UTC (rev 19639)
+++ short/3D/PyLith/branches/v1.6-stable/libsrc/pylith/utils/constdefs.h 2012-02-16 01:36:53 UTC (rev 19640)
@@ -26,7 +26,7 @@
#define pylith_utils_constdefs_h
namespace pylith {
- static const double PYLITH_MAXDOUBLE = 1.0e+30;
+ static const double PYLITH_MAXDOUBLE = 1.0e+99;
static const float PYLITH_MAXFLOAT = 1.0e+30;
}
Modified: short/3D/PyLith/branches/v1.6-stable/modulesrc/friction/FrictionModel.i
===================================================================
--- short/3D/PyLith/branches/v1.6-stable/modulesrc/friction/FrictionModel.i 2012-02-15 22:48:27 UTC (rev 19639)
+++ short/3D/PyLith/branches/v1.6-stable/modulesrc/friction/FrictionModel.i 2012-02-16 01:36:53 UTC (rev 19640)
@@ -147,21 +147,6 @@
const double normalTraction);
- /** Compute change in friction for a change in slip (Jacobian).
- *
- * @pre Must call retrievePropsAndVars for cell before calling
- * calcFriction().
- *
- * @param slip Current slip at location.
- * @param slipRate Current slip rate at location.
- * @param normalTraction Normal traction at location.
- *
- * @returns Change in friction for a chance in slip (dT/dD).
- */
- double calcFrictionSlope(const double slip,
- const double slipRate,
- const double normalTraction);
-
/** Compute friction at vertex.
*
* @pre Must call retrievePropsAndVars for cell before calling
@@ -255,27 +240,6 @@
const double* stateVars,
const int numStateVars) = 0;
- /** Compute change in friction for a change in slip (Jacobian).
- *
- * @param slip Current slip at location.
- * @param slipRate Current slip rate at location.
- * @param normalTraction Normal traction at location.
- * @param properties Properties at location.
- * @param numProperties Number of properties.
- * @param stateVars State variables at location.
- * @param numStateVars Number of state variables.
- *
- * @returns Change in friction for a chance in slip (dT/dD).
- */
- virtual
- double _calcFrictionSlope(const double slip,
- const double slipRate,
- const double normalTraction,
- const double* properties,
- const int numProperties,
- const double* stateVars,
- const int numStateVars) = 0;
-
/** Update state variables (for next time step).
*
* @param stateVars State variables at location.
Modified: short/3D/PyLith/branches/v1.6-stable/modulesrc/friction/RateStateAgeing.i
===================================================================
--- short/3D/PyLith/branches/v1.6-stable/modulesrc/friction/RateStateAgeing.i 2012-02-15 22:48:27 UTC (rev 19639)
+++ short/3D/PyLith/branches/v1.6-stable/modulesrc/friction/RateStateAgeing.i 2012-02-16 01:36:53 UTC (rev 19640)
@@ -88,26 +88,6 @@
const double* stateVars,
const int numStateVars);
- /** Compute change in friction for a change in slip (Jacobian).
- *
- * @param slip Current slip at location.
- * @param slipRate Current slip rate at location.
- * @param normalTraction Normal traction at location.
- * @param properties Properties at location.
- * @param numProperties Number of properties.
- * @param stateVars State variables at location.
- * @param numStateVars Number of state variables.
- *
- * @returns Change in friction for a chance in slip (dT/dD).
- */
- double _calcFrictionSlope(const double slip,
- const double slipRate,
- const double normalTraction,
- const double* properties,
- const int numProperties,
- const double* stateVars,
- const int numStateVars);
-
/** Update state variables (for next time step).
*
* @param slip Current slip at location.
Modified: short/3D/PyLith/branches/v1.6-stable/modulesrc/friction/SlipWeakening.i
===================================================================
--- short/3D/PyLith/branches/v1.6-stable/modulesrc/friction/SlipWeakening.i 2012-02-15 22:48:27 UTC (rev 19639)
+++ short/3D/PyLith/branches/v1.6-stable/modulesrc/friction/SlipWeakening.i 2012-02-16 01:36:53 UTC (rev 19640)
@@ -87,26 +87,6 @@
const double* stateVars,
const int numStateVars);
- /** Compute change in friction for a change in slip (Jacobian).
- *
- * @param slip Current slip at location.
- * @param slipRate Current slip rate at location.
- * @param normalTraction Normal traction at location.
- * @param properties Properties at location.
- * @param numProperties Number of properties.
- * @param stateVars State variables at location.
- * @param numStateVars Number of state variables.
- *
- * @returns Change in friction for a chance in slip (dT/dD).
- */
- double _calcFrictionSlope(const double slip,
- const double slipRate,
- const double normalTraction,
- const double* properties,
- const int numProperties,
- const double* stateVars,
- const int numStateVars);
-
/** Update state variables (for next time step).
*
* @param slip Current slip at location.
Modified: short/3D/PyLith/branches/v1.6-stable/modulesrc/friction/StaticFriction.i
===================================================================
--- short/3D/PyLith/branches/v1.6-stable/modulesrc/friction/StaticFriction.i 2012-02-15 22:48:27 UTC (rev 19639)
+++ short/3D/PyLith/branches/v1.6-stable/modulesrc/friction/StaticFriction.i 2012-02-16 01:36:53 UTC (rev 19640)
@@ -81,26 +81,6 @@
const double* stateVars,
const int numStateVars);
- /** Compute change in friction for a change in slip (Jacobian).
- *
- * @param slip Current slip at location.
- * @param slipRate Current slip rate at location.
- * @param normalTraction Normal traction at location.
- * @param properties Properties at location.
- * @param numProperties Number of properties.
- * @param stateVars State variables at location.
- * @param numStateVars Number of state variables.
- *
- * @returns Change in friction for a chance in slip (dT/dD).
- */
- double _calcFrictionSlope(const double slip,
- const double slipRate,
- const double normalTraction,
- const double* properties,
- const int numProperties,
- const double* stateVars,
- const int numStateVars);
-
}; // class StaticFriction
} // friction
Modified: short/3D/PyLith/branches/v1.6-stable/modulesrc/friction/TimeWeakening.i
===================================================================
--- short/3D/PyLith/branches/v1.6-stable/modulesrc/friction/TimeWeakening.i 2012-02-15 22:48:27 UTC (rev 19639)
+++ short/3D/PyLith/branches/v1.6-stable/modulesrc/friction/TimeWeakening.i 2012-02-16 01:36:53 UTC (rev 19640)
@@ -81,26 +81,6 @@
const double* stateVars,
const int numStateVars);
- /** Compute change in friction for a change in slip (Jacobian).
- *
- * @param slip Current slip at location.
- * @param slipRate Current slip rate at location.
- * @param normalTraction Normal traction at location.
- * @param properties Properties at location.
- * @param numProperties Number of properties.
- * @param stateVars State variables at location.
- * @param numStateVars Number of state variables.
- *
- * @returns Change in friction for a chance in slip (dT/dD).
- */
- double _calcFrictionSlope(const double slip,
- const double slipRate,
- const double normalTraction,
- const double* properties,
- const int numProperties,
- const double* stateVars,
- const int numStateVars);
-
/** Update state variables (for next time step).
*
* @param slip Current slip at location.
Modified: short/3D/PyLith/branches/v1.6-stable/tests/2d/frictionslide/plot_friction.py
===================================================================
--- short/3D/PyLith/branches/v1.6-stable/tests/2d/frictionslide/plot_friction.py 2012-02-15 22:48:27 UTC (rev 19639)
+++ short/3D/PyLith/branches/v1.6-stable/tests/2d/frictionslide/plot_friction.py 2012-02-16 01:36:53 UTC (rev 19640)
@@ -1,6 +1,6 @@
#!/usr/bin/env python
-sim = "ratestate_weak"
+sim = "ratestate_stable"
# ======================================================================
import tables
Modified: short/3D/PyLith/branches/v1.6-stable/tests/2d/frictionslide/pylithapp.cfg
===================================================================
--- short/3D/PyLith/branches/v1.6-stable/tests/2d/frictionslide/pylithapp.cfg 2012-02-15 22:48:27 UTC (rev 19639)
+++ short/3D/PyLith/branches/v1.6-stable/tests/2d/frictionslide/pylithapp.cfg 2012-02-16 01:36:53 UTC (rev 19640)
@@ -78,7 +78,7 @@
[pylithapp.timedependent.interfaces]
fault = pylith.faults.FaultCohesiveDyn
-fault.zero_tolerance = 1.0e-12
+fault.zero_tolerance = 1.0e-13
[pylithapp.timedependent.interfaces.fault]
id = 100
Modified: short/3D/PyLith/branches/v1.6-stable/tests/3d/cyclicfriction/pylithapp.cfg
===================================================================
--- short/3D/PyLith/branches/v1.6-stable/tests/3d/cyclicfriction/pylithapp.cfg 2012-02-15 22:48:27 UTC (rev 19639)
+++ short/3D/PyLith/branches/v1.6-stable/tests/3d/cyclicfriction/pylithapp.cfg 2012-02-16 01:36:53 UTC (rev 19640)
@@ -114,7 +114,7 @@
[pylithapp.timedependent.interfaces.fault]
label = fault
-zero_tolerance = 1.0e-10
+zero_tolerance = 1.0e-12
friction = pylith.friction.SlipWeakening
friction.label = Slip weakening
@@ -171,7 +171,7 @@
sub_pc_factor_shift_type = nonzero
# Convergence parameters.
-ksp_rtol = 1.0e-10
+ksp_rtol = 1.0e-12
ksp_atol = 1.0e-12
ksp_max_it = 100
ksp_gmres_restart = 50
@@ -182,11 +182,12 @@
ksp_converged_reason = true
# Nonlinear solver monitoring options.
-snes_rtol = 1.0e-10
-snes_atol = 1.0e-9
+snes_rtol = 1.0e-12
+snes_atol = 1.0e-10
snes_max_it = 100
snes_monitor = true
snes_ls_monitor = true
+#snes_steptol = 1.0e-20
#snes_view = true
snes_converged_reason = true
snes_monitor_solution_update = true
More information about the CIG-COMMITS
mailing list