[cig-commits] r19144 - 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
Wed Nov 2 18:09:17 PDT 2011
Author: brad
Date: 2011-11-02 18:09:17 -0700 (Wed, 02 Nov 2011)
New Revision: 19144
Added:
short/3D/PyLith/branches/v1.6-stable/tests/3d/cyclicfriction/fieldsplit.cfg
Modified:
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/problems/Formulation.cc
short/3D/PyLith/branches/v1.6-stable/libsrc/pylith/problems/Formulation.hh
short/3D/PyLith/branches/v1.6-stable/libsrc/pylith/problems/SolverNonlinear.cc
short/3D/PyLith/branches/v1.6-stable/tests/3d/cyclicfriction/initial_tractions.spatialdb
Log:
Adjusted fault friction implementation so that residual doesn't include fault terms at locations where fault opens. Added output of diagnostic info when line search fails.
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 2011-11-02 22:01:24 UTC (rev 19143)
+++ short/3D/PyLith/branches/v1.6-stable/libsrc/pylith/faults/FaultCohesiveDyn.cc 2011-11-03 01:09:17 UTC (rev 19144)
@@ -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,35 @@
// Get sections associated with cohesive cells
double_array residualVertexN(spaceDim);
double_array residualVertexP(spaceDim);
+ double_array residualVertexL(spaceDim);
const ALE::Obj<RealSection>& residualSection = residual.section();
assert(!residualSection.isNull());
- const ALE::Obj<RealSection>& areaSection = _fields->get("area").section();
- assert(!areaSection.isNull());
+ const ALE::Obj<RealSection>& dispTSection = fields->get("disp(t)").section();
+ assert(!dispTSection.isNull());
- const ALE::Obj<RealSection>& initialTractionsSection =
- _fields->get("initial traction").section();
- assert(!initialTractionsSection.isNull());
+ const ALE::Obj<RealSection>& dispTIncrSection =
+ fields->get("dispIncr(t->t+dt)").section();
+ assert(!dispTIncrSection.isNull());
+ double_array dispTpdtVertexN(spaceDim);
+ double_array dispTpdtVertexP(spaceDim);
+ double_array dispTpdtVertexL(spaceDim);
+
const ALE::Obj<RealSection>& dispRelSection =
_fields->get("relative disp").section();
assert(!dispRelSection.isNull());
+ double_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>& orientationSection =
_fields->get("orientation").section();
assert(!orientationSection.isNull());
@@ -243,21 +258,19 @@
_logger->eventBegin(restrictEvent);
#endif
- // Get initial tractions at fault vertex.
- assert(spaceDim == initialTractionsSection->getFiberDimension(v_fault));
- const double* initialTractionsVertex =
- initialTractionsSection->restrictPoint(v_fault);
- assert(initialTractionsVertex);
-
// Get relative dislplacement at fault vertex.
assert(spaceDim == dispRelSection->getFiberDimension(v_fault));
const double* dispRelVertex = dispRelSection->restrictPoint(v_fault);
assert(dispRelVertex);
- // Get area associated with fault vertex.
- assert(1 == areaSection->getFiberDimension(v_fault));
- assert(areaSection->restrictPoint(v_fault));
- const double areaVertex = *areaSection->restrictPoint(v_fault);
+ // Get initial tractions at fault vertex.
+ if (_dbInitialTract) {
+ initialTractionsSection->restrictPoint(v_fault,
+ &initialTractionsVertex[0],
+ initialTractionsVertex.size());
+ } else {
+ initialTractionsVertex = 0.0;
+ } // if/else
// Get orientation associated with fault vertex.
assert(spaceDim*spaceDim == orientationSection->getFiberDimension(v_fault));
@@ -265,32 +278,74 @@
orientationSection->restrictPoint(v_fault);
assert(orientationVertex);
+ // Get area associated with fault vertex.
+ assert(1 == areaSection->getFiberDimension(v_fault));
+ assert(areaSection->restrictPoint(v_fault));
+ const double areaVertex = *areaSection->restrictPoint(v_fault);
+
+ // Get disp(t) at conventional vertices and Lagrange vertex.
+ 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);
+
+ assert(spaceDim == dispTSection->getFiberDimension(v_lagrange));
+ const double* 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 double* dispTIncrVertexN =
+ dispTIncrSection->restrictPoint(v_negative);
+ assert(dispTIncrVertexN);
+
+ assert(spaceDim == dispTIncrSection->getFiberDimension(v_positive));
+ const double* dispTIncrVertexP =
+ dispTIncrSection->restrictPoint(v_positive);
+ assert(dispTIncrVertexP);
+
+ assert(spaceDim == dispTIncrSection->getFiberDimension(v_lagrange));
+ const double* 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 relative displacement.
double slipNormal = 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] * dispRelVertex[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;
+ for (int iDim=0; iDim < spaceDim; ++iDim) {
+ residualVertexL[iDim] = -areaVertex *
+ (dispTpdtVertexP[iDim] - dispTpdtVertexN[iDim] - dispRelVertex[iDim]);
+ } // for
} // if
+ residualVertexP = -residualVertexN;
#if defined(DETAILED_EVENT_LOGGING)
_logger->eventEnd(computeEvent);
@@ -306,11 +361,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 +386,7 @@
assert(0 != fields);
assert(0 != _fields);
- _updateVelRel(*fields);
+ _updateRelMotion(*fields);
const int spaceDim = _quadrature->spaceDim();
@@ -460,7 +519,7 @@
assert(0 != _fields);
assert(0 != _friction);
- _updateVelRel(*fields);
+ _updateRelMotion(*fields);
_sensitivitySetup(jacobian);
// Update time step in friction (can vary).
@@ -1494,16 +1553,32 @@
} // _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
+#if 0 // TEMPORARY test case
+ 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());
+
+ double_array dispRelVertex(spaceDim);
+ const ALE::Obj<RealSection>& dispRelSection =
+ _fields->get("relative disp").section();
+ assert(!dispRelSection.isNull());
+#endif
+
const ALE::Obj<RealSection>& velocitySection =
fields.get("velocity(t)").section();
assert(!velocitySection.isNull());
@@ -1511,22 +1586,56 @@
double_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
+#if 0 // TEMPORARY test case
+ // 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);
+
+ assert(spaceDim == dispTIncrSection->getFiberDimension(v_negative));
+ const double* dispTIncrVertexN =
+ dispTIncrSection->restrictPoint(v_negative);
+ assert(dispTIncrVertexN);
+
+ assert(spaceDim == dispTIncrSection->getFiberDimension(v_positive));
+ const double* dispTIncrVertexP =
+ dispTIncrSection->restrictPoint(v_positive);
+ assert(dispTIncrVertexP);
+
+ // Compute relative displacememt
+ for (int iDim=0; iDim < spaceDim; ++iDim) {
+ const double value =
+ dispTVertexP[iDim] + dispTIncrVertexP[iDim]
+ - dispTVertexN[iDim] - dispTIncrVertexN[iDim];
+ dispRelVertex[iDim] = value > _zeroTolerance ? value : 0.0;
+ } // for
+
+ // Update relative displacement field.
+ assert(dispRelVertex.size() ==
+ dispRelSection->getFiberDimension(v_fault));
+ dispRelSection->updatePoint(v_fault, &dispRelVertex[0]);
+#endif
+
+ // Get velocity values
+ assert(spaceDim == velocitySection->getFiberDimension(v_negative));
const double* velocityVertexN = velocitySection->restrictPoint(v_negative);
- assert(0 != velocityVertexN);
- assert(spaceDim == velocitySection->getFiberDimension(v_negative));
+ assert(velocityVertexN);
+ assert(spaceDim == velocitySection->getFiberDimension(v_positive));
const double* 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 +1643,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.
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 2011-11-02 22:01:24 UTC (rev 19143)
+++ short/3D/PyLith/branches/v1.6-stable/libsrc/pylith/faults/FaultCohesiveDyn.hh 2011-11-03 01:09:17 UTC (rev 19144)
@@ -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/branches/v1.6-stable/libsrc/pylith/problems/Formulation.cc
===================================================================
--- short/3D/PyLith/branches/v1.6-stable/libsrc/pylith/problems/Formulation.cc 2011-11-02 22:01:24 UTC (rev 19143)
+++ short/3D/PyLith/branches/v1.6-stable/libsrc/pylith/problems/Formulation.cc 2011-11-03 01:09:17 UTC (rev 19144)
@@ -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/branches/v1.6-stable/libsrc/pylith/problems/Formulation.hh
===================================================================
--- short/3D/PyLith/branches/v1.6-stable/libsrc/pylith/problems/Formulation.hh 2011-11-02 22:01:24 UTC (rev 19143)
+++ short/3D/PyLith/branches/v1.6-stable/libsrc/pylith/problems/Formulation.hh 2011-11-03 01:09:17 UTC (rev 19144)
@@ -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/branches/v1.6-stable/libsrc/pylith/problems/SolverNonlinear.cc
===================================================================
--- short/3D/PyLith/branches/v1.6-stable/libsrc/pylith/problems/SolverNonlinear.cc 2011-11-02 22:01:24 UTC (rev 19143)
+++ short/3D/PyLith/branches/v1.6-stable/libsrc/pylith/problems/SolverNonlinear.cc 2011-11-03 01:09:17 UTC (rev 19144)
@@ -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);
Added: short/3D/PyLith/branches/v1.6-stable/tests/3d/cyclicfriction/fieldsplit.cfg
===================================================================
--- short/3D/PyLith/branches/v1.6-stable/tests/3d/cyclicfriction/fieldsplit.cfg (rev 0)
+++ short/3D/PyLith/branches/v1.6-stable/tests/3d/cyclicfriction/fieldsplit.cfg 2011-11-03 01:09:17 UTC (rev 19144)
@@ -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/branches/v1.6-stable/tests/3d/cyclicfriction/initial_tractions.spatialdb
===================================================================
--- short/3D/PyLith/branches/v1.6-stable/tests/3d/cyclicfriction/initial_tractions.spatialdb 2011-11-02 22:01:24 UTC (rev 19143)
+++ short/3D/PyLith/branches/v1.6-stable/tests/3d/cyclicfriction/initial_tractions.spatialdb 2011-11-03 01:09:17 UTC (rev 19144)
@@ -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 -30.0
+0.0 0.0 -6.0 0.0 0.0 -25.0
More information about the CIG-COMMITS
mailing list