[cig-commits] r15686 - in short/3D/PyLith/branches/pylith-friction: libsrc/faults unittests/libtests/faults unittests/pytests/utils
brad at geodynamics.org
brad at geodynamics.org
Sun Sep 20 16:59:15 PDT 2009
Author: brad
Date: 2009-09-20 16:59:15 -0700 (Sun, 20 Sep 2009)
New Revision: 15686
Modified:
short/3D/PyLith/branches/pylith-friction/libsrc/faults/FaultCohesiveDynL.cc
short/3D/PyLith/branches/pylith-friction/libsrc/faults/FaultCohesiveDynL.hh
short/3D/PyLith/branches/pylith-friction/libsrc/faults/FaultCohesiveKin.cc
short/3D/PyLith/branches/pylith-friction/libsrc/faults/FaultCohesiveKin.hh
short/3D/PyLith/branches/pylith-friction/unittests/libtests/faults/TestFaultCohesiveKin.cc
short/3D/PyLith/branches/pylith-friction/unittests/pytests/utils/TestEventLogger.py
Log:
More work on revised friction implementation.
Modified: short/3D/PyLith/branches/pylith-friction/libsrc/faults/FaultCohesiveDynL.cc
===================================================================
--- short/3D/PyLith/branches/pylith-friction/libsrc/faults/FaultCohesiveDynL.cc 2009-09-20 20:04:25 UTC (rev 15685)
+++ short/3D/PyLith/branches/pylith-friction/libsrc/faults/FaultCohesiveDynL.cc 2009-09-20 23:59:15 UTC (rev 15686)
@@ -38,6 +38,8 @@
#include <sstream> // USES std::ostringstream
#include <stdexcept> // USES std::runtime_error
+// Precomputing geometry significantly increases storage but gives a
+// slight speed improvement.
//#define PRECOMPUTE_GEOMETRY
// ----------------------------------------------------------------------
@@ -65,7 +67,7 @@
{ // deallocate
FaultCohesive::deallocate();
- // :TODO: Use shared pointers for earthquake sources
+ // :TODO: Use shared pointers for initial database
} // deallocate
// ----------------------------------------------------------------------
@@ -114,13 +116,6 @@
slip.vectorFieldType(topology::FieldBase::VECTOR);
slip.scale(_normalizer->lengthScale());
- // Allocate cumulative slip field
- _fields->add("cumulative slip", "cumulative_slip");
- topology::Field<topology::SubMesh>& cumSlip = _fields->get("cumulative slip");
- cumSlip.cloneSection(slip);
- cumSlip.vectorFieldType(topology::FieldBase::VECTOR);
- cumSlip.scale(_normalizer->lengthScale());
-
const ALE::Obj<SieveSubMesh::label_sequence>& cells =
faultSieveMesh->heightStratum(0);
assert(!cells.isNull());
@@ -134,6 +129,9 @@
// Compute orientation at vertices in fault mesh.
_calcOrientation(upDir, normalDir);
+ // Compute tributary area for each vertex in fault mesh.
+ _calcArea();
+
// Initialize quadrature geometry.
_quadrature->initializeGeometry();
@@ -141,7 +139,7 @@
_getInitialTractions();
// Setup fault constitutive model.
- //_initConstitutiveModel();
+ _initConstitutiveModel();
//logger.stagePop();
} // initialize
@@ -210,8 +208,8 @@
// vertex k make 2 contributions to the residual:
//
// * DOF i and j: internal forces in soln field associated with
- // slip
- // * DOF k: slip values
+ // slip -[C]^T{L(t)+dL(t)}
+ // * DOF k: slip values -[C]{u(t)+dt(t)}
// Get cell information and setup storage for cell data
const int spaceDim = _quadrature->spaceDim();
@@ -419,11 +417,9 @@
assert(0 != _fields);
// Cohesive cells with normal vertices i and j, and constraint
- // vertex k make 2 contributions to the residual:
+ // vertex k make contributions to the assembled residual:
//
- // * DOF i and j: internal forces in soln field associated with
- // slip
- // * DOF k: slip values
+ // * DOF k: slip values {D(t+dt)}
topology::Field<topology::SubMesh>& slip = _fields->get("slip");
slip.zero();
@@ -764,7 +760,7 @@
_allocateBufferVertexVectorField();
topology::Field<topology::SubMesh>& buffer =
_fields->get("buffer (vector)");
- //_calcTractionsChange(&buffer, dispT);
+ _calcTractionsChange(&buffer, dispT);
return buffer;
} else {
@@ -991,8 +987,102 @@
//orientation.view("ORIENTATION");
} // _calcOrientation
-#if 0
// ----------------------------------------------------------------------
+void
+pylith::faults::FaultCohesiveDynL::_calcArea(void)
+{ // _calcArea
+ assert(0 != _faultMesh);
+ assert(0 != _fields);
+
+ // Containers for area information
+ const int cellDim = _quadrature->cellDim();
+ const int numBasis = _quadrature->numBasis();
+ const int numQuadPts = _quadrature->numQuadPts();
+ const int spaceDim = _quadrature->spaceDim();
+ const feassemble::CellGeometry& cellGeometry = _quadrature->refGeometry();
+ const double_array& quadWts = _quadrature->quadWts();
+ assert(quadWts.size() == numQuadPts);
+ double jacobianDet = 0;
+ double_array areaCell(numBasis);
+
+ // Get vertices in fault mesh.
+ const ALE::Obj<SieveSubMesh>& faultSieveMesh = _faultMesh->sieveMesh();
+ assert(!faultSieveMesh.isNull());
+ const ALE::Obj<SieveSubMesh::label_sequence>& vertices =
+ faultSieveMesh->depthStratum(0);
+ const SieveSubMesh::label_sequence::iterator verticesBegin = vertices->begin();
+ const SieveSubMesh::label_sequence::iterator verticesEnd = vertices->end();
+
+ // Allocate area field.
+ _fields->add("area", "area");
+
+ topology::Field<topology::SubMesh>& area = _fields->get("area");
+ const topology::Field<topology::SubMesh>& slip = _fields->get("slip");
+ area.newSection(slip, 1);
+ area.allocate();
+ area.zero();
+ const ALE::Obj<RealSection>& areaSection = area.section();
+ assert(!areaSection.isNull());
+ topology::Mesh::UpdateAddVisitor areaVisitor(*areaSection, &areaCell[0]);
+
+ double_array coordinatesCell(numBasis*spaceDim);
+ const ALE::Obj<RealSection>& coordinates =
+ faultSieveMesh->getRealSection("coordinates");
+ assert(!coordinates.isNull());
+ topology::Mesh::RestrictVisitor coordsVisitor(*coordinates,
+ coordinatesCell.size(),
+ &coordinatesCell[0]);
+
+ const ALE::Obj<SieveSubMesh::label_sequence>& cells =
+ faultSieveMesh->heightStratum(0);
+ assert(!cells.isNull());
+ const SieveSubMesh::label_sequence::iterator cellsBegin = cells->begin();
+ const SieveSubMesh::label_sequence::iterator cellsEnd = cells->end();
+
+ // Loop over cells in fault mesh, compute area
+ for (SieveSubMesh::label_sequence::iterator c_iter = cellsBegin;
+ c_iter != cellsEnd;
+ ++c_iter) {
+ areaCell = 0.0;
+
+ // Compute geometry information for current cell
+#if defined(PRECOMPUTE_GEOMETRY)
+ _quadrature->retrieveGeometry(*c_iter);
+#else
+ coordsVisitor.clear();
+ faultSieveMesh->restrictClosure(*c_iter, coordsVisitor);
+ _quadrature->computeGeometry(coordinatesCell, *c_iter);
+#endif
+
+ // Get cell geometry information that depends on cell
+ const double_array& basis = _quadrature->basis();
+ const double_array& jacobianDet = _quadrature->jacobianDet();
+
+ // Compute area
+ for (int iQuad=0; iQuad < numQuadPts; ++iQuad) {
+ const double wt = quadWts[iQuad] * jacobianDet[iQuad];
+ for (int iBasis=0; iBasis < numBasis; ++iBasis) {
+ const double dArea = wt*basis[iQuad*numBasis+iBasis];
+ areaCell[iBasis] += dArea;
+ } // for
+ } // for
+ areaVisitor.clear();
+ faultSieveMesh->updateAdd(*c_iter, areaVisitor);
+
+ PetscLogFlops( numQuadPts*(1+numBasis*2) );
+ } // for
+
+ // Assemble area information
+ area.complete();
+
+#if 0 // DEBUGGING
+ area.view("AREA");
+ //_faultMesh->getSendOverlap()->view("Send fault overlap");
+ //_faultMesh->getRecvOverlap()->view("Receive fault overlap");
+#endif
+} // _calcArea
+
+// ----------------------------------------------------------------------
// Compute change in tractions on fault surface using solution.
// NOTE: We must convert vertex labels to fault vertex labels
void
@@ -1108,7 +1198,6 @@
tractions->view("TRACTIONS");
#endif
} // _calcTractionsChange
-#endif
// ----------------------------------------------------------------------
void
Modified: short/3D/PyLith/branches/pylith-friction/libsrc/faults/FaultCohesiveDynL.hh
===================================================================
--- short/3D/PyLith/branches/pylith-friction/libsrc/faults/FaultCohesiveDynL.hh 2009-09-20 20:04:25 UTC (rev 15685)
+++ short/3D/PyLith/branches/pylith-friction/libsrc/faults/FaultCohesiveDynL.hh 2009-09-20 23:59:15 UTC (rev 15686)
@@ -67,10 +67,25 @@
* {r(t+dt)} = {b(t+dt)} - [A(t+dt) C^T ]{ u(t)+du(t) }
* {D(t+dt)} [ C 0 ]{ L(t)+dL(t) }
*
- * The term D does not involve integration over cohesive cells. We
- * integrate the Lagrange multiplier terms over the cohesive cells
- * because this introduces weighting of the orientation of the fault
- * for the direction of slip at the vertices of the cohesive cells.
+ * The terms in the residual contributing to the DOF at the Lagrange
+ * vertices are
+ *
+ * {r(t+dt)} = {D(t+dt)} - [C]{u(t)+dt(t)}
+ *
+ * The first term, {D(t+dt)}, does not involve integration over the
+ * cohesive cells, so it does not require assembling over cohesive
+ * cells or processors. We compute the term in
+ * integrateResidualAssembled().
+ *
+ * The term in the residual contributing to the DOF at the
+ * non-Lagrange vertices of the cohesive cells is
+ *
+ * {r(t+dt)} = -[C]^T{L(t)+dL(t)}
+ *
+ * We integrate the Lagrange multiplier term and the relative
+ * displacement term over the cohesive cells, because this introduces
+ * weighting of the orientation of the fault for the direction of slip
+ * at the vertices of the cohesive cells.
*/
class pylith::faults::FaultCohesiveDynL : public FaultCohesive
{ // class FaultCohesiveDynL
@@ -111,7 +126,8 @@
const double upDir[3],
const double normalDir[3]);
- /** Split solution field for separate preconditioning.
+ /** Split solution field for separate preconditioning of normal DOF
+ * from DOF associated with Lagrange multipliers.
*
* @param field Solution field.
*/
@@ -216,6 +232,9 @@
void _calcOrientation(const double upDir[3],
const double normalDir[3]);
+ /// Calculate fault area field.
+ void _calcArea(void);
+
/** Get initial tractions using a spatial database.
*/
void _getInitialTractions(void);
@@ -224,6 +243,15 @@
*/
void _initConstitutiveModel(void);
+ /** Compute change in tractions on fault surface using solution.
+ *
+ * @param tractions Field for tractions.
+ * @param mesh Finite-element mesh for domain
+ * @param solution Solution over domain
+ */
+ void _calcTractionsChange(topology::Field<topology::SubMesh>* tractions,
+ const topology::Field<topology::Mesh>& solution);
+
/// Allocate buffer for vector field.
void _allocateBufferVertexVectorField(void);
Modified: short/3D/PyLith/branches/pylith-friction/libsrc/faults/FaultCohesiveKin.cc
===================================================================
--- short/3D/PyLith/branches/pylith-friction/libsrc/faults/FaultCohesiveKin.cc 2009-09-20 20:04:25 UTC (rev 15685)
+++ short/3D/PyLith/branches/pylith-friction/libsrc/faults/FaultCohesiveKin.cc 2009-09-20 23:59:15 UTC (rev 15686)
@@ -135,13 +135,6 @@
slip.vectorFieldType(topology::FieldBase::VECTOR);
slip.scale(_normalizer->lengthScale());
- // Allocate cumulative slip field
- _fields->add("cumulative slip", "cumulative_slip");
- topology::Field<topology::SubMesh>& cumSlip = _fields->get("cumulative slip");
- cumSlip.cloneSection(slip);
- cumSlip.vectorFieldType(topology::FieldBase::VECTOR);
- cumSlip.scale(_normalizer->lengthScale());
-
const ALE::Obj<SieveSubMesh::label_sequence>& cells =
faultSieveMesh->heightStratum(0);
assert(!cells.isNull());
@@ -225,8 +218,8 @@
// vertex k make 2 contributions to the residual:
//
// * DOF i and j: internal forces in soln field associated with
- // slip
- // * DOF k: slip values
+ // slip -[C]^T{L(t)+dL(t)}
+ // * DOF k: slip values -[C]{u(t)+dt(t)}
// Get cell information and setup storage for cell data
const int spaceDim = _quadrature->spaceDim();
@@ -434,11 +427,9 @@
assert(0 != _fields);
// Cohesive cells with normal vertices i and j, and constraint
- // vertex k make 2 contributions to the residual:
+ // vertex k make contributions to the assembled residual:
//
- // * DOF i and j: internal forces in soln field associated with
- // slip
- // * DOF k: slip values
+ // * DOF k: slip values {D(t+dt)}
topology::Field<topology::SubMesh>& slip = _fields->get("slip");
slip.zero();
@@ -647,12 +638,6 @@
assert(0 != fields);
assert(0 != _fields);
- // Update cumulative slip
- topology::Field<topology::SubMesh>& cumSlip = _fields->get("cumulative slip");
- topology::Field<topology::SubMesh>& slip = _fields->get("slip");
- if (!_useSolnIncr)
- cumSlip.zero();
- cumSlip += slip;
} // updateStateVars
// ----------------------------------------------------------------------
@@ -727,9 +712,9 @@
double scale = 0.0;
int fiberDim = 0;
if (0 == strcasecmp("slip", name)) {
- const topology::Field<topology::SubMesh>& cumSlip =
- _fields->get("cumulative slip");
- return cumSlip;
+ const topology::Field<topology::SubMesh>& slip =
+ _fields->get("slip");
+ return slip;
} else if (cohesiveDim > 0 && 0 == strcasecmp("strike_dir", name)) {
const ALE::Obj<RealSection>& orientationSection =
@@ -1247,13 +1232,13 @@
ALE::MemoryLogger& logger = ALE::MemoryLogger::singleton();
logger.stagePush("Output");
- // Create vector field; use same shape/chart as cumulative slip field.
+ // Create vector field; use same shape/chart as slip field.
assert(0 != _faultMesh);
_fields->add("buffer (vector)", "buffer");
topology::Field<topology::SubMesh>& buffer =
_fields->get("buffer (vector)");
const topology::Field<topology::SubMesh>& slip =
- _fields->get("cumulative slip");
+ _fields->get("slip");
buffer.cloneSection(slip);
buffer.zero();
Modified: short/3D/PyLith/branches/pylith-friction/libsrc/faults/FaultCohesiveKin.hh
===================================================================
--- short/3D/PyLith/branches/pylith-friction/libsrc/faults/FaultCohesiveKin.hh 2009-09-20 20:04:25 UTC (rev 15685)
+++ short/3D/PyLith/branches/pylith-friction/libsrc/faults/FaultCohesiveKin.hh 2009-09-20 23:59:15 UTC (rev 15686)
@@ -67,10 +67,25 @@
* {r(t+dt)} = {b(t+dt)} - [A(t+dt) C^T ]{ u(t)+du(t) }
* {D(t+dt)} [ C 0 ]{ L(t)+dL(t) }
*
- * The term D does not involve integration over cohesive cells. We
- * integrate the Lagrange multiplier terms over the cohesive cells
- * because this introduces weighting of the orientation of the fault
- * for the direction of slip at the vertices of the cohesive cells.
+ * The terms in the residual contributing to the DOF at the Lagrange
+ * vertices are
+ *
+ * {r(t+dt)} = {D(t+dt)} - [C]{u(t)+dt(t)}
+ *
+ * The first term, {D(t+dt)}, does not involve integration over the
+ * cohesive cells, so it does not require assembling over cohesive
+ * cells or processors. We compute the term in
+ * integrateResidualAssembled().
+ *
+ * The term in the residual contributing to the DOF at the
+ * non-Lagrange vertices of the cohesive cells is
+ *
+ * {r(t+dt)} = -[C]^T{L(t)+dL(t)}
+ *
+ * We integrate the Lagrange multiplier term and the relative
+ * displacement term over the cohesive cells, because this introduces
+ * weighting of the orientation of the fault for the direction of slip
+ * at the vertices of the cohesive cells.
*/
class pylith::faults::FaultCohesiveKin : public FaultCohesive
{ // class FaultCohesiveKin
Modified: short/3D/PyLith/branches/pylith-friction/unittests/libtests/faults/TestFaultCohesiveKin.cc
===================================================================
--- short/3D/PyLith/branches/pylith-friction/unittests/libtests/faults/TestFaultCohesiveKin.cc 2009-09-20 20:04:25 UTC (rev 15685)
+++ short/3D/PyLith/branches/pylith-friction/unittests/libtests/faults/TestFaultCohesiveKin.cc 2009-09-20 23:59:15 UTC (rev 15686)
@@ -622,16 +622,16 @@
const SieveSubMesh::label_sequence::iterator verticesEnd = vertices->end();
SieveSubMesh::renumbering_type& renumbering = faultSieveMesh->getRenumbering();
- // Compute expected cumulative slip using eqsrcs
- topology::Field<topology::SubMesh> cumSlipE(*fault._faultMesh);
- cumSlipE.newSection(topology::FieldBase::VERTICES_FIELD, spaceDim);
- cumSlipE.allocate();
- const ALE::Obj<RealSection> cumSlipESection = cumSlipE.section();
- CPPUNIT_ASSERT(!cumSlipESection.isNull());
+ // Compute expected slip using eqsrcs
+ topology::Field<topology::SubMesh> slipE(*fault._faultMesh);
+ slipE.newSection(topology::FieldBase::VERTICES_FIELD, spaceDim);
+ slipE.allocate();
+ const ALE::Obj<RealSection> slipESection = slipE.section();
+ CPPUNIT_ASSERT(!slipESection.isNull());
- const ALE::Obj<RealSection> cumSlipSection =
- fault._fields->get("cumulative slip").section();
- CPPUNIT_ASSERT(!cumSlipSection.isNull());
+ const ALE::Obj<RealSection> slipSection =
+ fault._fields->get("slip").section();
+ CPPUNIT_ASSERT(!slipSection.isNull());
const FaultCohesiveKin::srcs_type::const_iterator srcsEnd = fault._eqSrcs.end();
for (FaultCohesiveKin::srcs_type::iterator s_iter=fault._eqSrcs.begin();
@@ -640,7 +640,7 @@
EqKinSrc* src = s_iter->second;
assert(0 != src);
if (t >= src->originTime())
- src->slip(&cumSlipE, t);
+ src->slip(&slipE, t);
} // for
int iVertex = 0;
@@ -660,13 +660,13 @@
} // for
CPPUNIT_ASSERT(found);
- // Check _cumSlip
- int fiberDim = cumSlipSection->getFiberDimension(*v_iter);
+ // Check _slip
+ int fiberDim = slipSection->getFiberDimension(*v_iter);
CPPUNIT_ASSERT_EQUAL(spaceDim, fiberDim);
- const double* slipV = cumSlipSection->restrictPoint(*v_iter);
+ const double* slipV = slipSection->restrictPoint(*v_iter);
CPPUNIT_ASSERT(0 != slipV);
- const double* slipE = cumSlipESection->restrictPoint(*v_iter);
+ const double* slipE = slipESection->restrictPoint(*v_iter);
CPPUNIT_ASSERT(0 != slipE);
for (int iDim=0; iDim < spaceDim; ++iDim) {
Modified: short/3D/PyLith/branches/pylith-friction/unittests/pytests/utils/TestEventLogger.py
===================================================================
--- short/3D/PyLith/branches/pylith-friction/unittests/pytests/utils/TestEventLogger.py 2009-09-20 20:04:25 UTC (rev 15685)
+++ short/3D/PyLith/branches/pylith-friction/unittests/pytests/utils/TestEventLogger.py 2009-09-20 23:59:15 UTC (rev 15686)
@@ -108,7 +108,7 @@
logger.className("logging A")
logger.initialize()
- stages = ["stage 1" , "stage 2" , "stage 3"]
+ stages = ["stage 1a" , "stage 2a" , "stage 3a"]
id = {}
for stage in stages:
id[stage] = logger.registerStage(stage)
@@ -125,18 +125,18 @@
logger = EventLogger()
logger.className("logging A")
logger.initialize()
- stages = ["stage 1" , "stage 2" , "stage 3"]
+ stages = ["stage 1b" , "stage 2b" , "stage 3b"]
for stage in stages:
logger.registerStage(stage)
- logger.stagePush("stage 2")
+ logger.stagePush("stage 2b")
logger.stagePop()
- logger.stagePush("stage 1")
+ logger.stagePush("stage 1b")
logger.stagePop()
- logger.stagePush("stage 3")
- logger.stagePush("stage 1")
+ logger.stagePush("stage 3b")
+ logger.stagePush("stage 1b")
logger.stagePop()
logger.stagePop()
return
More information about the CIG-COMMITS
mailing list