[cig-commits] r19645 - in short/3D/PyLith/branches/v1.6-stable: libsrc/pylith/faults libsrc/pylith/problems tests/3d/cyclicfriction
brad at geodynamics.org
brad at geodynamics.org
Fri Feb 17 17:19:16 PST 2012
Author: brad
Date: 2012-02-17 17:19:15 -0800 (Fri, 17 Feb 2012)
New Revision: 19645
Modified:
short/3D/PyLith/branches/v1.6-stable/libsrc/pylith/faults/FaultCohesiveDyn.cc
short/3D/PyLith/branches/v1.6-stable/libsrc/pylith/problems/SolverNonlinear.cc
short/3D/PyLith/branches/v1.6-stable/tests/3d/cyclicfriction/pylithapp.cfg
Log:
Only apply line search in constraining solution space if fault does not open.
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-17 01:53:04 UTC (rev 19644)
+++ short/3D/PyLith/branches/v1.6-stable/libsrc/pylith/faults/FaultCohesiveDyn.cc 2012-02-18 01:19:15 UTC (rev 19645)
@@ -813,6 +813,7 @@
// if residual is very small, we prefer the full step
break;
+#if 0
std::cout << "alphaL: " << pow(10.0, logAlphaL)
<< ", residuaL: " << residualL
<< ", alphaM: " << pow(10.0, logAlphaM)
@@ -820,26 +821,27 @@
<< ", alphaR: " << pow(10.0, logAlphaR)
<< ", residualR: " << residualR
<< std::endl;
+#endif
- if (residualL <= residualML && residualL <= residualM && residualL <= residualMR && residualL <= residualR) {
+ if (residualL < residualML && residualL < residualM && residualL < residualMR && residualL < residualR) {
logAlphaL = logAlphaL;
logAlphaR = logAlphaM;
residualL = residualL;
residualR = residualM;
residualM = residualML;
- } else if (residualML <= residualL && residualML <= residualM && residualML <= residualMR && residualML <= residualR) {
+ } else if (residualML <= residualL && residualML < residualM && residualML < residualMR && residualML < residualR) {
logAlphaL = logAlphaL;
logAlphaR = logAlphaM;
residualL = residualL;
residualR = residualM;
residualM = residualML;
- } else if (residualM <= residualL && residualM <= residualML && residualM <= residualMR && residualM <= residualR) {
+ } else if (residualM <= residualL && residualM <= residualML && residualM < residualMR && residualM < residualR) {
logAlphaL = logAlphaML;
logAlphaR = logAlphaMR;
residualL = residualML;
residualR = residualMR;
residualM = residualM;
- } else if (residualMR <= residualL && residualMR <= residualML && residualMR <= residualM && residualMR <= residualR) {
+ } else if (residualMR <= residualL && residualMR <= residualML && residualMR <= residualM && residualMR < residualR) {
logAlphaL = logAlphaM;
logAlphaR = logAlphaR;
residualL = residualM;
@@ -1014,6 +1016,7 @@
} // if/else
+#if 0 // UNNECESSARY????
// Step 5d: Prevent over-correction in shear slip resulting in backslip
double slipDot = 0.0;
double tractionSlipDot = 0.0;
@@ -1030,7 +1033,7 @@
sqrt(fabs(slipDot)) > _zeroTolerance &&
tractionSlipDot < 0.0) {
#if 1 // DEBUGGING
- std::cout << "STEP 5c CORRECTING BACKSLIP"
+ std::cout << "STEP 5d CORRECTING BACKSLIP"
<< ", v_fault: " << v_fault
<< ", slipDot: " << slipDot
<< ", tractionSlipDot: " << tractionSlipDot
@@ -1042,6 +1045,7 @@
dSlipTpdtVertex[iDim] = 0.5*(slipTVertex[iDim] - slipTpdtVertex[iDim]);
} // for
} // if/else
+#endif
// Update current estimate of slip from t to t+dt.
slipTpdtVertex += dSlipTpdtVertex;
@@ -2339,6 +2343,7 @@
fields->get("dispIncr(t->t+dt)").section();
assert(!dispIncrSection.isNull());
+ bool isOpening = false;
double norm2 = 0.0;
const int numVertices = _cohesiveVertices.size();
for (int iVertex=0; iVertex < numVertices; ++iVertex) {
@@ -2452,6 +2457,10 @@
slipTpdtVertex[indexN] = 0.0;
} // if
+ if (slipTpdtVertex[indexN] > _zeroTolerance) {
+ isOpening = true;
+ } // if
+
// Apply friction criterion to trial solution to get change in
// Lagrange multiplier (dLagrangeTpdtVertex) in fault coordinate
// system.
@@ -2496,6 +2505,10 @@
} // for
} // for
+ if (isOpening && alpha < 1.0) {
+ norm2 = PYLITH_MAXDOUBLE;
+ } // if
+
return sqrt(norm2) / numVertices;
} // _constrainSolnSpaceNorm
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-17 01:53:04 UTC (rev 19644)
+++ short/3D/PyLith/branches/v1.6-stable/libsrc/pylith/problems/SolverNonlinear.cc 2012-02-18 01:19:15 UTC (rev 19645)
@@ -329,10 +329,12 @@
ierr = PetscViewerASCIIPrintf(snes->ls_monitor," Line search: fnorm=%18.16e, gnorm=%18.16e, ynorm=%18.16e, minlambda=%18.16e, lambda=%18.16e, initial slope=%18.16e\n",(double)fnorm,(double)*gnorm,(double)*ynorm,(double)minlambda,(double)lambda,(double)initslope);CHKERRQ(ierr);
ierr = PetscViewerASCIISubtractTab(snes->ls_monitor,((PetscObject)snes)->tablevel);CHKERRQ(ierr);
}
+#if 0 // DEBUGGING
assert(lsctx);
Formulation* formulation = (Formulation*) lsctx;
assert(formulation);
formulation->printState(&w, &g, &x, &y);
+#endif
#if 0 // Original PETSc code
*flag = PETSC_FALSE; // DIVERGED_LINE_SEARCH
#else
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-17 01:53:04 UTC (rev 19644)
+++ short/3D/PyLith/branches/v1.6-stable/tests/3d/cyclicfriction/pylithapp.cfg 2012-02-18 01:19:15 UTC (rev 19645)
@@ -33,7 +33,7 @@
# problem
# ----------------------------------------------------------------------
[pylithapp.timedependent.formulation.time_step]
-total_time = 60.0*hour ; total time of simulation
+total_time = 90.0*hour ; total time of simulation
dt = 1.0*hour
[pylithapp.timedependent.normalizer]
@@ -74,7 +74,7 @@
db_change = spatialdata.spatialdb.UniformDB
db_change.label = Amplitude of Dirichlet BC on +x
db_change.values = [displacement-x,displacement-y,displacement-z,change-start-time]
-db_change.data = [-0.5*m,-0.2*m,0.0*m,0.0*hour]
+db_change.data = [-0.5*m,-0.3*m,0.0*m,0.0*hour]
th_change = spatialdata.spatialdb.TimeHistory
th_change.label = Time history of tidal load
@@ -89,7 +89,7 @@
db_change = spatialdata.spatialdb.UniformDB
db_change.label = Amplitude of Dirichlet BC on -x
db_change.values = [displacement-x,displacement-y,displacement-z,change-start-time]
-db_change.data = [+0.5*m,+0.2*m,0.0*m,0.0*hour]
+db_change.data = [+0.5*m,+0.3*m,0.0*m,0.0*hour]
th_change = spatialdata.spatialdb.TimeHistory
th_change.label = Time history of tidal load
More information about the CIG-COMMITS
mailing list