[cig-commits] [commit] baagaard/dynrup-new-lagrange: Fixed some bugs in computing residual in new FaultCohesiveDyn implementation. (ba6825e)
cig_noreply at geodynamics.org
cig_noreply at geodynamics.org
Wed Nov 5 15:40:36 PST 2014
Repository : https://github.com/geodynamics/pylith
On branch : baagaard/dynrup-new-lagrange
Link : https://github.com/geodynamics/pylith/compare/f33c75b19fd60eedb2a3405db76a1fee333bb1d7...5b6d812b1612809fea3bd331c4e5af98c25a536a
>---------------------------------------------------------------
commit ba6825e04f1851181ee6fadc84526954d3e896b8
Author: Brad Aagaard <baagaard at usgs.gov>
Date: Fri Jun 6 14:51:38 2014 -0700
Fixed some bugs in computing residual in new FaultCohesiveDyn implementation.
>---------------------------------------------------------------
ba6825e04f1851181ee6fadc84526954d3e896b8
libsrc/pylith/faults/FaultCohesiveDyn.cc | 39 +++++++++++++++++++---
libsrc/pylith/faults/FaultCohesiveLagrange.cc | 3 +-
libsrc/pylith/friction/FrictionModel.cc | 1 -
unittests/libtests/faults/TestFaultCohesiveDyn.cc | 38 ++++++++++-----------
.../libtests/faults/data/CohesiveDynDataTri3.cc | 12 +++----
5 files changed, 60 insertions(+), 33 deletions(-)
diff --git a/libsrc/pylith/faults/FaultCohesiveDyn.cc b/libsrc/pylith/faults/FaultCohesiveDyn.cc
index fd935ab..156c4d7 100644
--- a/libsrc/pylith/faults/FaultCohesiveDyn.cc
+++ b/libsrc/pylith/faults/FaultCohesiveDyn.cc
@@ -349,6 +349,9 @@ pylith::faults::FaultCohesiveDyn::integrateResidual(const topology::Field& resid
const PetscInt rloff = residualVisitor.sectionOffset(e_lagrange);
assert(spaceDim == residualVisitor.sectionDof(e_lagrange));
+ // Get friction properties and state variables.
+ _friction->retrievePropsStateVars(v_fault);
+
#if defined(DETAILED_EVENT_LOGGING)
_logger->eventEnd(restrictEvent);
_logger->eventBegin(computeEvent);
@@ -370,10 +373,10 @@ pylith::faults::FaultCohesiveDyn::integrateResidual(const topology::Field& resid
if (fabs(slipRateVertex[iDim]) < _zeroTolerance) {
slipRateVertex[iDim] = 0.0;
} // if
+ if (fabs(slipVertex[iDim]) < _zeroTolerance) {
+ slipVertex[iDim] = 0.0;
+ } // if
} // for
- if (fabs(slipVertex[indexN]) < _zeroTolerance) {
- slipVertex[indexN] = 0.0;
- } // if
#if defined(DETAILED_EVENT_LOGGING)
_logger->eventEnd(computeEvent);
@@ -403,8 +406,8 @@ pylith::faults::FaultCohesiveDyn::integrateResidual(const topology::Field& resid
tractionContactArray[coff+iDim] = tractionRheologyVertex[iDim] - tractionInternalVertex[iDim];
for (PetscInt jDim = 0; jDim < spaceDim; ++jDim) {
- tractionRheologyGlobalVertex += orientationArray[ooff+jDim*spaceDim+iDim] * tractionRheologyVertex[jDim];
- tractionInternalGlobalVertex += orientationArray[ooff+jDim*spaceDim+iDim] * tractionInternalVertex[jDim];
+ tractionRheologyGlobalVertex[iDim] += orientationArray[ooff+jDim*spaceDim+iDim] * tractionRheologyVertex[jDim];
+ tractionInternalGlobalVertex[iDim] += orientationArray[ooff+jDim*spaceDim+iDim] * tractionInternalVertex[jDim];
} // for
} // for
@@ -416,7 +419,29 @@ pylith::faults::FaultCohesiveDyn::integrateResidual(const topology::Field& resid
const PylithScalar dispRelVertex = dispTArray[dtpoff+iDim] + dispTIncrArray[dipoff+iDim] - dispTArray[dtnoff+iDim] - dispTIncrArray[dinoff+iDim];
residualArray[rloff+iDim] += areaArray[aoff] * tractionInternalGlobalVertex[iDim] * dispRelVertex;
} // for
+
+#if 1 // DEBUGGING
+ std::cout << "v_fault: " << v_fault
+ << ", v_neg: " << v_negative
+ << ", v_pos: " << v_positive
+ << ", e_lagrange: " << e_lagrange
+ << ", area: " << areaArray[aoff]
+ << ", T_f: (";
+ for (int iDim=0; iDim < spaceDim; ++iDim) { std::cout << " " << tractionRheologyVertex[iDim]; }
+ std::cout << "), T_i: (";
+ for (int iDim=0; iDim < spaceDim; ++iDim) { std::cout << " " << tractionInternalVertex[iDim]; }
+ std::cout << "), T_fg: (";
+ for (int iDim=0; iDim < spaceDim; ++iDim) { std::cout << " " << tractionRheologyGlobalVertex[iDim]; }
+ std::cout << "), T_ig: (";
+ for (int iDim=0; iDim < spaceDim; ++iDim) { std::cout << " " << tractionInternalGlobalVertex[iDim]; }
+ std::cout << "), residual_neg: (";
+ for (int iDim=0; iDim < spaceDim; ++iDim) { std::cout << " " << residualArray[rnoff+iDim]; }
+ std::cout << "), residual_pos: (";
+ for (int iDim=0; iDim < spaceDim; ++iDim) { std::cout << " " << residualArray[rpoff+iDim]; }
+ std::cout << std::endl;
+#endif
} else { // opening, normal traction should be zero
+ std::cout << "Opening" << std::endl;
for (PetscInt iDim = 0; iDim < spaceDim; ++iDim) {
tractionContactArray[coff+iDim] = 0.0;
} // for
@@ -1080,6 +1105,8 @@ pylith::faults::FaultCohesiveDyn::_calcRheologyTraction2D(scalar_array* traction
if (fabs(slip[1]) < _zeroTolerance && tractionNormal < -_zeroTolerance) {
// if in compression and no opening
PylithScalar frictionStress = _friction->calcFriction(t, slipMag, slipRateMag, tractionNormal);
+ std::cout << " slipMag: " << slipMag << ", tractionNormal: " << tractionNormal << ", friction: " << frictionStress << std::endl;
+
#if 1 // New Newton stuff
if (tractionShearMag > 0.0 && 0.0 != jacobianShear) {
@@ -1118,6 +1145,7 @@ pylith::faults::FaultCohesiveDyn::_calcRheologyTraction2D(scalar_array* traction
} else {
(*tractionRheology)[0] = -frictionStress;
} // if/else
+ (*tractionRheology)[1] = tractionNormal;
// Determine if flipping between sliding and locked (means need new Jacobian)
if (slipMag < _zeroTolerance && frictionStress > 0.0 && tractionShearMag < _zeroTolerance) {
@@ -1200,6 +1228,7 @@ pylith::faults::FaultCohesiveDyn::_calcRheologyTraction3D(scalar_array* traction
(*tractionRheology)[0] = -frictionStress * tractionInternal[0] / tractionShearMag;
(*tractionRheology)[1] = -frictionStress * tractionInternal[1] / tractionShearMag;
} // if/else
+ (*tractionRheology)[2] = tractionNormal;
// Determine if flipping between sliding and locked (means need new Jacobian)
if (slipMag < _zeroTolerance && frictionStress > 0.0 && tractionShearMag < _zeroTolerance) {
diff --git a/libsrc/pylith/faults/FaultCohesiveLagrange.cc b/libsrc/pylith/faults/FaultCohesiveLagrange.cc
index 59855bd..44a6d2c 100644
--- a/libsrc/pylith/faults/FaultCohesiveLagrange.cc
+++ b/libsrc/pylith/faults/FaultCohesiveLagrange.cc
@@ -222,7 +222,8 @@ pylith::faults::FaultCohesiveLagrange::setupSolnDof(topology::Field* field)
const int v_negative = _cohesiveVertices[iVertex].negative;
const int e_lagrange = _cohesiveVertices[iVertex].lagrange;
- if (e_lagrange < 0) continue;
+ if (e_lagrange < 0) continue; // skipped clamped edges
+
// Set DOF in section (all subfields)
err = PetscSectionSetDof(fieldSection, e_lagrange, spaceDim);PYLITH_CHECK_ERROR(err);
err = PetscSectionSetDof(fieldSection, v_positive, spaceDim);PYLITH_CHECK_ERROR(err);
diff --git a/libsrc/pylith/friction/FrictionModel.cc b/libsrc/pylith/friction/FrictionModel.cc
index d84e565..d477de0 100644
--- a/libsrc/pylith/friction/FrictionModel.cc
+++ b/libsrc/pylith/friction/FrictionModel.cc
@@ -233,7 +233,6 @@ pylith::friction::FrictionModel::initialize(const topology::Mesh& faultMesh,
PetscScalar* stateVarArray = stateVarVisitor.localArray();
const PetscInt off = stateVarVisitor.sectionOffset(v);
const PetscInt dof = stateVarVisitor.sectionDof(v);
- std::cout << "v: " << v << ", dof: " << dof << ", stateVarsVetex: " << stateVarsVertex.size() << std::endl;
assert(stateVarsVertex.size() == dof);
for(PetscInt d = 0; d < dof; ++d, ++iOff) {
stateVarArray[off+d] += stateVarsVertex[iOff];
diff --git a/unittests/libtests/faults/TestFaultCohesiveDyn.cc b/unittests/libtests/faults/TestFaultCohesiveDyn.cc
index 1185ae9..81418c3 100644
--- a/unittests/libtests/faults/TestFaultCohesiveDyn.cc
+++ b/unittests/libtests/faults/TestFaultCohesiveDyn.cc
@@ -240,32 +240,30 @@ pylith::faults::TestFaultCohesiveDyn::testIntegrateResidualStick(void)
const topology::Field& residual = fields.get("residual");
fault.integrateResidual(residual, t, &fields);
- fault.updateStateVars(t, &fields);
-
- //residual.view("RESIDUAL"); // DEBUGGING
+ residual.view("RESIDUAL"); // DEBUGGING
{ // Check residual values
- PetscDM dmMesh = mesh.dmMesh();CPPUNIT_ASSERT(dmMesh);
- topology::Stratum verticesStratum(dmMesh, topology::Stratum::DEPTH, 0);
- const PetscInt vStart = verticesStratum.begin();
- const PetscInt vEnd = verticesStratum.end();
-
- topology::VecVisitorMesh residualVisitor(fields.get("residual"));
+ PetscInt pStart=0, pEnd=0;
+ topology::VecVisitorMesh residualVisitor(residual);
+ PetscErrorCode err = PetscSectionGetChart(residual.petscSection(), &pStart, &pEnd);CPPUNIT_ASSERT(!err);
const PetscScalar* residualArray = residualVisitor.localArray();CPPUNIT_ASSERT(residualArray);
const PylithScalar* residualE = _data->residualStickE;CPPUNIT_ASSERT(residualE);
const int fiberDimE = spaceDim; // number of values per point
const PylithScalar tolerance = 1.0e-06;
- for (PetscInt v = vStart, iVertex = 0; v < vEnd; ++v, ++iVertex) {
- const PetscInt off = residualVisitor.sectionOffset(v);
- CPPUNIT_ASSERT_EQUAL(fiberDimE, residualVisitor.sectionDof(v));
- for(PetscInt d = 0; d < fiberDimE; ++d) {
- const PylithScalar valE = residualE[iVertex*spaceDim+d];
- if (fabs(valE) > tolerance) {
- CPPUNIT_ASSERT_DOUBLES_EQUAL(1.0, residualArray[off+d]/valE, tolerance);
- } else {
- CPPUNIT_ASSERT_DOUBLES_EQUAL(valE, residualArray[off+d], tolerance);
- } // if/else
- } // for
+
+ for (PetscInt p=pStart, iPoint=0; p < pEnd; ++p, ++iPoint) {
+ if (residualVisitor.sectionDof(p) > 0) {
+ const PetscInt off = residualVisitor.sectionOffset(p);
+ CPPUNIT_ASSERT_EQUAL(fiberDimE, residualVisitor.sectionDof(p));
+ for (PetscInt iDim=0; iDim < fiberDimE; ++iDim) {
+ const PylithScalar valE = residualE[iPoint*fiberDimE+iDim];
+ if (fabs(valE) > tolerance) {
+ CPPUNIT_ASSERT_DOUBLES_EQUAL(1.0, residualArray[off+iDim]/valE, tolerance);
+ } else {
+ CPPUNIT_ASSERT_DOUBLES_EQUAL(valE, residualArray[off+iDim], tolerance);
+ } // if/else
+ } // for
+ } // if
} // for
} // Check residual values
diff --git a/unittests/libtests/faults/data/CohesiveDynDataTri3.cc b/unittests/libtests/faults/data/CohesiveDynDataTri3.cc
index 67a5f75..5f1e7de 100644
--- a/unittests/libtests/faults/data/CohesiveDynDataTri3.cc
+++ b/unittests/libtests/faults/data/CohesiveDynDataTri3.cc
@@ -282,12 +282,12 @@ const PylithScalar pylith::faults::CohesiveDynDataTri3::_fieldIncrStick[] = {
// Output
const PylithScalar pylith::faults::CohesiveDynDataTri3::_residualStickE[] = {
- 1.1, 2.1,
- 1.2, 2.2, // 3
- 1.3, 2.3, // 4
- 1.4, 2.4,
- 1.2, 2.2, // 6
- 1.3, 2.3, // 7
+ 0.0, 0.0,
+ 0.0, 0.001*(0.6*(21.6+8.6)-(9.6+2.6)), // 3
+ 0.0, 0.001*(0.6*(21.8+8.8)-(9.8+2.8)), // 4
+ 0.0, 0.0,
+ 0.0, -0.001*(0.6*(21.6+8.6)-(9.6+2.6)), // 6
+ 0.0, -0.001*(0.6*(21.8+8.8)-(9.8+2.8)), // 7
-21.6, 2.6, // 15
-21.8, 2.8, // 16
};
More information about the CIG-COMMITS
mailing list