[cig-commits] r19154 - in short/3D/PyLith/trunk: libsrc/pylith/faults libsrc/pylith/problems tests/2d/frictionslide tests/3d/cyclicfriction
brad at geodynamics.org
brad at geodynamics.org
Mon Nov 7 10:43:31 PST 2011
Author: brad
Date: 2011-11-07 10:43:31 -0800 (Mon, 07 Nov 2011)
New Revision: 19154
Added:
short/3D/PyLith/trunk/tests/3d/cyclicfriction/fieldsplit.cfg
Modified:
short/3D/PyLith/trunk/libsrc/pylith/faults/FaultCohesiveDyn.cc
short/3D/PyLith/trunk/libsrc/pylith/faults/FaultCohesiveDyn.hh
short/3D/PyLith/trunk/libsrc/pylith/problems/Formulation.cc
short/3D/PyLith/trunk/libsrc/pylith/problems/Formulation.hh
short/3D/PyLith/trunk/libsrc/pylith/problems/SolverNonlinear.cc
short/3D/PyLith/trunk/tests/2d/frictionslide/plot_friction.py
short/3D/PyLith/trunk/tests/2d/frictionslide/pylithapp.cfg
short/3D/PyLith/trunk/tests/2d/frictionslide/ratestate.cfg
short/3D/PyLith/trunk/tests/3d/cyclicfriction/geometry.jou
short/3D/PyLith/trunk/tests/3d/cyclicfriction/hex8_1000m.exo
short/3D/PyLith/trunk/tests/3d/cyclicfriction/hex8_1000m.jou
short/3D/PyLith/trunk/tests/3d/cyclicfriction/initial_tractions.spatialdb
short/3D/PyLith/trunk/tests/3d/cyclicfriction/pylithapp.cfg
Log:
Merge from stable.
Modified: short/3D/PyLith/trunk/libsrc/pylith/faults/FaultCohesiveDyn.cc
===================================================================
--- short/3D/PyLith/trunk/libsrc/pylith/faults/FaultCohesiveDyn.cc 2011-11-06 22:16:51 UTC (rev 19153)
+++ short/3D/PyLith/trunk/libsrc/pylith/faults/FaultCohesiveDyn.cc 2011-11-07 18:43:31 UTC (rev 19154)
@@ -173,15 +173,15 @@
assert(0 != _fields);
assert(0 != _logger);
- // Initial fault tractions have been assembled, so they do not need
- // assembling across processors.
+ // Cohesive cells with conventional vertices N and P, and constraint
+ // vertex L make contributions to the assembled residual:
+ //
+ // DOF P: \int_{S_f^+} \tensor{N}_m^T \cdot \tensor{N}_p \cdot \vec{l}_p dS
+ // DOF N: -\int_{S_f^+} \tensor{N}_m^T \cdot \tensor{N}_p \cdot \vec{l}_p dS
+ // DOF L: \int_S_f \tensor{N}_p^T ( \tensor{R} \cdot \vec{d}
+ // -\tensor{N}_{n^+} \cdot \vec{u}_{n^+}
+ // +\tensor{N}_{n^-} \cdot \vec{u}_{n^-} dS
- FaultCohesiveLagrange::integrateResidual(residual, t, fields);
-
- // No contribution if no initial tractions are specified.
- if (0 == _dbInitialTract)
- return;
-
const int setupEvent = _logger->eventId("FaIR setup");
const int geometryEvent = _logger->eventId("FaIR geometry");
const int computeEvent = _logger->eventId("FaIR compute");
@@ -196,20 +196,31 @@
// Get sections associated with cohesive cells
scalar_array residualVertexN(spaceDim);
scalar_array residualVertexP(spaceDim);
+ scalar_array residualVertexL(spaceDim);
const ALE::Obj<RealSection>& residualSection = residual.section();
assert(!residualSection.isNull());
+ const ALE::Obj<RealSection>& dispTSection = fields->get("disp(t)").section();
+ assert(!dispTSection.isNull());
+
+ const ALE::Obj<RealSection>& dispTIncrSection =
+ fields->get("dispIncr(t->t+dt)").section();
+ assert(!dispTIncrSection.isNull());
+
+ scalar_array dispTpdtVertexN(spaceDim);
+ scalar_array dispTpdtVertexP(spaceDim);
+ scalar_array dispTpdtVertexL(spaceDim);
+
+ scalar_array initialTractionsVertex(spaceDim);
+ ALE::Obj<RealSection> initialTractionsSection;
+ if (_dbInitialTract) {
+ initialTractionsSection = _fields->get("initial traction").section();
+ assert(!initialTractionsSection.isNull());
+ } // if
+
const ALE::Obj<RealSection>& areaSection = _fields->get("area").section();
assert(!areaSection.isNull());
- const ALE::Obj<RealSection>& initialTractionsSection =
- _fields->get("initial traction").section();
- assert(!initialTractionsSection.isNull());
-
- const ALE::Obj<RealSection>& dispRelSection =
- _fields->get("relative disp").section();
- assert(!dispRelSection.isNull());
-
const ALE::Obj<RealSection>& orientationSection =
_fields->get("orientation").section();
assert(!orientationSection.isNull());
@@ -244,53 +255,89 @@
#endif
// Get initial tractions at fault vertex.
- assert(spaceDim == initialTractionsSection->getFiberDimension(v_fault));
- const PylithScalar* initialTractionsVertex =
- initialTractionsSection->restrictPoint(v_fault);
- assert(initialTractionsVertex);
+ if (_dbInitialTract) {
+ initialTractionsSection->restrictPoint(v_fault,
+ &initialTractionsVertex[0],
+ initialTractionsVertex.size());
+ } else {
+ initialTractionsVertex = 0.0;
+ } // if/else
- // Get relative dislplacement at fault vertex.
- assert(spaceDim == dispRelSection->getFiberDimension(v_fault));
- const PylithScalar* dispRelVertex = dispRelSection->restrictPoint(v_fault);
- assert(dispRelVertex);
+ // Get orientation associated with fault vertex.
+ assert(spaceDim*spaceDim == orientationSection->getFiberDimension(v_fault));
+ const PylithScalar* orientationVertex =
+ orientationSection->restrictPoint(v_fault);
+ assert(orientationVertex);
// Get area associated with fault vertex.
assert(1 == areaSection->getFiberDimension(v_fault));
assert(areaSection->restrictPoint(v_fault));
const PylithScalar areaVertex = *areaSection->restrictPoint(v_fault);
- // Get orientation associated with fault vertex.
- assert(spaceDim*spaceDim == orientationSection->getFiberDimension(v_fault));
- const PylithScalar* orientationVertex =
- orientationSection->restrictPoint(v_fault);
- assert(orientationVertex);
+ // Get disp(t) at conventional vertices and Lagrange vertex.
+ assert(spaceDim == dispTSection->getFiberDimension(v_negative));
+ const PylithScalar* dispTVertexN = dispTSection->restrictPoint(v_negative);
+ assert(dispTVertexN);
+ assert(spaceDim == dispTSection->getFiberDimension(v_positive));
+ const PylithScalar* dispTVertexP = dispTSection->restrictPoint(v_positive);
+ assert(dispTVertexP);
+
+ assert(spaceDim == dispTSection->getFiberDimension(v_lagrange));
+ const PylithScalar* dispTVertexL = dispTSection->restrictPoint(v_lagrange);
+ assert(dispTVertexL);
+
+ // Get dispIncr(t->t+dt) at conventional vertices and Lagrange vertex.
+ assert(spaceDim == dispTIncrSection->getFiberDimension(v_negative));
+ const PylithScalar* dispTIncrVertexN =
+ dispTIncrSection->restrictPoint(v_negative);
+ assert(dispTIncrVertexN);
+
+ assert(spaceDim == dispTIncrSection->getFiberDimension(v_positive));
+ const PylithScalar* dispTIncrVertexP =
+ dispTIncrSection->restrictPoint(v_positive);
+ assert(dispTIncrVertexP);
+
+ assert(spaceDim == dispTIncrSection->getFiberDimension(v_lagrange));
+ const PylithScalar* dispTIncrVertexL =
+ dispTIncrSection->restrictPoint(v_lagrange);
+ assert(dispTIncrVertexL);
+
#if defined(DETAILED_EVENT_LOGGING)
_logger->eventEnd(restrictEvent);
_logger->eventBegin(computeEvent);
#endif
- // Initial (external) tractions oppose (internal) tractions
- // associated with Lagrange multiplier, so these terms have the
- // opposite sign as the integration of the Lagrange multipliers in
- // FaultCohesiveLagrange.
+ // Compute current estimate of displacement at time t+dt using
+ // solution increment.
for (int iDim=0; iDim < spaceDim; ++iDim) {
- residualVertexP[iDim] = areaVertex * initialTractionsVertex[iDim];
+ dispTpdtVertexN[iDim] = dispTVertexN[iDim] + dispTIncrVertexN[iDim];
+ dispTpdtVertexP[iDim] = dispTVertexP[iDim] + dispTIncrVertexP[iDim];
+ dispTpdtVertexL[iDim] = dispTVertexL[iDim] + dispTIncrVertexL[iDim];
} // for
- residualVertexN = -residualVertexP;
-
- // Only apply initial tractions if there is no opening.
- // If there is opening, zero out initial tractions
+
+ // Compute slip (in fault coordinates system) from displacements.
PylithScalar slipNormal = 0.0;
+ PylithScalar tractionNormal = 0.0;
const int indexN = spaceDim - 1;
for (int jDim=0; jDim < spaceDim; ++jDim) {
- slipNormal += orientationVertex[indexN*spaceDim+jDim]*dispRelVertex[jDim];
+ slipNormal += orientationVertex[indexN*spaceDim+jDim] *
+ (dispTpdtVertexP[jDim] - dispTpdtVertexN[jDim]);
+ tractionNormal +=
+ orientationVertex[indexN*spaceDim+jDim] * dispTpdtVertexL[jDim];
} // for
+
+ residualVertexN = 0.0;
+ residualVertexL = 0.0;
+ if (slipNormal < _zeroTolerance) { // if no opening
+ // Initial (external) tractions oppose (internal) tractions
+ // associated with Lagrange multiplier.
+ residualVertexN = areaVertex * (dispTpdtVertexL - initialTractionsVertex);
- if (slipNormal > 0.0) {
- residualVertexN = 0.0;
- residualVertexP = 0.0;
- } // if
+ } else { // opening, normal traction should be zero
+ assert(fabs(tractionNormal) < _zeroTolerance);
+ } // if/else
+ residualVertexP = -residualVertexN;
#if defined(DETAILED_EVENT_LOGGING)
_logger->eventEnd(computeEvent);
@@ -306,11 +353,15 @@
residualSection->getFiberDimension(v_positive));
residualSection->updateAddPoint(v_positive, &residualVertexP[0]);
+ assert(residualVertexL.size() ==
+ residualSection->getFiberDimension(v_lagrange));
+ residualSection->updateAddPoint(v_lagrange, &residualVertexL[0]);
+
#if defined(DETAILED_EVENT_LOGGING)
_logger->eventEnd(updateEvent);
#endif
} // for
- PetscLogFlops(numVertices*spaceDim*(2+spaceDim*2));
+ PetscLogFlops(numVertices*spaceDim*8);
#if !defined(DETAILED_EVENT_LOGGING)
_logger->eventEnd(computeEvent);
@@ -327,7 +378,7 @@
assert(0 != fields);
assert(0 != _fields);
- _updateVelRel(*fields);
+ _updateRelMotion(*fields);
const int spaceDim = _quadrature->spaceDim();
@@ -460,16 +511,19 @@
assert(0 != _fields);
assert(0 != _friction);
- _updateVelRel(*fields);
_sensitivitySetup(jacobian);
// Update time step in friction (can vary).
_friction->timeStep(_dt);
+ const PylithScalar dt = _dt;
const int spaceDim = _quadrature->spaceDim();
+ const int indexN = spaceDim - 1;
// Allocate arrays for vertex values
scalar_array tractionTpdtVertex(spaceDim);
+ scalar_array dTractionTpdtVertex(spaceDim);
+ scalar_array dDispRelVertex(spaceDim);
// Get sections
scalar_array slipVertex(spaceDim);
@@ -491,9 +545,9 @@
scalar_array dDispTIncrVertexN(spaceDim);
scalar_array dDispTIncrVertexP(spaceDim);
- const ALE::Obj<RealSection>& dispTIncrSection =
+ const ALE::Obj<RealSection>& dispIncrSection =
fields->get("dispIncr(t->t+dt)").section();
- assert(!dispTIncrSection.isNull());
+ assert(!dispIncrSection.isNull());
scalar_array dLagrangeTpdtVertex(spaceDim);
scalar_array dLagrangeTpdtVertexGlobal(spaceDim);
@@ -524,24 +578,36 @@
#if 0 // DEBUGGING
dispRelSection->view("BEFORE RELATIVE DISPLACEMENT");
- dispTIncrSection->view("BEFORE DISP INCR (t->t+dt)");
+ dispIncrSection->view("BEFORE DISP INCR (t->t+dt)");
#endif
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 relative displacement
- assert(spaceDim == dispRelSection->getFiberDimension(v_fault));
- const PylithScalar* dispRelVertex = dispRelSection->restrictPoint(v_fault);
- assert(dispRelVertex);
+ // Get displacement values
+ assert(spaceDim == dispTSection->getFiberDimension(v_negative));
+ const PylithScalar* dispTVertexN = dispTSection->restrictPoint(v_negative);
+ assert(dispTVertexN);
- // Get relative velocity
- assert(spaceDim == velRelSection->getFiberDimension(v_fault));
- const PylithScalar* velRelVertex = velRelSection->restrictPoint(v_fault);
- assert(velRelVertex);
+ assert(spaceDim == dispTSection->getFiberDimension(v_positive));
+ const PylithScalar* dispTVertexP = dispTSection->restrictPoint(v_positive);
+ assert(dispTVertexP);
+ // Get displacement increment values.
+ assert(spaceDim == dispIncrSection->getFiberDimension(v_negative));
+ const PylithScalar* dispIncrVertexN =
+ dispIncrSection->restrictPoint(v_negative);
+ assert(dispIncrVertexN);
+
+ assert(spaceDim == dispIncrSection->getFiberDimension(v_positive));
+ const PylithScalar* dispIncrVertexP =
+ dispIncrSection->restrictPoint(v_positive);
+ assert(dispIncrVertexP);
+
// Get orientation
assert(spaceDim*spaceDim == orientationSection->getFiberDimension(v_fault));
const PylithScalar* orientationVertex =
@@ -552,11 +618,15 @@
const PylithScalar* lagrangeTVertex = dispTSection->restrictPoint(v_lagrange);
assert(lagrangeTVertex);
- assert(spaceDim == dispTIncrSection->getFiberDimension(v_lagrange));
+ assert(spaceDim == dispIncrSection->getFiberDimension(v_lagrange));
const PylithScalar* lagrangeTIncrVertex =
- dispTIncrSection->restrictPoint(v_lagrange);
+ dispIncrSection->restrictPoint(v_lagrange);
assert(lagrangeTIncrVertex);
+ // 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.
slipVertex = 0.0;
@@ -565,14 +635,60 @@
for (int iDim=0; iDim < spaceDim; ++iDim) {
for (int jDim=0; jDim < spaceDim; ++jDim) {
slipVertex[iDim] += orientationVertex[iDim*spaceDim+jDim] *
- dispRelVertex[jDim];
+ (dispTVertexP[jDim] + dispIncrVertexP[jDim]
+ - dispTVertexN[jDim] - dispIncrVertexN[jDim]);
slipRateVertex[iDim] += orientationVertex[iDim*spaceDim+jDim] *
- velRelVertex[jDim];
+ (dispIncrVertexP[jDim] - dispIncrVertexN[jDim]) / dt;
tractionTpdtVertex[iDim] += orientationVertex[iDim*spaceDim+jDim] *
(lagrangeTVertex[jDim] + lagrangeTIncrVertex[jDim]);
} // for
+ if (fabs(slipRateVertex[iDim]) < _zeroTolerance) {
+ slipRateVertex[iDim] = 0.0;
+ } // if
} // for
+ if (fabs(slipVertex[indexN]) < _zeroTolerance) {
+ slipVertex[indexN] = 0.0;
+ } // if
+ PylithScalar dSlipVertexNormal = 0.0;
+ PylithScalar dTractionTpdtVertexNormal = 0.0;
+ if (slipVertex[indexN]*tractionTpdtVertex[indexN] < 0.0) {
+#if 0 // DEBUGGING
+ std::cout << "STEP 1 CORRECTING NONPHYSICAL SLIP/TRACTIONS"
+ << ", v_fault: " << v_fault
+ << ", slipNormal: " << slipVertex[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(slipVertex[indexN]) > fabs(tractionTpdtVertex[indexN])) {
+ // slip is bigger, so force normal traction back to zero
+ dTractionTpdtVertexNormal = -tractionTpdtVertex[indexN];
+ tractionTpdtVertex[indexN] = 0.0;
+ } else {
+ // traction is bigger, so force slip back to zero
+ dSlipVertexNormal = -slipVertex[indexN];
+ slipVertex[indexN] = 0.0;
+ } // if/else
+ } // if
+ if (slipVertex[indexN] < 0.0) {
+#if 0 // DEBUGGING
+ std::cout << "STEP 1 CORRECTING INTERPENETRATION"
+ << ", v_fault: " << v_fault
+ << ", slipNormal: " << slipVertex[indexN]
+ << ", tractionNormal: " << tractionTpdtVertex[indexN]
+ << std::endl;
+#endif
+ dSlipVertexNormal = -slipVertex[indexN];
+ slipVertex[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);
@@ -586,13 +702,18 @@
tractionTpdtVertex,
iterating);
- // Rotate traction back to global coordinate system.
+ // Rotate increment in traction back to global coordinate system.
dLagrangeTpdtVertexGlobal = 0.0;
for (int iDim=0; iDim < spaceDim; ++iDim) {
for (int jDim=0; jDim < spaceDim; ++jDim) {
dLagrangeTpdtVertexGlobal[iDim] +=
orientationVertex[jDim*spaceDim+iDim] * dLagrangeTpdtVertex[jDim];
} // for
+
+ // Add in potential contribution from adjusting Lagrange
+ // multiplier for fault normal DOF of trial solution in Step 1.
+ dLagrangeTpdtVertexGlobal[iDim] +=
+ orientationVertex[indexN*spaceDim+iDim] * dTractionTpdtVertexNormal;
} // for
#if 0 // debugging
@@ -620,11 +741,39 @@
std::cout << std::endl;
#endif
+ // Set change in Lagrange multiplier
assert(dLagrangeTpdtVertexGlobal.size() ==
dLagrangeTpdtSection->getFiberDimension(v_fault));
dLagrangeTpdtSection->updatePoint(v_fault, &dLagrangeTpdtVertexGlobal[0]);
+
+ // Update displacement in trial solution (if necessary) so that it
+ // conforms to physical constraints.
+ if (0.0 != dSlipVertexNormal) {
+ // Compute relative displacement from slip.
+ dDispRelVertex = 0.0;
+ for (int iDim=0; iDim < spaceDim; ++iDim) {
+ dDispRelVertex[iDim] +=
+ orientationVertex[indexN*spaceDim+iDim] * dSlipVertexNormal;
+
+ dDispTIncrVertexN[iDim] = -0.5*dDispRelVertex[iDim];
+ dDispTIncrVertexP[iDim] = +0.5*dDispRelVertex[iDim];
+ } // for
+
+ // Update displacement field
+ assert(dDispTIncrVertexN.size() ==
+ dispIncrSection->getFiberDimension(v_negative));
+ dispIncrSection->updateAddPoint(v_negative, &dDispTIncrVertexN[0]);
+
+ assert(dDispTIncrVertexP.size() ==
+ dispIncrSection->getFiberDimension(v_positive));
+ dispIncrSection->updateAddPoint(v_positive, &dDispTIncrVertexP[0]);
+ } // if
+
} // for
+ // Step 3: Calculate change in displacement field corresponding to
+ // change in Lagrange multipliers imposed by friction criterion.
+
// Solve sensitivity problem for negative side of the fault.
bool negativeSide = true;
_sensitivityUpdateJacobian(negativeSide, jacobian, *fields);
@@ -639,13 +788,14 @@
_sensitivitySolve();
_sensitivityUpdateSoln(negativeSide);
+ // 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).
scalar_array dSlipVertex(spaceDim);
- scalar_array dDispRelVertex(spaceDim);
scalar_array dispRelVertex(spaceDim);
- // Update slip field based on solution of sensitivity problem and
- // increment in Lagrange multipliers.
const ALE::Obj<RealSection>& sensDispRelSection =
_fields->get("sensitivity relative disp").section();
for (int iVertex=0; iVertex < numVertices; ++iVertex) {
@@ -654,6 +804,19 @@
const int v_negative = _cohesiveVertices[iVertex].negative;
const int v_positive = _cohesiveVertices[iVertex].positive;
+ // Get change in Lagrange multiplier computed from friction criterion.
+ dLagrangeTpdtSection->restrictPoint(v_fault, &dLagrangeTpdtVertex[0],
+ dLagrangeTpdtVertex.size());
+
+ // Solution only changes if Lagrange multiplier changed, so skip
+ // updates if change in Lagrange multiplier is zero.
+ PylithScalar dLagrangeMag = 0.0;
+ for (int iDim=0; iDim < spaceDim; ++iDim)
+ dLagrangeMag += dLagrangeTpdtVertex[iDim]*dLagrangeTpdtVertex[iDim];
+ if (0.0 == dLagrangeMag) {
+ continue; // No change, so continue
+ } // if
+
// Get change in relative displacement from sensitivity solve.
assert(spaceDim == sensDispRelSection->getFiberDimension(v_fault));
const PylithScalar* sensDispRelVertex =
@@ -670,7 +833,7 @@
orientationSection->restrictPoint(v_fault);
assert(orientationVertex);
- // Get displacements from disp(t) and dispIncr(t).
+ // Get displacement.
assert(spaceDim == dispTSection->getFiberDimension(v_negative));
const PylithScalar* dispTVertexN = dispTSection->restrictPoint(v_negative);
assert(dispTVertexN);
@@ -679,94 +842,122 @@
const PylithScalar* dispTVertexP = dispTSection->restrictPoint(v_positive);
assert(dispTVertexP);
- assert(spaceDim == dispTIncrSection->getFiberDimension(v_negative));
- const PylithScalar* dispTIncrVertexN =
- dispTIncrSection->restrictPoint(v_negative);
- assert(dispTIncrVertexN);
+ // Get displacement increment (trial solution).
+ assert(spaceDim == dispIncrSection->getFiberDimension(v_negative));
+ const PylithScalar* dispIncrVertexN =
+ dispIncrSection->restrictPoint(v_negative);
+ assert(dispIncrVertexN);
- assert(spaceDim == dispTIncrSection->getFiberDimension(v_positive));
- const PylithScalar* dispTIncrVertexP =
- dispTIncrSection->restrictPoint(v_positive);
- assert(dispTIncrVertexP);
+ assert(spaceDim == dispIncrSection->getFiberDimension(v_positive));
+ const PylithScalar* dispIncrVertexP =
+ dispIncrSection->restrictPoint(v_positive);
+ assert(dispIncrVertexP);
// Get Lagrange multiplier at time t
assert(spaceDim == dispTSection->getFiberDimension(v_lagrange));
const PylithScalar* lagrangeTVertex = dispTSection->restrictPoint(v_lagrange);
assert(lagrangeTVertex);
- // Get Lagrange multiplier increment at time t
- assert(spaceDim == dispTIncrSection->getFiberDimension(v_lagrange));
+ // Get Lagrange multiplier increment (trial solution)
+ assert(spaceDim == dispIncrSection->getFiberDimension(v_lagrange));
const PylithScalar* lagrangeTIncrVertex =
- dispTIncrSection->restrictPoint(v_lagrange);
+ dispIncrSection->restrictPoint(v_lagrange);
assert(lagrangeTIncrVertex);
- // Get change in Lagrange multiplier.
- dLagrangeTpdtSection->restrictPoint(v_fault, &dLagrangeTpdtVertex[0],
- dLagrangeTpdtVertex.size());
+ // 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).
- // Only update slip, disp, etc if Lagrange multiplier is
- // changing. If Lagrange multiplier is the same, there is no slip
- // so disp, etc should be unchanged.
- PylithScalar dLagrangeMag = 0.0;
- for (int iDim=0; iDim < spaceDim; ++iDim)
- dLagrangeMag += dLagrangeTpdtVertex[iDim]*dLagrangeTpdtVertex[iDim];
- if (0.0 == dLagrangeMag) {
- continue; // No change, so continue
- } // if
-
- // Compute slip and change in slip in fault coordinates.
+ // Compute slip, change in slip, and tractions in fault coordinates.
dSlipVertex = 0.0;
slipVertex = 0.0;
+ tractionTpdtVertex = 0.0;
+ dTractionTpdtVertex = 0.0;
for (int iDim=0; iDim < spaceDim; ++iDim) {
for (int jDim=0; jDim < spaceDim; ++jDim) {
dSlipVertex[iDim] += orientationVertex[iDim*spaceDim+jDim] *
sensDispRelVertex[jDim];
slipVertex[iDim] += orientationVertex[iDim*spaceDim+jDim] *
(dispTVertexP[jDim] - dispTVertexN[jDim] +
- dispTIncrVertexP[jDim] - dispTIncrVertexN[jDim]);
+ dispIncrVertexP[jDim] - dispIncrVertexN[jDim]);
+ tractionTpdtVertex[iDim] += orientationVertex[iDim*spaceDim+jDim] *
+ (lagrangeTVertex[jDim] + lagrangeTIncrVertex[jDim]);
+ dTractionTpdtVertex[iDim] += orientationVertex[iDim*spaceDim+jDim] *
+ dLagrangeTpdtVertex[jDim];
} // for
} // for
+ if (fabs(slipVertex[indexN]) < _zeroTolerance) {
+ slipVertex[indexN] = 0.0;
+ } // if
+ if (fabs(dSlipVertex[indexN]) < _zeroTolerance) {
+ dSlipVertex[indexN] = 0.0;
+ } // if
- // Compute normal traction in fault coordinates.
- PylithScalar tractionNormal = 0.0;
- const int indexN = spaceDim - 1;
- for (int jDim=0; jDim < spaceDim; ++jDim) {
- tractionNormal += orientationVertex[indexN*spaceDim+jDim] *
- (lagrangeTVertex[jDim] + lagrangeTIncrVertex[jDim] +
- dLagrangeTpdtVertex[jDim]);
- } // for
+ if ((slipVertex[indexN] + dSlipVertex[indexN]) *
+ (tractionTpdtVertex[indexN] + dTractionTpdtVertex[indexN])
+ < 0.0) {
+#if 0 // DEBUGGING
+ std::cout << "STEP 4a CORRECTING NONPHYSICAL SLIP/TRACTIONS"
+ << ", v_fault: " << v_fault
+ << ", slipNormal: " << slipVertex[indexN] + dSlipVertex[indexN]
+ << ", tractionNormal: " << tractionTpdtVertex[indexN] + dTractionTpdtVertex[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(slipVertex[indexN] + dSlipVertex[indexN]) >
+ fabs(tractionTpdtVertex[indexN] + dTractionTpdtVertex[indexN])) {
+ // slip is bigger, so force normal traction back to zero
+ dTractionTpdtVertex[indexN] = -tractionTpdtVertex[indexN];
+ } else {
+ // traction is bigger, so force slip back to zero
+ dSlipVertex[indexN] = -slipVertex[indexN];
+ } // if/else
- // Do not allow fault interpenetration and set fault opening to
- // zero if fault is under compression.
- if (tractionNormal < -_zeroTolerance ||
- slipVertex[indexN] + dSlipVertex[indexN] < 0.0) {
+ } // if
+ // Do not allow fault interpenetration.
+ if (slipVertex[indexN] + dSlipVertex[indexN] < 0.0) {
+#if 0 // DEBUGGING
+ std::cout << "STEP 4a CORRECTING INTERPENETATION"
+ << ", v_fault: " << v_fault
+ << ", slipNormal: " << slipVertex[indexN] + dSlipVertex[indexN]
+ << std::endl;
+#endif
dSlipVertex[indexN] = -slipVertex[indexN];
} // if
- // Compute current estimate of slip.
- for (int iDim=0; iDim < spaceDim; ++iDim) {
- slipVertex[iDim] = slipVertex[iDim] + dSlipVertex[iDim];
- } // for
+ // Update current estimate of slip from t to t+dt.
+ slipVertex += dSlipVertex;
- // Update relative displacement from slip.
+ // Compute relative displacement from slip.
dispRelVertex = 0.0;
dDispRelVertex = 0.0;
+ dLagrangeTpdtVertex = 0.0;
for (int iDim=0; iDim < spaceDim; ++iDim) {
for (int jDim=0; jDim < spaceDim; ++jDim) {
dispRelVertex[iDim] += orientationVertex[jDim*spaceDim+iDim] *
slipVertex[jDim];
dDispRelVertex[iDim] += orientationVertex[jDim*spaceDim+iDim] *
dSlipVertex[jDim];
+ dLagrangeTpdtVertex[iDim] += orientationVertex[jDim*spaceDim+iDim] *
+ dTractionTpdtVertex[jDim];
} // for
dDispTIncrVertexN[iDim] = -0.5*dDispRelVertex[iDim];
dDispTIncrVertexP[iDim] = +0.5*dDispRelVertex[iDim];
-
} // for
#if 0 // debugging
- std::cout << "dLagrangeTpdtVertex: ";
+ std::cout << "v_fault: " << v_fault;
+ std::cout << ", tractionTpdtVertex: ";
for (int iDim=0; iDim < spaceDim; ++iDim)
+ std::cout << " " << tractionTpdtVertex[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 << ", slipVertex: ";
for (int iDim=0; iDim < spaceDim; ++iDim)
@@ -784,27 +975,27 @@
assert(dispRelVertex.size() ==
dispRelSection->getFiberDimension(v_fault));
dispRelSection->updatePoint(v_fault, &dispRelVertex[0]);
-
+
// Update Lagrange multiplier increment.
assert(dLagrangeTpdtVertex.size() ==
- dispTIncrSection->getFiberDimension(v_lagrange));
- dispTIncrSection->updateAddPoint(v_lagrange, &dLagrangeTpdtVertex[0]);
+ dispIncrSection->getFiberDimension(v_lagrange));
+ dispIncrSection->updateAddPoint(v_lagrange, &dLagrangeTpdtVertex[0]);
// Update displacement field
assert(dDispTIncrVertexN.size() ==
- dispTIncrSection->getFiberDimension(v_negative));
- dispTIncrSection->updateAddPoint(v_negative, &dDispTIncrVertexN[0]);
+ dispIncrSection->getFiberDimension(v_negative));
+ dispIncrSection->updateAddPoint(v_negative, &dDispTIncrVertexN[0]);
assert(dDispTIncrVertexP.size() ==
- dispTIncrSection->getFiberDimension(v_positive));
- dispTIncrSection->updateAddPoint(v_positive, &dDispTIncrVertexP[0]);
+ dispIncrSection->getFiberDimension(v_positive));
+ dispIncrSection->updateAddPoint(v_positive, &dDispTIncrVertexP[0]);
} // for
#if 0 // DEBUGGING
//dLagrangeTpdtSection->view("AFTER dLagrange");
dispRelSection->view("AFTER RELATIVE DISPLACEMENT");
- dispTIncrSection->view("AFTER DISP INCR (t->t+dt)");
+ dispIncrSection->view("AFTER DISP INCR (t->t+dt)");
#endif
} // constrainSolnSpace
@@ -883,16 +1074,16 @@
const ALE::Obj<RealSection>& dispTSection = fields->get("disp(t)").section();
assert(!dispTSection.isNull());
- scalar_array dispTIncrVertexN(spaceDim);
- scalar_array dispTIncrVertexP(spaceDim);
+ scalar_array dispIncrVertexN(spaceDim);
+ scalar_array dispIncrVertexP(spaceDim);
scalar_array lagrangeTIncrVertex(spaceDim);
- const ALE::Obj<RealSection>& dispTIncrSection =
+ const ALE::Obj<RealSection>& dispIncrSection =
fields->get("dispIncr(t->t+dt)").section();
- assert(!dispTIncrSection.isNull());
+ assert(!dispIncrSection.isNull());
- const ALE::Obj<RealSection>& dispTIncrAdjSection = fields->get(
+ const ALE::Obj<RealSection>& dispIncrAdjSection = fields->get(
"dispIncr adjust").section();
- assert(!dispTIncrAdjSection.isNull());
+ assert(!dispIncrAdjSection.isNull());
const ALE::Obj<RealSection>& jacobianSection = jacobian.section();
assert(!jacobianSection.isNull());
@@ -970,11 +1161,11 @@
assert(lagrangeTVertex);
// Get dispIncr(t) at cohesive cell's vertices.
- dispTIncrSection->restrictPoint(v_negative, &dispTIncrVertexN[0],
- dispTIncrVertexN.size());
- dispTIncrSection->restrictPoint(v_positive, &dispTIncrVertexP[0],
- dispTIncrVertexP.size());
- dispTIncrSection->restrictPoint(v_lagrange, &lagrangeTIncrVertex[0],
+ dispIncrSection->restrictPoint(v_negative, &dispIncrVertexN[0],
+ dispIncrVertexN.size());
+ dispIncrSection->restrictPoint(v_positive, &dispIncrVertexP[0],
+ dispIncrVertexP.size());
+ dispIncrSection->restrictPoint(v_lagrange, &lagrangeTIncrVertex[0],
lagrangeTIncrVertex.size());
// Get relative displacement at fault vertex.
@@ -1008,14 +1199,14 @@
assert(S > 0.0);
lagrangeTIncrVertex[iDim] = 1.0/S *
(-residualVertexL[iDim] +
- areaVertex * (dispTIncrVertexP[iDim] - dispTIncrVertexN[iDim]));
+ areaVertex * (dispIncrVertexP[iDim] - dispIncrVertexN[iDim]));
assert(jacobianVertexN[iDim] > 0.0);
- dispTIncrVertexN[iDim] =
+ dispIncrVertexN[iDim] =
+areaVertex / jacobianVertexN[iDim]*lagrangeTIncrVertex[iDim];
assert(jacobianVertexP[iDim] > 0.0);
- dispTIncrVertexP[iDim] =
+ dispIncrVertexP[iDim] =
-areaVertex / jacobianVertexP[iDim]*lagrangeTIncrVertex[iDim];
} // for
@@ -1059,12 +1250,12 @@
} // for
#if 0 // debugging
- std::cout << "dispTIncrP: ";
+ std::cout << "dispIncrP: ";
for (int iDim=0; iDim < spaceDim; ++iDim)
- std::cout << " " << dispTIncrVertexP[iDim];
- std::cout << ", dispTIncrN: ";
+ std::cout << " " << dispIncrVertexP[iDim];
+ std::cout << ", dispIncrN: ";
for (int iDim=0; iDim < spaceDim; ++iDim)
- std::cout << " " << dispTIncrVertexN[iDim];
+ std::cout << " " << dispIncrVertexN[iDim];
std::cout << ", slipVertex: ";
for (int iDim=0; iDim < spaceDim; ++iDim)
std::cout << " " << slipVertex[iDim];
@@ -1097,9 +1288,9 @@
assert(jacobianVertexP[iDim] > 0.0);
assert(jacobianVertexN[iDim] > 0.0);
- dispTIncrVertexN[iDim] +=
+ dispIncrVertexN[iDim] +=
areaVertex * dLagrangeTpdtVertexGlobal[iDim] / jacobianVertexN[iDim];
- dispTIncrVertexP[iDim] -=
+ dispIncrVertexP[iDim] -=
areaVertex * dLagrangeTpdtVertexGlobal[iDim] / jacobianVertexP[iDim];
// Set increment in relative displacement.
@@ -1120,13 +1311,13 @@
if (globalOrder->isLocal(v_lagrange)) {
// Adjust displacements to account for Lagrange multiplier values
// (assumed to be zero in preliminary solve).
- assert(dispTIncrVertexN.size() ==
- dispTIncrAdjSection->getFiberDimension(v_negative));
- dispTIncrAdjSection->updateAddPoint(v_negative, &dispTIncrVertexN[0]);
+ assert(dispIncrVertexN.size() ==
+ dispIncrAdjSection->getFiberDimension(v_negative));
+ dispIncrAdjSection->updateAddPoint(v_negative, &dispIncrVertexN[0]);
- assert(dispTIncrVertexP.size() ==
- dispTIncrAdjSection->getFiberDimension(v_positive));
- dispTIncrAdjSection->updateAddPoint(v_positive, &dispTIncrVertexP[0]);
+ assert(dispIncrVertexP.size() ==
+ dispIncrAdjSection->getFiberDimension(v_positive));
+ dispIncrAdjSection->updateAddPoint(v_positive, &dispIncrVertexP[0]);
} // if
// The Lagrange multiplier and relative displacement are NOT
@@ -1135,8 +1326,8 @@
// Set Lagrange multiplier value. Value from preliminary solve is
// bogus due to artificial diagonal entry in Jacobian of 1.0.
assert(lagrangeTIncrVertex.size() ==
- dispTIncrSection->getFiberDimension(v_lagrange));
- dispTIncrSection->updatePoint(v_lagrange, &lagrangeTIncrVertex[0]);
+ dispIncrSection->getFiberDimension(v_lagrange));
+ dispIncrSection->updatePoint(v_lagrange, &lagrangeTIncrVertex[0]);
// Update the relative displacement estimate based on adjustment
// to the Lagrange multiplier values.
@@ -1152,14 +1343,13 @@
9 + // updates
spaceDim*9));
-
#if !defined(DETAILED_EVENT_LOGGING)
_logger->eventEnd(computeEvent);
#endif
#if 0 // DEBUGGING
//dLagrangeTpdtSection->view("AFTER dLagrange");
- //dispTIncrSection->view("AFTER DISP INCR (t->t+dt)");
+ //dispIncrSection->view("AFTER DISP INCR (t->t+dt)");
dispRelSection->view("AFTER RELATIVE DISPLACEMENT");
//velRelSection->view("AFTER RELATIVE VELOCITY");
#endif
@@ -1494,16 +1684,30 @@
} // _calcTractions
// ----------------------------------------------------------------------
-// Update slip rate associated with Lagrange vertex k corresponding
-// to diffential velocity between conventional vertices i and j.
+// Update relative displacement and velocity (slip and slip rate)
+// associated with Lagrange vertex k corresponding to diffential
+// velocity between conventional vertices i and j.
void
-pylith::faults::FaultCohesiveDyn::_updateVelRel(const topology::SolutionFields& fields)
-{ // _updateVelRel
+pylith::faults::FaultCohesiveDyn::_updateRelMotion(const topology::SolutionFields& fields)
+{ // _updateRelMotion
assert(0 != _fields);
const int spaceDim = _quadrature->spaceDim();
// Get section information
+ 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());
+
+ scalar_array dispRelVertex(spaceDim);
+ const ALE::Obj<RealSection>& dispRelSection =
+ _fields->get("relative disp").section();
+ assert(!dispRelSection.isNull());
+
const ALE::Obj<RealSection>& velocitySection =
fields.get("velocity(t)").section();
assert(!velocitySection.isNull());
@@ -1511,22 +1715,54 @@
scalar_array velRelVertex(spaceDim);
const ALE::Obj<RealSection>& velRelSection =
_fields->get("relative velocity").section();
+ assert(!velRelSection.isNull());
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 values
+ // Get displacement values
+ assert(spaceDim == dispTSection->getFiberDimension(v_negative));
+ const PylithScalar* dispTVertexN = dispTSection->restrictPoint(v_negative);
+ assert(dispTVertexN);
+
+ assert(spaceDim == dispTSection->getFiberDimension(v_positive));
+ const PylithScalar* dispTVertexP = dispTSection->restrictPoint(v_positive);
+ assert(dispTVertexP);
+
+ assert(spaceDim == dispIncrSection->getFiberDimension(v_negative));
+ const PylithScalar* dispIncrVertexN =
+ dispIncrSection->restrictPoint(v_negative);
+ assert(dispIncrVertexN);
+
+ assert(spaceDim == dispIncrSection->getFiberDimension(v_positive));
+ const PylithScalar* dispIncrVertexP =
+ dispIncrSection->restrictPoint(v_positive);
+ assert(dispIncrVertexP);
+
+ // Compute relative displacememt
+ for (int iDim=0; iDim < spaceDim; ++iDim) {
+ const PylithScalar value =
+ dispTVertexP[iDim] + dispIncrVertexP[iDim]
+ - dispTVertexN[iDim] - dispIncrVertexN[iDim];
+ dispRelVertex[iDim] = fabs(value) > _zeroTolerance ? value : 0.0;
+ } // for
+
+ // Update relative displacement field.
+ assert(dispRelVertex.size() ==
+ dispRelSection->getFiberDimension(v_fault));
+ dispRelSection->updatePoint(v_fault, &dispRelVertex[0]);
+
+ // Get velocity values
+ assert(spaceDim == velocitySection->getFiberDimension(v_negative));
const PylithScalar* velocityVertexN = velocitySection->restrictPoint(v_negative);
- assert(0 != velocityVertexN);
- assert(spaceDim == velocitySection->getFiberDimension(v_negative));
+ assert(velocityVertexN);
+ assert(spaceDim == velocitySection->getFiberDimension(v_positive));
const PylithScalar* velocityVertexP = velocitySection->restrictPoint(v_positive);
- assert(0 != velocityVertexP);
- assert(spaceDim == velocitySection->getFiberDimension(v_positive));
+ assert(velocityVertexP);
// Compute relative velocity
for (int iDim=0; iDim < spaceDim; ++iDim) {
@@ -1534,14 +1770,14 @@
velRelVertex[iDim] = fabs(value) > _zeroTolerance ? value : 0.0;
} // for
- // Update slip rate field.
+ // Update relative velocity field.
assert(velRelVertex.size() ==
velRelSection->getFiberDimension(v_fault));
velRelSection->updatePoint(v_fault, &velRelVertex[0]);
} // for
PetscLogFlops(numVertices*spaceDim*spaceDim*4);
-} // _updateVelRel
+} // _updateRelMotion
// ----------------------------------------------------------------------
// Setup sensitivity problem to compute change in slip given change in Lagrange multipliers.
@@ -1611,7 +1847,7 @@
int maxIters = 0;
err = KSPGetTolerances(_ksp, &rtol, &atol, &dtol, &maxIters);
CHECK_PETSC_ERROR(err);
- rtol = 1.0e-2*_zeroTolerance;
+ rtol = 1.0e-3*_zeroTolerance;
atol = 1.0e-5*_zeroTolerance;
err = KSPSetTolerances(_ksp, rtol, atol, dtol, maxIters);
CHECK_PETSC_ERROR(err);
@@ -1932,16 +2168,16 @@
{ // _constrainSolnSpace1D
assert(0 != dLagrangeTpdt);
- if (tractionTpdt[0] < 0) {
- // if compression, then no changes to solution
- } else {
- // if tension, then traction is zero.
-
- const PylithScalar dlp = -tractionTpdt[0];
- (*dLagrangeTpdt)[0] = dlp;
- } // else
-
- PetscLogFlops(2);
+ if (fabs(slip[0]) < _zeroTolerance) {
+ // if compression, then no changes to solution
+ } else {
+ // if tension, then traction is zero.
+
+ const PylithScalar dlp = -tractionTpdt[0];
+ (*dLagrangeTpdt)[0] = dlp;
+ } // else
+
+ PetscLogFlops(2);
} // _constrainSolnSpace1D
// ----------------------------------------------------------------------
@@ -1961,7 +2197,7 @@
const PylithScalar tractionNormal = tractionTpdt[1];
const PylithScalar tractionShearMag = fabs(tractionTpdt[0]);
- if (tractionNormal < 0 && 0.0 == slip[1]) {
+ if (fabs(slip[1]) < _zeroTolerance && tractionNormal < -_zeroTolerance) {
// if in compression and no opening
const PylithScalar frictionStress = _friction->calcFriction(slipMag, slipRateMag,
tractionNormal);
@@ -2015,7 +2251,7 @@
sqrt(tractionTpdt[0] * tractionTpdt[0] +
tractionTpdt[1] * tractionTpdt[1]);
- if (tractionNormal < 0.0 && 0.0 == slip[2]) {
+ if (fabs(slip[2]) < _zeroTolerance && tractionNormal < -_zeroTolerance) {
// if in compression and no opening
const PylithScalar frictionStress =
_friction->calcFriction(slipShearMag, slipRateMag, tractionNormal);
Modified: short/3D/PyLith/trunk/libsrc/pylith/faults/FaultCohesiveDyn.hh
===================================================================
--- short/3D/PyLith/trunk/libsrc/pylith/faults/FaultCohesiveDyn.hh 2011-11-06 22:16:51 UTC (rev 19153)
+++ short/3D/PyLith/trunk/libsrc/pylith/faults/FaultCohesiveDyn.hh 2011-11-07 18:43:31 UTC (rev 19154)
@@ -152,13 +152,13 @@
void _calcTractions(topology::Field<topology::SubMesh>* tractions,
const topology::Field<topology::Mesh>& solution);
- /** Update relative velocity associated with Lagrange vertex k
- * corresponding to diffential velocity between conventional
- * vertices i and j.
+ /** Update relative displacement and velocity associated with
+ * Lagrange vertex k corresponding to diffential velocity between
+ * conventional vertices i and j.
*
* @param fields Solution fields.
*/
- void _updateVelRel(const topology::SolutionFields& fields);
+ void _updateRelMotion(const topology::SolutionFields& fields);
/** Setup sensitivity problem to compute change in slip given change
* in Lagrange multipliers.
Modified: short/3D/PyLith/trunk/libsrc/pylith/problems/Formulation.cc
===================================================================
--- short/3D/PyLith/trunk/libsrc/pylith/problems/Formulation.cc 2011-11-06 22:16:51 UTC (rev 19153)
+++ short/3D/PyLith/trunk/libsrc/pylith/problems/Formulation.cc 2011-11-07 18:43:31 UTC (rev 19154)
@@ -425,5 +425,52 @@
solution += adjust;
} // adjustSolnLumped
+#include "pylith/meshio/DataWriterHDF5.hh"
+// ----------------------------------------------------------------------
+void
+pylith::problems::Formulation::printState(PetscVec* solutionVec,
+ PetscVec* residualVec,
+ PetscVec* solution0Vec,
+ PetscVec* searchDirVec)
+{ // printState
+ assert(solutionVec);
+ assert(residualVec);
+ assert(searchDirVec);
+ meshio::DataWriterHDF5<topology::Mesh,topology::Field<topology::Mesh> > writer;
+
+ const topology::Mesh& mesh = _fields->mesh();
+
+ writer.filename("state.h5");
+ const int numTimeSteps = 1;
+ writer.open(mesh, numTimeSteps);
+
+
+ topology::Field<topology::Mesh>& solution = _fields->solution();
+ solution.scatterVectorToSection(*solutionVec);
+ writer.writeVertexField(0.0, solution, mesh);
+ solution.view("DIVERGED_SOLUTION");
+ const char* label = solution.label();
+
+ solution.label("solution_0");
+ solution.scatterVectorToSection(*solution0Vec);
+ writer.writeVertexField(0.0, solution, mesh);
+ solution.view("DIVERGED_SOLUTION0");
+ solution.label(label);
+
+ topology::Field<topology::Mesh>& residual = _fields->get("residual");
+ residual.scatterVectorToSection(*residualVec);
+ writer.writeVertexField(0.0, residual, mesh);
+ residual.view("DIVERGED_RESIDUAL");
+
+ residual.label("search_dir");
+ residual.scatterVectorToSection(*searchDirVec);
+ writer.writeVertexField(0.0, residual, mesh);
+ residual.view("DIVERGED_SEARCHDIR");
+
+ writer.close();
+} // printState
+
+
+
// End of file
Modified: short/3D/PyLith/trunk/libsrc/pylith/problems/Formulation.hh
===================================================================
--- short/3D/PyLith/trunk/libsrc/pylith/problems/Formulation.hh 2011-11-06 22:16:51 UTC (rev 19153)
+++ short/3D/PyLith/trunk/libsrc/pylith/problems/Formulation.hh 2011-11-07 18:43:31 UTC (rev 19154)
@@ -186,6 +186,12 @@
virtual
void calcRateFields(void) = 0;
+ /// Write state of system
+ void printState(PetscVec* solutionVec,
+ PetscVec* residualVec,
+ PetscVec* solution0Vec,
+ PetscVec* searchDirVec);
+
// PROTECTED MEMBERS ////////////////////////////////////////////////////
protected :
Modified: short/3D/PyLith/trunk/libsrc/pylith/problems/SolverNonlinear.cc
===================================================================
--- short/3D/PyLith/trunk/libsrc/pylith/problems/SolverNonlinear.cc 2011-11-06 22:16:51 UTC (rev 19153)
+++ short/3D/PyLith/trunk/libsrc/pylith/problems/SolverNonlinear.cc 2011-11-07 18:43:31 UTC (rev 19154)
@@ -362,6 +362,10 @@
ierr = PetscViewerASCIISubtractTab(snes->ls_monitor,((PetscObject)snes)->tablevel);CHKERRQ(ierr);
}
*flag = PETSC_FALSE; // DIVERGED_LINE_SEARCH
+ assert(lsctx);
+ Formulation* formulation = (Formulation*) lsctx;
+ assert(formulation);
+ formulation->printState(&w, &g, &x, &y);
break;
}
t1 = .5*((*gnorm)*(*gnorm) - fnorm*fnorm) - lambda*initslope;
@@ -380,11 +384,11 @@
// range. Necessary due to underflow and overflow of a, b, c,
// and d.
#if 0
- lambdatemp = (-b + sqrt(d))/(3.0*a);
+ lambdatemp = (-b + PetscSqrtScalar(d))/(3.0*a);
#else
- if ((-b + sqrt(d) > 0.0 && a > 0.0) ||
- (-b + sqrt(d) < 0.0 && a < 0.0)) {
- lambdatemp = (-b + sqrt(d))/(3.0*a);
+ if ((-b + PetscSqrtScalar(d) > 0.0 && a > 0.0) ||
+ (-b + PetscSqrtScalar(d) < 0.0 && a < 0.0)) {
+ lambdatemp = (-b + PetscSqrtScalar(d))/(3.0*a);
} else {
lambdatemp = 0.05*lambda;
} // else
@@ -454,9 +458,9 @@
// ======================================================================
// Code to constrain solution space.
- assert(0 != lsctx);
+ assert(lsctx);
Formulation* formulation = (Formulation*) lsctx;
- assert(0 != formulation);
+ assert(formulation);
formulation->constrainSolnSpace(&w);
ierr = SNESComputeFunction(snes,w,g);CHKERRQ(ierr);
Modified: short/3D/PyLith/trunk/tests/2d/frictionslide/plot_friction.py
===================================================================
--- short/3D/PyLith/trunk/tests/2d/frictionslide/plot_friction.py 2011-11-06 22:16:51 UTC (rev 19153)
+++ short/3D/PyLith/trunk/tests/2d/frictionslide/plot_friction.py 2011-11-07 18:43:31 UTC (rev 19154)
@@ -1,6 +1,6 @@
#!/usr/bin/env python
-sim = "ratestate_stable"
+sim = "ratestate_weak"
# ======================================================================
import tables
Modified: short/3D/PyLith/trunk/tests/2d/frictionslide/pylithapp.cfg
===================================================================
--- short/3D/PyLith/trunk/tests/2d/frictionslide/pylithapp.cfg 2011-11-06 22:16:51 UTC (rev 19153)
+++ short/3D/PyLith/trunk/tests/2d/frictionslide/pylithapp.cfg 2011-11-07 18:43:31 UTC (rev 19154)
@@ -78,7 +78,7 @@
[pylithapp.timedependent.interfaces]
fault = pylith.faults.FaultCohesiveDyn
-fault.zero_tolerance = 1.0e-14
+fault.zero_tolerance = 1.0e-12
[pylithapp.timedependent.interfaces.fault]
id = 100
Modified: short/3D/PyLith/trunk/tests/2d/frictionslide/ratestate.cfg
===================================================================
--- short/3D/PyLith/trunk/tests/2d/frictionslide/ratestate.cfg 2011-11-06 22:16:51 UTC (rev 19153)
+++ short/3D/PyLith/trunk/tests/2d/frictionslide/ratestate.cfg 2011-11-07 18:43:31 UTC (rev 19154)
@@ -46,6 +46,7 @@
db_initial_tractions.values = [traction-shear,traction-normal]
db_initial_tractions.data = [0.0*MPa, -10.0*MPa]
+zero_tolerance = 1.0e-14
# ----------------------------------------------------------------------
# PETSc
Copied: short/3D/PyLith/trunk/tests/3d/cyclicfriction/fieldsplit.cfg (from rev 19153, short/3D/PyLith/branches/v1.6-stable/tests/3d/cyclicfriction/fieldsplit.cfg)
===================================================================
--- short/3D/PyLith/trunk/tests/3d/cyclicfriction/fieldsplit.cfg (rev 0)
+++ short/3D/PyLith/trunk/tests/3d/cyclicfriction/fieldsplit.cfg 2011-11-07 18:43:31 UTC (rev 19154)
@@ -0,0 +1,18 @@
+[pylithapp.timedependent.formulation]
+split_fields = True
+matrix_type = aij
+use_custom_constraint_pc = True
+
+[pylithapp.petsc]
+ksp_gmres_restart = 500
+fs_pc_type = fieldsplit
+fs_pc_fieldsplit_real_diagonal =
+fs_pc_fieldsplit_type = multiplicative
+fs_fieldsplit_0_pc_type = ml
+fs_fieldsplit_1_pc_type = ml
+fs_fieldsplit_2_pc_type = ml
+fs_fieldsplit_3_pc_type = ml
+fs_fieldsplit_0_ksp_type = preonly
+fs_fieldsplit_1_ksp_type = preonly
+fs_fieldsplit_2_ksp_type = preonly
+fs_fieldsplit_3_ksp_type = preonly
Modified: short/3D/PyLith/trunk/tests/3d/cyclicfriction/geometry.jou
===================================================================
--- short/3D/PyLith/trunk/tests/3d/cyclicfriction/geometry.jou 2011-11-06 22:16:51 UTC (rev 19153)
+++ short/3D/PyLith/trunk/tests/3d/cyclicfriction/geometry.jou 2011-11-07 18:43:31 UTC (rev 19154)
@@ -25,6 +25,7 @@
# Create interface surfaces
# ----------------------------------------------------------------------
create planar surface with plane xplane offset 0
+rotate surface 7 angle -30 about z include_merged
create sheet extended from surface 7 intersecting volume 1
delete body 2
surface 8 name "fault_surface"
@@ -52,3 +53,4 @@
imprint all with volume all
merge all
+
Modified: short/3D/PyLith/trunk/tests/3d/cyclicfriction/hex8_1000m.exo
===================================================================
(Binary files differ)
Modified: short/3D/PyLith/trunk/tests/3d/cyclicfriction/hex8_1000m.jou
===================================================================
--- short/3D/PyLith/trunk/tests/3d/cyclicfriction/hex8_1000m.jou 2011-11-06 22:16:51 UTC (rev 19153)
+++ short/3D/PyLith/trunk/tests/3d/cyclicfriction/hex8_1000m.jou 2011-11-07 18:43:31 UTC (rev 19154)
@@ -6,6 +6,7 @@
# ----------------------------------------------------------------------
# Set discretization size
# ----------------------------------------------------------------------
+surface 10 17 16 12 scheme pave
volume all size 1000
# ----------------------------------------------------------------------
@@ -101,3 +102,4 @@
# Export exodus file
# ----------------------------------------------------------------------
export mesh "hex8_1000m.exo" dimension 3 overwrite
+
Modified: short/3D/PyLith/trunk/tests/3d/cyclicfriction/initial_tractions.spatialdb
===================================================================
--- short/3D/PyLith/trunk/tests/3d/cyclicfriction/initial_tractions.spatialdb 2011-11-06 22:16:51 UTC (rev 19153)
+++ short/3D/PyLith/trunk/tests/3d/cyclicfriction/initial_tractions.spatialdb 2011-11-07 18:43:31 UTC (rev 19154)
@@ -12,4 +12,4 @@
}
}
0.0 0.0 0.0 0.0 0.0 0.0
-0.0 0.0 -6.0 0.0 0.0 -60.0
+0.0 0.0 -6.0 0.0 0.0 -25.0
Modified: short/3D/PyLith/trunk/tests/3d/cyclicfriction/pylithapp.cfg
===================================================================
--- short/3D/PyLith/trunk/tests/3d/cyclicfriction/pylithapp.cfg 2011-11-06 22:16:51 UTC (rev 19153)
+++ short/3D/PyLith/trunk/tests/3d/cyclicfriction/pylithapp.cfg 2011-11-07 18:43:31 UTC (rev 19154)
@@ -38,6 +38,7 @@
[pylithapp.timedependent.normalizer]
relaxation_time = 1.0*hour
+length_scale = 1.0*km
[pylithapp.timedependent.implicit]
solver = pylith.problems.SolverNonlinear
@@ -113,6 +114,7 @@
[pylithapp.timedependent.interfaces.fault]
label = fault
+zero_tolerance = 1.0e-8
friction = pylith.friction.SlipWeakening
friction.label = Slip weakening
@@ -170,8 +172,8 @@
# Convergence parameters.
ksp_rtol = 1.0e-8
ksp_atol = 1.0e-10
-ksp_max_it = 200
-ksp_gmres_restart = 70
+ksp_max_it = 100
+ksp_gmres_restart = 50
# Linear solver monitoring options.
ksp_monitor = true
@@ -180,13 +182,26 @@
# Nonlinear solver monitoring options.
snes_rtol = 1.0e-8
-snes_atol = 1.0e-8
-snes_max_it = 200
+snes_atol = 1.0e-7
+snes_max_it = 50
snes_monitor = true
snes_ls_monitor = true
#snes_view = true
snes_converged_reason = true
-#snes_monitor_solution_update = true
+snes_monitor_solution_update = true
+snes_monitor_residual = true
#info =
+# Friction sensitivity solve used to compute the increment in slip
+# associated with changes in the Lagrange multiplier imposed by the
+# fault constitutive model.
+friction_pc_type = asm
+friction_sub_pc_factor_shift_type = nonzero
+friction_ksp_max_it = 200
+friction_ksp_gmres_restart = 50
+# Uncomment to view details of friction sensitivity solve.
+#friction_ksp_monitor = true
+#friction_ksp_view = true
+friction_ksp_converged_reason = true
+
log_summary = true
More information about the CIG-COMMITS
mailing list