[cig-commits] r19617 - in short/3D/PyLith/trunk: examples/3d/hex8 examples/3d/tet4 libsrc/pylith/faults libsrc/pylith/feassemble libsrc/pylith/friction libsrc/pylith/materials libsrc/pylith/meshio libsrc/pylith/topology pylith/bc pylith/faults pylith/meshio pylith/perf pylith/topology pylith/utils tests/topology unittests/libtests/topology
brad at geodynamics.org
brad at geodynamics.org
Fri Feb 10 21:38:06 PST 2012
Author: brad
Date: 2012-02-10 21:38:06 -0800 (Fri, 10 Feb 2012)
New Revision: 19617
Modified:
short/3D/PyLith/trunk/examples/3d/hex8/pylithapp.cfg
short/3D/PyLith/trunk/examples/3d/hex8/test.cfg
short/3D/PyLith/trunk/examples/3d/tet4/pylithapp.cfg
short/3D/PyLith/trunk/libsrc/pylith/faults/CohesiveTopology.cc
short/3D/PyLith/trunk/libsrc/pylith/faults/Fault.cc
short/3D/PyLith/trunk/libsrc/pylith/faults/Fault.hh
short/3D/PyLith/trunk/libsrc/pylith/faults/Fault.icc
short/3D/PyLith/trunk/libsrc/pylith/faults/FaultCohesive.cc
short/3D/PyLith/trunk/libsrc/pylith/faults/FaultCohesive.hh
short/3D/PyLith/trunk/libsrc/pylith/faults/FaultCohesiveDyn.cc
short/3D/PyLith/trunk/libsrc/pylith/faults/FaultCohesiveDyn.hh
short/3D/PyLith/trunk/libsrc/pylith/faults/FaultCohesiveLagrange.cc
short/3D/PyLith/trunk/libsrc/pylith/feassemble/IntegratorElasticity.cc
short/3D/PyLith/trunk/libsrc/pylith/friction/FrictionModel.cc
short/3D/PyLith/trunk/libsrc/pylith/friction/FrictionModel.hh
short/3D/PyLith/trunk/libsrc/pylith/friction/RateStateAgeing.cc
short/3D/PyLith/trunk/libsrc/pylith/friction/RateStateAgeing.hh
short/3D/PyLith/trunk/libsrc/pylith/friction/SlipWeakening.cc
short/3D/PyLith/trunk/libsrc/pylith/friction/SlipWeakening.hh
short/3D/PyLith/trunk/libsrc/pylith/friction/SlipWeakeningTime.cc
short/3D/PyLith/trunk/libsrc/pylith/friction/SlipWeakeningTime.hh
short/3D/PyLith/trunk/libsrc/pylith/friction/StaticFriction.cc
short/3D/PyLith/trunk/libsrc/pylith/friction/StaticFriction.hh
short/3D/PyLith/trunk/libsrc/pylith/friction/TimeWeakening.cc
short/3D/PyLith/trunk/libsrc/pylith/friction/TimeWeakening.hh
short/3D/PyLith/trunk/libsrc/pylith/materials/ElasticMaterial.cc
short/3D/PyLith/trunk/libsrc/pylith/materials/Material.cc
short/3D/PyLith/trunk/libsrc/pylith/meshio/CellFilterAvg.cc
short/3D/PyLith/trunk/libsrc/pylith/meshio/DataWriterHDF5.cc
short/3D/PyLith/trunk/libsrc/pylith/meshio/DataWriterHDF5Ext.cc
short/3D/PyLith/trunk/libsrc/pylith/meshio/DataWriterVTK.cc
short/3D/PyLith/trunk/libsrc/pylith/meshio/MeshBuilder.cc
short/3D/PyLith/trunk/libsrc/pylith/meshio/VertexFilterVecNorm.cc
short/3D/PyLith/trunk/libsrc/pylith/topology/MeshOrder.cc
short/3D/PyLith/trunk/libsrc/pylith/topology/SubMesh.hh
short/3D/PyLith/trunk/libsrc/pylith/topology/SubMesh.icc
short/3D/PyLith/trunk/pylith/bc/AbsorbingDampers.py
short/3D/PyLith/trunk/pylith/bc/DirichletBC.py
short/3D/PyLith/trunk/pylith/bc/DirichletBoundary.py
short/3D/PyLith/trunk/pylith/bc/Neumann.py
short/3D/PyLith/trunk/pylith/bc/PointForce.py
short/3D/PyLith/trunk/pylith/faults/FaultCohesiveDyn.py
short/3D/PyLith/trunk/pylith/faults/FaultCohesiveKin.py
short/3D/PyLith/trunk/pylith/meshio/OutputFaultDyn.py
short/3D/PyLith/trunk/pylith/meshio/OutputFaultKin.py
short/3D/PyLith/trunk/pylith/perf/Fault.py
short/3D/PyLith/trunk/pylith/perf/MemoryLogger.py
short/3D/PyLith/trunk/pylith/perf/Mesh.py
short/3D/PyLith/trunk/pylith/topology/Mesh.py
short/3D/PyLith/trunk/pylith/topology/MeshGenerator.py
short/3D/PyLith/trunk/pylith/utils/profiling.py
short/3D/PyLith/trunk/tests/topology/test_meshmem.py
short/3D/PyLith/trunk/unittests/libtests/topology/TestRefineUniform.cc
Log:
Merge from stable.
Modified: short/3D/PyLith/trunk/examples/3d/hex8/pylithapp.cfg
===================================================================
--- short/3D/PyLith/trunk/examples/3d/hex8/pylithapp.cfg 2012-02-11 05:37:58 UTC (rev 19616)
+++ short/3D/PyLith/trunk/examples/3d/hex8/pylithapp.cfg 2012-02-11 05:38:06 UTC (rev 19617)
@@ -89,7 +89,7 @@
# Linear solver monitoring options.
ksp_monitor = true
-ksp_view = true
+#ksp_view = true
ksp_converged_reason = true
# Nonlinear solver monitoring options.
@@ -98,7 +98,7 @@
snes_max_it = 100
snes_monitor = true
snes_ls_monitor = true
-snes_view = true
+#snes_view = true
snes_converged_reason = true
# PETSc summary -- useful for performance information.
Modified: short/3D/PyLith/trunk/examples/3d/hex8/test.cfg
===================================================================
--- short/3D/PyLith/trunk/examples/3d/hex8/test.cfg 2012-02-11 05:37:58 UTC (rev 19616)
+++ short/3D/PyLith/trunk/examples/3d/hex8/test.cfg 2012-02-11 05:38:06 UTC (rev 19617)
@@ -119,12 +119,12 @@
[pylithapp.timedependent.interfaces.fault]
# The label corresponds to the name of the nodeset in CUBIT.
label = fault
-zero_tolerance = 1.0e-15
+zero_tolerance = 1.0e-12
# Use the rate-and-state aging friction model.
friction = pylith.friction.RateStateAgeing
friction.label = Rate and state
-friction.min_slip_rate = 1.0e-10
+friction.min_slip_rate = 1.0e-9
# We must define the quadrature information for fault cells.
# The fault cells are 2D (surface).
@@ -142,7 +142,7 @@
friction.db_properties = spatialdata.spatialdb.UniformDB
friction.db_properties.label = Rate Stete Ageing
friction.db_properties.values = [reference-friction-coefficient,reference-slip-rate,characteristic-slip-distance,constitutive-parameter-a,constitutive-parameter-b,cohesion]
-friction.db_properties.data = [0.4, 1.0e-9*m/s, 0.05*m, 0.01, 0.08, 0.0*Pa]
+friction.db_properties.data = [0.4, 1.0e-7*m/s, 2.0*m, 0.01, 0.02, 0.0*Pa]
# Set spatial database for the initial value of the state variable.
friction.db_initial_state = spatialdata.spatialdb.UniformDB
@@ -162,7 +162,7 @@
# fault constitutive model.
friction_pc_type = asm
friction_sub_pc_factor_shift_type = nonzero
-friction_ksp_max_it = 25
+friction_ksp_max_it = 50
friction_ksp_gmres_restart = 30
# Uncomment to view details of friction sensitivity solve.
#friction_ksp_monitor = true
@@ -170,15 +170,15 @@
friction_ksp_converged_reason = true
# Reduce convergence tolerances.
-ksp_rtol = 1.0e-18
-ksp_atol = 1.0e-18
+ksp_rtol = 1.0e-16
+ksp_atol = 1.0e-15
-snes_rtol = 1.0e-16
-snes_atol = 1.0e-16
+snes_rtol = 1.0e-14
+snes_atol = 1.0e-12
snes_max_it = 200
+info = 1
-
# ----------------------------------------------------------------------
# output
# ----------------------------------------------------------------------
Modified: short/3D/PyLith/trunk/examples/3d/tet4/pylithapp.cfg
===================================================================
--- short/3D/PyLith/trunk/examples/3d/tet4/pylithapp.cfg 2012-02-11 05:37:58 UTC (rev 19616)
+++ short/3D/PyLith/trunk/examples/3d/tet4/pylithapp.cfg 2012-02-11 05:38:06 UTC (rev 19617)
@@ -18,6 +18,11 @@
meshimporter = 1
#quadrature3d = 1
+[pylithapp.journal.debug]
+pylithapp = 1
+problem = 1
+implicit = 1
+
# ----------------------------------------------------------------------
# mesh_generator
# ----------------------------------------------------------------------
Modified: short/3D/PyLith/trunk/libsrc/pylith/faults/CohesiveTopology.cc
===================================================================
--- short/3D/PyLith/trunk/libsrc/pylith/faults/CohesiveTopology.cc 2012-02-11 05:37:58 UTC (rev 19616)
+++ short/3D/PyLith/trunk/libsrc/pylith/faults/CohesiveTopology.cc 2012-02-11 05:38:06 UTC (rev 19617)
@@ -51,7 +51,7 @@
// Memory logging
ALE::MemoryLogger& logger = ALE::MemoryLogger::singleton();
- logger.stagePush("FaultCreation");
+ logger.stagePush("SerialFaultCreation");
faultMesh->coordsys(mesh);
@@ -98,10 +98,10 @@
fault->setSieve(faultSieve);
logger.stagePop();
- logger.stagePush("FaultStratification");
+ logger.stagePush("SerialFaultStratification");
fault->stratify();
logger.stagePop();
- logger.stagePush("FaultCreation");
+ logger.stagePush("SerialFaultCreation");
if (debug)
fault->view("Fault mesh");
@@ -143,7 +143,7 @@
// Memory logging
ALE::MemoryLogger& logger = ALE::MemoryLogger::singleton();
- logger.stagePush("FaultCreation");
+ logger.stagePush("SerialFaultCreation");
const ALE::Obj<SieveMesh>& sieveMesh = mesh->sieveMesh();
assert(!sieveMesh.isNull());
@@ -212,7 +212,7 @@
<< vertexRenumber[*v_iter] << std::endl;
logger.stagePop();
- logger.stagePush("FaultStratification");
+ logger.stagePush("SerialFaultStratification");
// Add shadow and constraint vertices (if they exist) to group
// associated with fault
groupField->addPoint(firstFaultVertex, 1);
@@ -232,7 +232,7 @@
++firstLagrangeVertex;
} // if
logger.stagePop();
- logger.stagePush("FaultCreation");
+ logger.stagePush("SerialFaultCreation");
// Add shadow vertices to other groups, don't add constraint
// vertices (if they exist) because we don't want BC, etc to act
@@ -247,12 +247,16 @@
group->addPoint(firstFaultVertex, 1);
} // for
} // for
+ logger.stagePop();
+ logger.stagePush("SerialFaultCreation");
const std::set<std::string>::const_iterator namesEnd = groupNames->end();
for(std::set<std::string>::const_iterator name = groupNames->begin();
name != namesEnd;
++name) {
sieveMesh->reallocate(sieveMesh->getIntSection(*name));
} // for
+ logger.stagePop();
+ logger.stagePush("SerialFault");
// Split the mesh along the fault sieve and create cohesive elements
const ALE::Obj<SieveSubMesh::label_sequence>& faces =
@@ -410,14 +414,14 @@
// TODO: Need to reform the material label when sieve is reallocated
sieveMesh->setValue(material, firstFaultCell, materialId);
logger.stagePop();
- logger.stagePush("FaultStratification");
+ logger.stagePush("SerialFaultStratification");
#if defined(FAST_STRATIFY)
// OPTIMIZATION
sieveMesh->setHeight(firstFaultCell, 0);
sieveMesh->setDepth(firstFaultCell, 1);
#endif
logger.stagePop();
- logger.stagePush("FaultCreation");
+ logger.stagePush("SerialFaultCreation");
sV2.clear();
cV2.clear();
} // for
@@ -606,10 +610,10 @@
delete [] indices;
#if !defined(FAST_STRATIFY)
logger.stagePop();
- logger.stagePush("FaultStratification");
+ logger.stagePush("SerialFaultStratification");
sieveMesh->stratify();
logger.stagePop();
- logger.stagePush("FaultCreation");
+ logger.stagePush("SerialFaultCreation");
#endif
const std::string labelName("censored depth");
Modified: short/3D/PyLith/trunk/libsrc/pylith/faults/Fault.cc
===================================================================
--- short/3D/PyLith/trunk/libsrc/pylith/faults/Fault.cc 2012-02-11 05:37:58 UTC (rev 19616)
+++ short/3D/PyLith/trunk/libsrc/pylith/faults/Fault.cc 2012-02-11 05:38:06 UTC (rev 19617)
@@ -54,4 +54,44 @@
return *_faultMesh;
} // faultMesh
+
+// ----------------------------------------------------------------------
+// Get dimension of mesh.
+int
+pylith::faults::Fault::dimension(void) const
+{ // dimension
+ return (_faultMesh) ? _faultMesh->dimension() : 0;
+} // dimension
+
+
+// ----------------------------------------------------------------------
+// Get representative cone size for mesh.
+int
+pylith::faults::Fault::coneSize(void) const
+{ // coneSize
+
+ return (_faultMesh && numCells() > 0) ?
+ _faultMesh->sieveMesh()->getSieve()->getConeSize(*_faultMesh->sieveMesh()->heightStratum(1)->begin()) : 0;
+} // coneSize
+
+
+// ----------------------------------------------------------------------
+// Get number of vertices in mesh.
+int
+pylith::faults::Fault::numVertices(void) const
+{ // numVertices
+ return (_faultMesh) ? _faultMesh->numVertices() : 0;
+} // numVertices
+
+
+// ----------------------------------------------------------------------
+// Get number of cells in mesh.
+int
+pylith::faults::Fault::numCells(void) const
+{ // numCells
+ return (_faultMesh && !_faultMesh->sieveMesh().isNull() && _faultMesh->sieveMesh()->height() > 0) ?
+ _faultMesh->sieveMesh()->heightStratum(1)->size() : 0;
+} // numCells
+
+
// End of file
Modified: short/3D/PyLith/trunk/libsrc/pylith/faults/Fault.hh
===================================================================
--- short/3D/PyLith/trunk/libsrc/pylith/faults/Fault.hh 2012-02-11 05:37:58 UTC (rev 19616)
+++ short/3D/PyLith/trunk/libsrc/pylith/faults/Fault.hh 2012-02-11 05:38:06 UTC (rev 19617)
@@ -85,13 +85,38 @@
*/
const char* label(void) const;
- /** Get the number of vertices on the fault.
+ /** Get dimension of mesh.
*
+ * @returns Dimension of mesh.
+ */
+ int dimension(void) const;
+
+ /** Get representative cone size for mesh.
+ *
+ * @returns Representative cone size for mesh.
+ */
+ int coneSize(void) const;
+
+ /** Get number of vertices in mesh.
+ *
+ * @returns Number of vertices in mesh.
+ */
+ int numVertices(void) const;
+
+ /** Get number of cells in mesh.
+ *
+ * @returns Number of cells in mesh.
+ */
+ int numCells(void) const;
+
+ /** Get the number of vertices associated with the fault (before
+ * fault mesh exists).
+ *
* @param mesh PETSc mesh
* @return Number of vertices on the fault.
*/
virtual
- int numVertices(const topology::Mesh& mesh) const = 0;
+ int numVerticesNoMesh(const topology::Mesh& mesh) const = 0;
/** Adjust mesh topology for fault implementation.
*
Modified: short/3D/PyLith/trunk/libsrc/pylith/faults/Fault.icc
===================================================================
--- short/3D/PyLith/trunk/libsrc/pylith/faults/Fault.icc 2012-02-11 05:37:58 UTC (rev 19616)
+++ short/3D/PyLith/trunk/libsrc/pylith/faults/Fault.icc 2012-02-11 05:38:06 UTC (rev 19617)
@@ -48,5 +48,4 @@
return _label.c_str();
}
-
// End of file
Modified: short/3D/PyLith/trunk/libsrc/pylith/faults/FaultCohesive.cc
===================================================================
--- short/3D/PyLith/trunk/libsrc/pylith/faults/FaultCohesive.cc 2012-02-11 05:37:58 UTC (rev 19616)
+++ short/3D/PyLith/trunk/libsrc/pylith/faults/FaultCohesive.cc 2012-02-11 05:38:06 UTC (rev 19617)
@@ -78,10 +78,11 @@
} // faultMeshFilename
// ----------------------------------------------------------------------
-// Get number of vertices in fault.
+// Get the number of vertices associated with the fault (before
+// fault mesh exists).
int
-pylith::faults::FaultCohesive::numVertices(const topology::Mesh& mesh) const
-{ // numVertices
+pylith::faults::FaultCohesive::numVerticesNoMesh(const topology::Mesh& mesh) const
+{ // numVerticesNoMesh
int nvertices = 0;
if (!_useFaultMesh) {
@@ -105,7 +106,7 @@
} // else
return nvertices;
-} // numVertices
+} // numVerticesNoMesh
// ----------------------------------------------------------------------
// Adjust mesh topology for fault implementation.
Modified: short/3D/PyLith/trunk/libsrc/pylith/faults/FaultCohesive.hh
===================================================================
--- short/3D/PyLith/trunk/libsrc/pylith/faults/FaultCohesive.hh 2012-02-11 05:37:58 UTC (rev 19616)
+++ short/3D/PyLith/trunk/libsrc/pylith/faults/FaultCohesive.hh 2012-02-11 05:38:06 UTC (rev 19617)
@@ -74,12 +74,13 @@
*/
void faultMeshFilename(const char* filename);
- /** Get the number of vertices on the fault.
+ /** Get the number of vertices associated with the fault (before
+ * fault mesh exists).
*
* @param mesh PETSc mesh
* @return Number of vertices on the fault.
*/
- int numVertices(const topology::Mesh& mesh) const;
+ int numVerticesNoMesh(const topology::Mesh& mesh) const;
/** Adjust mesh topology for fault implementation.
*
Modified: short/3D/PyLith/trunk/libsrc/pylith/faults/FaultCohesiveDyn.cc
===================================================================
--- short/3D/PyLith/trunk/libsrc/pylith/faults/FaultCohesiveDyn.cc 2012-02-11 05:37:58 UTC (rev 19616)
+++ short/3D/PyLith/trunk/libsrc/pylith/faults/FaultCohesiveDyn.cc 2012-02-11 05:38:06 UTC (rev 19617)
@@ -149,6 +149,9 @@
const spatialdata::geocoords::CoordSys* cs = mesh.coordsys();
assert(0 != cs);
+ ALE::MemoryLogger& logger = ALE::MemoryLogger::singleton();
+ logger.stagePush("FaultFields");
+
// Create field for relative velocity associated with Lagrange vertex k
_fields->add("relative velocity", "relative_velocity");
topology::Field<topology::SubMesh>& velRel =
@@ -158,7 +161,7 @@
velRel.vectorFieldType(topology::FieldBase::VECTOR);
velRel.scale(_normalizer->lengthScale() / _normalizer->timeScale());
- //logger.stagePop();
+ logger.stagePop();
} // initialize
// ----------------------------------------------------------------------
@@ -511,6 +514,13 @@
const scalar_array&,
const scalar_array&,
const bool);
+ /// Member prototype for _constrainSolnSpaceImproveXD()
+ typedef void (pylith::faults::FaultCohesiveDyn::*constrainSolnSpaceImprove_fn_type)
+ (scalar_array*,
+ scalar_array*,
+ const scalar_array&,
+ const scalar_array&,
+ const scalar_array&);
assert(0 != fields);
assert(0 != _quadrature);
@@ -562,18 +572,25 @@
assert(!dLagrangeTpdtSection.isNull());
constrainSolnSpace_fn_type constrainSolnSpaceFn;
+ constrainSolnSpaceImprove_fn_type constrainSolnSpaceImproveFn;
switch (spaceDim) { // switch
case 1:
constrainSolnSpaceFn =
&pylith::faults::FaultCohesiveDyn::_constrainSolnSpace1D;
+ constrainSolnSpaceImproveFn =
+ &pylith::faults::FaultCohesiveDyn::_constrainSolnSpaceImprove1D;
break;
case 2:
constrainSolnSpaceFn =
&pylith::faults::FaultCohesiveDyn::_constrainSolnSpace2D;
+ constrainSolnSpaceImproveFn =
+ &pylith::faults::FaultCohesiveDyn::_constrainSolnSpaceImprove2D;
break;
case 3:
constrainSolnSpaceFn =
&pylith::faults::FaultCohesiveDyn::_constrainSolnSpace3D;
+ constrainSolnSpaceImproveFn =
+ &pylith::faults::FaultCohesiveDyn::_constrainSolnSpaceImprove3D;
break;
default :
assert(0);
@@ -700,25 +717,25 @@
// Use fault constitutive model to compute traction associated with
// friction.
- dLagrangeTpdtVertex = 0.0;
+ dTractionTpdtVertex = 0.0;
const bool iterating = true; // Iterating to get friction
CALL_MEMBER_FN(*this,
- constrainSolnSpaceFn)(&dLagrangeTpdtVertex,
+ constrainSolnSpaceFn)(&dTractionTpdtVertex,
t, slipTpdtVertex, slipRateVertex,
tractionTpdtVertex,
iterating);
// Rotate increment in traction back to global coordinate system.
- dLagrangeTpdtVertexGlobal = 0.0;
+ dLagrangeTpdtVertex = 0.0;
for (int iDim=0; iDim < spaceDim; ++iDim) {
for (int jDim=0; jDim < spaceDim; ++jDim) {
- dLagrangeTpdtVertexGlobal[iDim] +=
- orientationVertex[jDim*spaceDim+iDim] * dLagrangeTpdtVertex[jDim];
+ dLagrangeTpdtVertex[iDim] +=
+ orientationVertex[jDim*spaceDim+iDim] * dTractionTpdtVertex[jDim];
} // for
// Add in potential contribution from adjusting Lagrange
// multiplier for fault normal DOF of trial solution in Step 1.
- dLagrangeTpdtVertexGlobal[iDim] +=
+ dLagrangeTpdtVertex[iDim] +=
orientationVertex[indexN*spaceDim+iDim] * dTractionTpdtVertexNormal;
} // for
@@ -730,7 +747,7 @@
std::cout << ", slipRateVertex: ";
for (int iDim=0; iDim < spaceDim; ++iDim)
std::cout << " " << slipRateVertex[iDim];
- std::cout << ", tractionVertex: ";
+ std::cout << ", tractionTpdtVertex: ";
for (int iDim=0; iDim < spaceDim; ++iDim)
std::cout << " " << tractionTpdtVertex[iDim];
std::cout << ", lagrangeTVertex: ";
@@ -739,19 +756,19 @@
std::cout << ", lagrangeTIncrVertex: ";
for (int iDim=0; iDim < spaceDim; ++iDim)
std::cout << " " << lagrangeTIncrVertex[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 << ", dLagrangeTpdtVertexGlobal: ";
- for (int iDim=0; iDim < spaceDim; ++iDim)
- std::cout << " " << dLagrangeTpdtVertexGlobal[iDim];
std::cout << std::endl;
#endif
// Set change in Lagrange multiplier
- assert(dLagrangeTpdtVertexGlobal.size() ==
+ assert(dLagrangeTpdtVertex.size() ==
dLagrangeTpdtSection->getFiberDimension(v_fault));
- dLagrangeTpdtSection->updatePoint(v_fault, &dLagrangeTpdtVertexGlobal[0]);
+ dLagrangeTpdtSection->updatePoint(v_fault, &dLagrangeTpdtVertex[0]);
// Update displacement in trial solution (if necessary) so that it
// conforms to physical constraints.
@@ -818,7 +835,7 @@
// If no change in the Lagrange multiplier computed from friction
// criterion, there are no updates, so continue.
- double dLagrangeTpdtMag = 0.0;
+ PylithScalar dLagrangeTpdtMag = 0.0;
for (int i=0; i < spaceDim; ++i)
dLagrangeTpdtMag += dLagrangeTpdtVertex[i]*dLagrangeTpdtVertex[i];
if (0.0 == dLagrangeTpdtMag)
@@ -937,28 +954,40 @@
dSlipTpdtVertex[indexN] = -slipTpdtVertex[indexN];
} // if
- // Prevent over-correction in slip resulting in backslip.
- double slipDot = 0.0;
- double tractionSlipDot = 0.0;
- // :TODO:
- //
- // Need slipTVertex, slipTpdtVertex, dSlipTpdtVertex
- // slipDot = (slipTpdtVertex - slipTVertex) dot dSlipVertex
- //
- // if slipDot < 0, dSlipVertex = -0.5*(slipTpdtVertex - slipTVertex)
-
+ // Improve estimate of slip and change in traction using dFriction/dD.
+ // Get friction properties and state variables.
+ _friction->retrievePropsStateVars(v_fault);
+
+ CALL_MEMBER_FN(*this,
+ constrainSolnSpaceImproveFn)(&dTractionTpdtVertex, &dSlipTpdtVertex,
+ slipTVertex, slipTpdtVertex,
+ tractionTpdtVertex);
+
+#if 0 // OBSOLETE? Move to ImproveFn?
+ // Prevent over-correction in slip resulting in backslip
+ PylithScalar slipDot = 0.0;
+ PylithScalar tractionSlipDot = 0.0;
for (int iDim=0; iDim < spaceDim-1; ++iDim) { // :TODO: Optimize this
- slipDot += (slipTpdtVertex[iDim] - slipTVertex[iDim]) * (slipTpdtVertex[iDim] + dSlipTpdtVertex[iDim] - slipTVertex[iDim]);
+ // Compute dot product between slip and increment in slip (want positive)
+ slipDot +=
+ (slipTpdtVertex[iDim] - slipTVertex[iDim]) *
+ (slipTpdtVertex[iDim] + dSlipTpdtVertex[iDim] - slipTVertex[iDim]);
+ // Compute dot product of traction and slip
tractionSlipDot += (tractionTpdtVertex[iDim] + dTractionTpdtVertex[iDim])
* (slipTpdtVertex[iDim] + dSlipTpdtVertex[iDim]);
} // for
- if (slipDot < 0.0 && tractionSlipDot < 0.0) {
+ if (slipDot < 0.0 &&
+ sqrt(fabs(slipDot)) > _zeroTolerance &&
+ tractionSlipDot < 0.0) {
+ // Correct backslip
dTractionTpdtVertex *= 0.5; // Use bisection as guess for traction
for (int iDim=0; iDim < spaceDim-1; ++iDim) {
- dSlipTpdtVertex[iDim] *= -0.5*(slipTpdtVertex[iDim] - slipTVertex[iDim]);
+ // Use bisection for slip
+ dSlipTpdtVertex[iDim] = 0.5*(slipTVertex[iDim] - slipTpdtVertex[iDim]);
} // for
- } // if
+ } // if/else
+#endif
// Update current estimate of slip from t to t+dt.
slipTpdtVertex += dSlipTpdtVertex;
@@ -989,20 +1018,12 @@
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 << ", slipTpdtVertex: ";
for (int iDim=0; iDim < spaceDim; ++iDim)
- std::cout << " " << slipTpdtVertex[iDim];
+ std::cout << " " << slipTpdtVertex[iDim]-dSlipTpdtVertex[iDim];
std::cout << ", dSlipTpdtVertex: ";
for (int iDim=0; iDim < spaceDim; ++iDim)
std::cout << " " << dSlipTpdtVertex[iDim];
- //std::cout << ", dDispRelVertex: ";
- //for (int iDim=0; iDim < spaceDim; ++iDim)
- // std::cout << " " << dDispRelVertex[iDim];
- std::cout << ", slipDot: " << slipDot
- << ", tractionSlipDot: " << tractionSlipDot;
std::cout << std::endl;
#endif
@@ -1671,13 +1692,13 @@
const ALE::Obj<RealSection>& tractionsSection = tractions->section();
if (tractionsSection.isNull()) {
ALE::MemoryLogger& logger = ALE::MemoryLogger::singleton();
- //logger.stagePush("Fault");
+ logger.stagePush("FaultFields");
const topology::Field<topology::SubMesh>& dispRel =
_fields->get("relative disp");
tractions->cloneSection(dispRel);
- //logger.stagePop();
+ logger.stagePop();
} // if
const PylithScalar pressureScale = _normalizer->pressureScale();
tractions->label("traction");
@@ -2198,14 +2219,14 @@
// ----------------------------------------------------------------------
// Constrain solution space in 1-D.
void
-pylith::faults::FaultCohesiveDyn::_constrainSolnSpace1D(scalar_array* dLagrangeTpdt,
+pylith::faults::FaultCohesiveDyn::_constrainSolnSpace1D(scalar_array* dTractionTpdt,
const PylithScalar t,
const scalar_array& slip,
const scalar_array& sliprate,
const scalar_array& tractionTpdt,
const bool iterating)
{ // _constrainSolnSpace1D
- assert(0 != dLagrangeTpdt);
+ assert(0 != dTractionTpdt);
if (fabs(slip[0]) < _zeroTolerance) {
// if compression, then no changes to solution
@@ -2213,7 +2234,7 @@
// if tension, then traction is zero.
const PylithScalar dlp = -tractionTpdt[0];
- (*dLagrangeTpdt)[0] = dlp;
+ (*dTractionTpdt)[0] = dlp;
} // else
PetscLogFlops(2);
@@ -2222,14 +2243,14 @@
// ----------------------------------------------------------------------
// Constrain solution space in 2-D.
void
-pylith::faults::FaultCohesiveDyn::_constrainSolnSpace2D(scalar_array* dLagrangeTpdt,
+pylith::faults::FaultCohesiveDyn::_constrainSolnSpace2D(scalar_array* dTractionTpdt,
const PylithScalar t,
const scalar_array& slip,
const scalar_array& slipRate,
const scalar_array& tractionTpdt,
const bool iterating)
{ // _constrainSolnSpace2D
- assert(0 != dLagrangeTpdt);
+ assert(dTractionTpdt);
const PylithScalar slipMag = fabs(slip[0]);
const PylithScalar slipRateMag = fabs(slipRate[0]);
@@ -2250,11 +2271,9 @@
// versus friction
const PylithScalar dlp = -(tractionShearMag - frictionStress) *
tractionTpdt[0] / tractionShearMag;
- (*dLagrangeTpdt)[0] = dlp;
- (*dLagrangeTpdt)[1] = 0.0;
+ (*dTractionTpdt)[0] = dlp;
} else {
- (*dLagrangeTpdt)[0] = -(*dLagrangeTpdt)[0];
- (*dLagrangeTpdt)[1] = 0.0;
+ // No shear stress and no friction.
} // if/else
} else {
// friction exceeds value necessary to stick
@@ -2265,8 +2284,8 @@
} // if/else
} else {
// if in tension, then traction is zero.
- (*dLagrangeTpdt)[0] = -tractionTpdt[0];
- (*dLagrangeTpdt)[1] = -tractionTpdt[1];
+ (*dTractionTpdt)[0] = -tractionTpdt[0];
+ (*dTractionTpdt)[1] = -tractionTpdt[1];
} // else
PetscLogFlops(8);
@@ -2275,14 +2294,14 @@
// ----------------------------------------------------------------------
// Constrain solution space in 3-D.
void
-pylith::faults::FaultCohesiveDyn::_constrainSolnSpace3D(scalar_array* dLagrangeTpdt,
+pylith::faults::FaultCohesiveDyn::_constrainSolnSpace3D(scalar_array* dTractionTpdt,
const PylithScalar t,
const scalar_array& slip,
const scalar_array& slipRate,
const scalar_array& tractionTpdt,
const bool iterating)
{ // _constrainSolnSpace3D
- assert(0 != dLagrangeTpdt);
+ assert(dTractionTpdt);
const PylithScalar slipShearMag = sqrt(slip[0] * slip[0] +
slip[1] * slip[1]);
@@ -2311,13 +2330,10 @@
const PylithScalar dlq = -(tractionShearMag - frictionStress) *
tractionTpdt[1] / tractionShearMag;
- (*dLagrangeTpdt)[0] = dlp;
- (*dLagrangeTpdt)[1] = dlq;
- (*dLagrangeTpdt)[2] = 0.0;
+ (*dTractionTpdt)[0] = dlp;
+ (*dTractionTpdt)[1] = dlq;
} else {
- (*dLagrangeTpdt)[0] = -(*dLagrangeTpdt)[0];
- (*dLagrangeTpdt)[0] = -(*dLagrangeTpdt)[0];
- (*dLagrangeTpdt)[2] = 0.0;
+ // No shear stress and no friction.
} // if/else
} else {
@@ -2329,13 +2345,329 @@
} // if/else
} else {
// if in tension, then traction is zero.
- (*dLagrangeTpdt)[0] = -tractionTpdt[0];
- (*dLagrangeTpdt)[1] = -tractionTpdt[1];
- (*dLagrangeTpdt)[2] = -tractionTpdt[2];
+ (*dTractionTpdt)[0] = -tractionTpdt[0];
+ (*dTractionTpdt)[1] = -tractionTpdt[1];
+ (*dTractionTpdt)[2] = -tractionTpdt[2];
} // else
PetscLogFlops(22);
} // _constrainSolnSpace3D
+// ----------------------------------------------------------------------
+// Constrain solution space in 1-D.
+void
+pylith::faults::FaultCohesiveDyn::_constrainSolnSpaceImprove1D(
+ scalar_array* dTractionTpdt,
+ scalar_array* dSlipTpdt,
+ const scalar_array& slipT,
+ const scalar_array& slipTpdt,
+ const scalar_array& tractionTpdt)
+{ // _constrainSolnSpaceImprove3D
+ assert(dTractionTpdt);
+ assert(dSlipTpdt);
+
+ // Improving slip estimate only applies in shear. Do nothing.
+
+ PetscLogFlops(0); // :TODO: Update this.
+} // _constrainSolnSpaceImprove1D
+
+
+// ----------------------------------------------------------------------
+// Constrain solution space in 3-D.
+void
+pylith::faults::FaultCohesiveDyn::_constrainSolnSpaceImprove2D(
+ scalar_array* dTractionTpdt,
+ scalar_array* dSlipTpdt,
+ const scalar_array& slipT,
+ const scalar_array& slipTpdt,
+ const scalar_array& tractionTpdt)
+{ // _constrainSolnSpaceImprove2D
+ assert(dTractionTpdt);
+ assert(dSlipTpdt);
+
+#if 0 // DEBUGGING
+ std::cout << "BEFORE improvement"
+ << ", dTractionTpdt:";
+ for (int i=0; i < 2; ++i)
+ std::cout << " " << (*dTractionTpdt)[i];
+ std::cout << ", dSlipTpdt:";
+ for (int i=0; i < 2; ++i)
+ std::cout << " " << (*dSlipTpdt)[i];
+ std::cout << std::endl;
+#endif
+
+ // Compute magnitude of slip and slip rate (with current increment).
+ const PylithScalar slipShearMag = fabs(slipTpdt[0]);
+ const PylithScalar slipShearNewMag = fabs(slipTpdt[0]+(*dSlipTpdt)[0]);
+ const PylithScalar slipRateMag = fabs(slipTpdt[0]-slipT[0]) / _dt;
+ const PylithScalar dSlipShearNewMag = fabs((*dSlipTpdt)[0]);
+
+ const PylithScalar tractionNormal = tractionTpdt[1];
+
+ // Friction stress for old estimate of slip is tractionTpdt +
+ // dTractionTpdt.
+ const PylithScalar frictionStress =
+ fabs(tractionTpdt[0]+(*dTractionTpdt)[0]);
+ const PylithScalar tractionShearMag = fabs(tractionTpdt[0]);
+
+ if (fabs(slipTpdt[1] + (*dSlipTpdt)[1]) < _zeroTolerance &&
+ tractionNormal < -_zeroTolerance &&
+ dSlipShearNewMag > 0.0) {
+ // if in compression and no opening, and changing slip
+
+ // Calculate slope (Jacobian) of friction at slip before adding
+ // new increment.
+ const PylithScalar slopeF =
+ _friction->calcFrictionSlope(slipShearMag, slipRateMag, tractionNormal);
+
+#if 0 // linear space
+ const PylithScalar slopeT =
+ (frictionStress - tractionShearMag) / dSlipShearNewMag;
+
+ // Set adjustments to increments to original increments as default.
+ PylithScalar dSlipShearNew2Mag = dSlipShearNewMag;
+ PylithScalar dTractionShearNew2Mag = frictionStress - tractionShearMag;
+ if (slopeF > 0.0 && tractionShearMag > frictionStress) {
+ // Strengthening, so reduce increment in slip
+ dSlipShearNew2Mag =
+ (tractionShearMag - frictionStress) / (slopeF-slopeT));
+ dTractionShearNew2Mag += slopeF * dSlipShearNew2Mag;
+ } // if
+ // Ignore other cases, because slip will exceed estimate based on
+ // elasticity
+
+#else // log space
+
+ const PylithScalar slopeT = (frictionStress - tractionShearMag) /
+ (log(slipShearNewMag) - log(slipShearMag));
+
+ // Set adjustments to increments to original increments as default.
+ PylithScalar dSlipShearNew2Mag = dSlipShearNewMag;
+ PylithScalar dTractionShearNew2Mag = frictionStress - tractionShearMag;
+ if (slopeF > 0.0 && tractionShearMag > frictionStress) {
+ // Strengthening, so reduce increment in slip
+ const PylithScalar slipShearMagEff = std::max(slipShearMag, _zeroTolerance);
+ dSlipShearNew2Mag = slipShearMagEff *
+ (-1.0 + exp((tractionShearMag - frictionStress) / (slopeF-slopeT)));
+ dTractionShearNew2Mag +=
+ slopeF * log((slipShearMagEff + dSlipShearNew2Mag)/slipShearMagEff);
+
+ } // if
+ // Ignore other cases, because slip will exceed estimate based on
+ // elasticity
+#endif
+
+ // Project slip and traction into vector components
+ (*dSlipTpdt)[0] *= dSlipShearNew2Mag / dSlipShearNewMag;
+
+ const PylithScalar dTractionTpdtMag = fabs((*dTractionTpdt)[0]);
+ assert(dTractionTpdtMag > 0.0);
+ (*dTractionTpdt)[0] *= fabs(dTractionShearNew2Mag) / dTractionTpdtMag;
+
+ // Prevent over-correction in slip resulting in backslip.
+ // Expect slip direction to match tractions
+
+ // Compute dot product between slip and increment in slip (want positive)
+ const PylithScalar slipDot =
+ (slipTpdt[0] - slipT[0]) * (slipTpdt[0] + (*dSlipTpdt)[0] - slipT[0]);
+ // Compute dot product of traction and slip
+ const PylithScalar tractionSlipDot =
+ (tractionTpdt[0] + (*dTractionTpdt)[0]) * (slipTpdt[0] + (*dSlipTpdt)[0]);
+ if (slipDot < 0.0 &&
+ sqrt(fabs(slipDot)) > _zeroTolerance &&
+ tractionSlipDot < 0.0) {
+ // Correct backslip
+ (*dTractionTpdt)[0] *= 0.5; // Use bisection as guess for traction
+ // Use bisection for slip
+#if 0 // linear space
+ (*dSlipTpdt)[0] = 0.5*(slipT[0] - slipTpdt[0]);
+#else // log space
+ assert(slipT[0] * slipTpdt[0] > 0.0);
+ (*dSlipTpdt)[0] *= sqrt(slipT[0] * slipTpdt[0]) / dSlipShearNew2Mag;
+#endif
+ } // if/else
+
+#if 0 // DEBUGGING
+ std::cout << "AFTER improvement"
+ << ", dTractionTpdt:";
+ for (int i=0; i < 2; ++i)
+ std::cout << " " << (*dTractionTpdt)[i];
+ std::cout << ", dSlipTpdt:";
+ for (int i=0; i < 2; ++i)
+ std::cout << " " << (*dSlipTpdt)[i];
+ std::cout << ", frictionStress: " << frictionStress
+ << ", tractionShearMag: " << tractionShearMag
+ << ", slopeF: " << slopeF
+ << ", slopeT: " << slopeT
+ << std::endl;
+#endif
+
+ } // if
+
+ PetscLogFlops(0); // :TODO: Update this.
+
+
+} // _constrainSolnSpaceImprove2D
+
+
+// ----------------------------------------------------------------------
+// Constrain solution space in 3-D.
+void
+pylith::faults::FaultCohesiveDyn::_constrainSolnSpaceImprove3D(
+ scalar_array* dTractionTpdt,
+ scalar_array* dSlipTpdt,
+ const scalar_array& slipT,
+ const scalar_array& slipTpdt,
+ const scalar_array& tractionTpdt)
+{ // _constrainSolnSpaceImprove3D
+ assert(dTractionTpdt);
+ assert(dSlipTpdt);
+
+#if 0 // DEBUGGING
+ std::cout << "BEFORE improvement"
+ << ", dTractionTpdt:";
+ for (int i=0; i < 2; ++i)
+ std::cout << " " << (*dTractionTpdt)[i];
+ std::cout << ", dSlipTpdt:";
+ for (int i=0; i < 2; ++i)
+ std::cout << " " << (*dSlipTpdt)[i];
+ std::cout << std::endl;
+#endif
+
+ // Compute magnitude of slip and slip rate (with current increment).
+ const PylithScalar slipShearMag =
+ sqrt(slipTpdt[0]*slipTpdt[0] +
+ slipTpdt[1]*slipTpdt[1]);
+ const PylithScalar slipShearNewMag =
+ sqrt(pow(slipTpdt[0]+(*dSlipTpdt)[0], 2) +
+ pow(slipTpdt[1]+(*dSlipTpdt)[1], 2));
+ const PylithScalar slipRateMag =
+ sqrt(pow(slipTpdt[0]-slipT[0], 2) +
+ pow(slipTpdt[1]-slipT[1], 2)) / _dt;
+ const PylithScalar dSlipShearNewMag =
+ sqrt(pow((*dSlipTpdt)[0], 2) +
+ pow((*dSlipTpdt)[1], 2));
+
+ // Friction stress for old estimate of slip is tractionTpdt +
+ // dTractionTpdt.
+ const PylithScalar frictionStress =
+ sqrt(pow(tractionTpdt[0]+(*dTractionTpdt)[0], 2) +
+ pow(tractionTpdt[1]+(*dTractionTpdt)[1], 2));
+ const PylithScalar tractionShearMag =
+ sqrt(tractionTpdt[0]*tractionTpdt[0] +
+ tractionTpdt[1]*tractionTpdt[1]);
+
+ const PylithScalar tractionNormal = tractionTpdt[2] + (*dTractionTpdt)[2];
+
+ if (fabs(slipTpdt[2] + (*dSlipTpdt)[2]) < _zeroTolerance &&
+ tractionNormal < -_zeroTolerance &&
+ dSlipShearNewMag > 0.0) {
+ // if in compression and no opening, and changing slip
+
+ // Calculate slope (Jacobian) of friction at slip before adding
+ // new increment.
+ const PylithScalar slopeF =
+ _friction->calcFrictionSlope(slipShearMag, slipRateMag, tractionNormal);
+
+#if 0 // linear space
+ const PylithScalar slopeT =
+ (frictionStress - tractionShearMag) / dSlipShearNewMag;
+
+ // Set adjustments to increments to original increments as default.
+ PylithScalar dSlipShearNew2Mag = dSlipShearNewMag;
+ PylithScalar dTractionShearNew2Mag = frictionStress - tractionShearMag;
+ if (slopeF > 0.0 && tractionShearMag > frictionStress) {
+ // Strengthening, so reduce increment in slip
+ dSlipShearNew2Mag =
+ (tractionShearMag - frictionStress) / (slopeF-slopeT));
+ dTractionShearNew2Mag += slopeF * dSlipShearNew2Mag;
+ } // if
+ // Ignore other cases, because slip will exceed estimate based on
+ // elasticity
+
+#else // log space
+
+ const PylithScalar slopeT = (frictionStress - tractionShearMag) /
+ (log(slipShearNewMag) - log(slipShearMag));
+
+ // Set adjustments to increments to original increments as default.
+ PylithScalar dSlipShearNew2Mag = dSlipShearNewMag;
+ PylithScalar dTractionShearNew2Mag = frictionStress - tractionShearMag;
+ if (slopeF > 0.0 && tractionShearMag > frictionStress) {
+ // Strengthening, so reduce increment in slip
+ const PylithScalar slipShearMagEff = std::max(slipShearMag, _zeroTolerance);
+ dSlipShearNew2Mag = slipShearMagEff *
+ (-1.0 + exp((tractionShearMag - frictionStress) / (slopeF-slopeT)));
+ dTractionShearNew2Mag +=
+ slopeF * log((slipShearMagEff + dSlipShearNew2Mag)/slipShearMagEff);
+
+ } // if
+ // Ignore other cases, because slip will exceed estimate based on
+ // elasticity
+#endif
+
+ // Project slip and traction into vector components
+ // keeping same direction as original
+ (*dSlipTpdt)[0] *= dSlipShearNew2Mag / dSlipShearNewMag;
+ (*dSlipTpdt)[1] *= dSlipShearNew2Mag / dSlipShearNewMag;
+
+ const PylithScalar dTractionTpdtMag =
+ sqrt(pow((*dTractionTpdt)[0], 2) +
+ pow((*dTractionTpdt)[1], 2));
+ assert(dTractionTpdtMag > 0.0);
+ (*dTractionTpdt)[0] *= fabs(dTractionShearNew2Mag) / dTractionTpdtMag;
+ (*dTractionTpdt)[1] *= fabs(dTractionShearNew2Mag) / dTractionTpdtMag;
+
+ // Prevent over-correction in slip resulting in backslip.
+ // Expect slip direction to match tractions
+
+ // Compute dot product between slip and increment in slip (want positive)
+ const PylithScalar slipDot =
+ (slipTpdt[0] - slipT[0]) * (slipTpdt[0] + (*dSlipTpdt)[0] - slipT[0]) +
+ (slipTpdt[1] - slipT[1]) * (slipTpdt[1] + (*dSlipTpdt)[1] - slipT[1]);
+ // Compute dot product of traction and slip
+ const PylithScalar tractionSlipDot =
+ (tractionTpdt[0] + (*dTractionTpdt)[0]) * (slipTpdt[0] + (*dSlipTpdt)[0])+
+ (tractionTpdt[1] + (*dTractionTpdt)[1]) * (slipTpdt[1] + (*dSlipTpdt)[1]);
+ if (slipDot < 0.0 &&
+ sqrt(fabs(slipDot)) > _zeroTolerance &&
+ tractionSlipDot < 0.0) {
+ // Correct backslip
+ (*dTractionTpdt)[0] *= 0.5; // Use bisection as guess for traction
+ (*dTractionTpdt)[1] *= 0.5;
+ // Use bisection for slip
+#if 0 // linear space
+ (*dSlipTpdt)[0] = 0.5*(slipT[0] - slipTpdt[0]);
+ (*dSlipTpdt)[1] = 0.5*(slipT[1] - slipTpdt[1]);
+#else // log space
+ assert(slipT[0] * slipTpdt[0] > 0.0);
+ assert(slipT[1] * slipTpdt[1] > 0.0);
+
+ (*dSlipTpdt)[0] *= sqrt(slipT[0] * slipTpdt[0]) / dSlipShearNew2Mag;
+ (*dSlipTpdt)[1] *= sqrt(slipT[1] * slipTpdt[1]) / dSlipShearNew2Mag;
+#endif
+ } // if
+
+#if 0 // DEBUGGING
+ std::cout << "AFTER improvement"
+ << ", dTractionTpdt:";
+ for (int i=0; i < 2; ++i)
+ std::cout << " " << (*dTractionTpdt)[i];
+ std::cout << ", dSlipTpdt:";
+ for (int i=0; i < 2; ++i)
+ std::cout << " " << (*dSlipTpdt)[i];
+ std::cout << ", frictionStress: " << frictionStress
+ << ", tractionShearMag: " << tractionShearMag
+ << ", slopeF: " << slopeF
+ << ", slopeT: " << slopeT
+ << std::endl;
+#endif
+
+} // if
+
+ PetscLogFlops(0); // :TODO: Update this.
+} // _constrainSolnSpaceImprove3D
+
+
// End of file
Modified: short/3D/PyLith/trunk/libsrc/pylith/faults/FaultCohesiveDyn.hh
===================================================================
--- short/3D/PyLith/trunk/libsrc/pylith/faults/FaultCohesiveDyn.hh 2012-02-11 05:37:58 UTC (rev 19616)
+++ short/3D/PyLith/trunk/libsrc/pylith/faults/FaultCohesiveDyn.hh 2012-02-11 05:38:06 UTC (rev 19617)
@@ -201,7 +201,7 @@
*/
void _sensitivityUpdateSoln(const bool negativeSide);
- /** Constrain solution space with lumped Jacobian in 1-D.
+ /** Constrain solution space in 1-D.
*
* @param dLagrangeTpdt Adjustment to Lagrange multiplier.
* @param t Current time.
@@ -216,7 +216,7 @@
const scalar_array& tractionTpdt,
const bool iterating =true);
- /** Constrain solution space with lumped Jacobian in 2-D.
+ /** Constrain solution space in 2-D.
*
* @param dLagrangeTpdt Adjustment to Lagrange multiplier.
* @param t Current time.
@@ -231,7 +231,7 @@
const scalar_array& tractionTpdt,
const bool iterating =true);
- /** Constrain solution space with lumped Jacobian in 3-D.
+ /** Constrain solution space in 3-D.
*
* @param dLagrangeTpdt Adjustment to Lagrange multiplier.
* @param t Current time.
@@ -246,6 +246,48 @@
const scalar_array& tractionTpdt,
const bool iterating =true);
+ /** Improve slip estimate when constraining solution space in 1-D.
+ *
+ * @param dTractionTpdt Adjustment to fault traction vector.
+ * @param dSlipTpdt Adjustment to fault slip vector.
+ * @param slipT Fault slip vector at time t.
+ * @param slipTpdt Fault slip vector at time t+dt (without adjustment).
+ * @param tractionTpdt Fault traction vector (without adjustment).
+ */
+ void _constrainSolnSpaceImprove1D(double_array* dTractionTpdt,
+ double_array* dSlipTpdt,
+ const double_array& slipT,
+ const double_array& slipTpdt,
+ const double_array& tractionTpdt);
+
+ /** Improve slip estimate when constraining solution space in 2-D.
+ *
+ * @param dTractionTpdt Adjustment to fault traction vector.
+ * @param dSlipTpdt Adjustment to fault slip vector.
+ * @param slipT Fault slip vector at time t.
+ * @param slipTpdt Fault slip vector at time t+dt (without adjustment).
+ * @param tractionTpdt Fault traction vector (without adjustment).
+ */
+ void _constrainSolnSpaceImprove2D(double_array* dTractionTpdt,
+ double_array* dSlipTpdt,
+ const double_array& slipT,
+ const double_array& slipTpdt,
+ const double_array& tractionTpdt);
+
+ /** Improve slip estimate when constraining solution space in 3-D.
+ *
+ * @param dTractionTpdt Adjustment to fault traction vector.
+ * @param dSlipTpdt Adjustment to fault slip vector.
+ * @param slipT Fault slip vector at time t.
+ * @param slipTpdt Fault slip vector at time t+dt (without adjustment).
+ * @param tractionTpdt Fault traction vector (without adjustment).
+ */
+ void _constrainSolnSpaceImprove3D(double_array* dTractionTpdt,
+ double_array* dSlipTpdt,
+ const double_array& slipT,
+ const double_array& slipTpdt,
+ const double_array& tractionTpdt);
+
// PRIVATE MEMBERS ////////////////////////////////////////////////////
private :
Modified: short/3D/PyLith/trunk/libsrc/pylith/faults/FaultCohesiveLagrange.cc
===================================================================
--- short/3D/PyLith/trunk/libsrc/pylith/faults/FaultCohesiveLagrange.cc 2012-02-11 05:37:58 UTC (rev 19616)
+++ short/3D/PyLith/trunk/libsrc/pylith/faults/FaultCohesiveLagrange.cc 2012-02-11 05:38:06 UTC (rev 19617)
@@ -105,7 +105,7 @@
*_faultMesh);
ALE::MemoryLogger& logger = ALE::MemoryLogger::singleton();
- //logger.stagePush("Fault");
+ logger.stagePush("FaultFields");
// Allocate dispRel field
const ALE::Obj<SieveSubMesh>& faultSieveMesh = _faultMesh->sieveMesh();
@@ -120,6 +120,8 @@
dispRel.vectorFieldType(topology::FieldBase::VECTOR);
dispRel.scale(_normalizer->lengthScale());
+ logger.stagePop();
+
const ALE::Obj<SieveSubMesh::label_sequence>& cells =
faultSieveMesh->heightStratum(0);
assert(!cells.isNull());
@@ -130,26 +132,12 @@
_quadrature->computeGeometry(*_faultMesh, cells);
#endif
- _fields->add("distribution", "distribution", pylith::topology::FieldBase::CELLS_FIELD, 1);
- topology::Field<topology::SubMesh>& dist = _fields->get("distribution");
- dist.allocate();
- const ALE::Obj<RealSection>& distSection = dist.section();
- assert(!distSection.isNull());
- const PylithScalar rank = (PylithScalar) distSection->commRank();
-
- // Loop over cells in fault mesh, set distribution
- for (SieveSubMesh::label_sequence::iterator c_iter = cellsBegin; c_iter
- != cellsEnd; ++c_iter) {
- distSection->updatePoint(*c_iter, &rank);
- } // for
-
// Compute orientation at vertices in fault mesh.
_calcOrientation(upDir);
// Compute tributary area for each vertex in fault mesh.
_calcArea();
- //logger.stagePop();
} // initialize
// ----------------------------------------------------------------------
@@ -1209,6 +1197,9 @@
PylithScalar jacobianDet = 0;
scalar_array refCoordsVertex(cohesiveDim);
+ ALE::MemoryLogger& logger = ALE::MemoryLogger::singleton();
+ logger.stagePush("FaultFields");
+
// Allocate orientation field.
scalar_array orientationVertex(orientationSize);
_fields->add("orientation", "orientation");
@@ -1226,6 +1217,8 @@
orientation.allocate();
orientation.zero();
+ logger.stagePop();
+
// Get fault cells.
const ALE::Obj<SieveSubMesh::label_sequence>& cells =
faultSieveMesh->heightStratum(0);
@@ -1455,6 +1448,9 @@
vertices->begin();
const SieveSubMesh::label_sequence::iterator verticesEnd = vertices->end();
+ ALE::MemoryLogger& logger = ALE::MemoryLogger::singleton();
+ logger.stagePush("FaultFields");
+
// Allocate area field.
_fields->add("area", "area");
topology::Field<topology::SubMesh>& area = _fields->get("area");
@@ -1468,6 +1464,8 @@
assert(!areaSection.isNull());
UpdateAddVisitor areaVisitor(*areaSection, &areaCell[0]);
+ logger.stagePop();
+
scalar_array coordinatesCell(numBasis * spaceDim);
const ALE::Obj<RealSection>& coordinates = faultSieveMesh->getRealSection(
"coordinates");
@@ -1554,13 +1552,13 @@
const ALE::Obj<RealSection>& tractionsSection = tractions->section();
if (tractionsSection.isNull()) {
ALE::MemoryLogger& logger = ALE::MemoryLogger::singleton();
- //logger.stagePush("Fault");
+ logger.stagePush("FaultFields");
const topology::Field<topology::SubMesh>& dispRel =
_fields->get("relative disp");
tractions->cloneSection(dispRel);
- //logger.stagePop();
+ logger.stagePop();
} // if
assert(!tractionsSection.isNull());
tractions->zero();
@@ -1716,7 +1714,7 @@
return;
ALE::MemoryLogger& logger = ALE::MemoryLogger::singleton();
- logger.stagePush("Output");
+ logger.stagePush("OutputFields");
// Create vector field; use same shape/chart as relative
// displacement field.
@@ -1742,7 +1740,7 @@
return;
ALE::MemoryLogger& logger = ALE::MemoryLogger::singleton();
- logger.stagePush("Output");
+ logger.stagePush("OutputFields");
// Create vector field; use same shape/chart as area field.
assert(0 != _faultMesh);
@@ -1857,10 +1855,41 @@
pylith::faults::FaultCohesiveLagrange::cellField(const char* name,
const topology::SolutionFields* fields)
{ // cellField
- if (0 == strcasecmp("distribution", name)) {
- const topology::Field<topology::SubMesh>& dist = _fields->get("distribution");
- return dist;
- }
+ if (0 == strcasecmp("partition", name)) {
+
+ const ALE::Obj<SieveSubMesh>& faultSieveMesh = _faultMesh->sieveMesh();
+ assert(!faultSieveMesh.isNull());
+ 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();
+
+ ALE::MemoryLogger& logger = ALE::MemoryLogger::singleton();
+ logger.stagePush("OutputFields");
+
+ const int fiberDim = 1;
+ _fields->add("partition", "partition",
+ pylith::topology::FieldBase::CELLS_FIELD, fiberDim);
+ topology::Field<topology::SubMesh>& partition = _fields->get("partition");
+ partition.allocate();
+ const ALE::Obj<RealSection>& partitionSection = partition.section();
+ assert(!partitionSection.isNull());
+
+ const double rank = (double) partitionSection->commRank();
+ // Loop over cells in fault mesh, set partition
+ for (SieveSubMesh::label_sequence::iterator c_iter = cellsBegin;
+ c_iter != cellsEnd;
+ ++c_iter) {
+ partitionSection->updatePoint(*c_iter, &rank);
+ } // for
+
+ logger.stagePop();
+
+ return partition;
+
+ } // if
+
// Should not reach this point if requested field was found
std::ostringstream msg;
msg << "Request for unknown cell field '" << name << "' for fault '"
Modified: short/3D/PyLith/trunk/libsrc/pylith/feassemble/IntegratorElasticity.cc
===================================================================
--- short/3D/PyLith/trunk/libsrc/pylith/feassemble/IntegratorElasticity.cc 2012-02-11 05:37:58 UTC (rev 19616)
+++ short/3D/PyLith/trunk/libsrc/pylith/feassemble/IntegratorElasticity.cc 2012-02-11 05:38:06 UTC (rev 19617)
@@ -464,7 +464,7 @@
if (!_outputFields->hasField("buffer (tensor)")) {
ALE::MemoryLogger& logger = ALE::MemoryLogger::singleton();
- logger.stagePush("Output");
+ logger.stagePush("OutputFields");
_outputFields->add("buffer (tensor)", "buffer");
topology::Field<topology::Mesh>& buffer =
_outputFields->get("buffer (tensor)");
Modified: short/3D/PyLith/trunk/libsrc/pylith/friction/FrictionModel.cc
===================================================================
--- short/3D/PyLith/trunk/libsrc/pylith/friction/FrictionModel.cc 2012-02-11 05:37:58 UTC (rev 19616)
+++ short/3D/PyLith/trunk/libsrc/pylith/friction/FrictionModel.cc 2012-02-11 05:38:06 UTC (rev 19617)
@@ -325,6 +325,28 @@
} // calcFriction
// ----------------------------------------------------------------------
+// Compute change in friction with change in slip.
+PylithScalar
+pylith::friction::FrictionModel::calcFrictionSlope(const PylithScalar slip,
+ const PylithScalar slipRate,
+ const PylithScalar normalTraction)
+{ // calcFrictionSlope
+ assert(_fieldsPropsStateVars);
+
+ assert(_propsFiberDim+_varsFiberDim == _propsStateVarsVertex.size());
+ const PylithScalar* propertiesVertex = &_propsStateVarsVertex[0];
+ const PylithScalar* stateVarsVertex = (_varsFiberDim > 0) ?
+ &_propsStateVarsVertex[_propsFiberDim] : 0;
+
+ const PylithScalar slope =
+ _calcFrictionSlope(slip, slipRate, normalTraction,
+ propertiesVertex, _propsFiberDim,
+ stateVarsVertex, _varsFiberDim);
+
+ return slope;
+} // calcFrictionSlope
+
+// ----------------------------------------------------------------------
// Update state variables (for next time step).
void
pylith::friction::FrictionModel::updateStateVars(const PylithScalar t,
Modified: short/3D/PyLith/trunk/libsrc/pylith/friction/FrictionModel.hh
===================================================================
--- short/3D/PyLith/trunk/libsrc/pylith/friction/FrictionModel.hh 2012-02-11 05:37:58 UTC (rev 19616)
+++ short/3D/PyLith/trunk/libsrc/pylith/friction/FrictionModel.hh 2012-02-11 05:38:06 UTC (rev 19617)
@@ -168,6 +168,21 @@
const PylithScalar slipRate,
const PylithScalar normalTraction);
+ /** Compute change in friction for a change in slip (Jacobian).
+ *
+ * @pre Must call retrievePropsAndVars for cell before calling
+ * calcFriction().
+ *
+ * @param slip Current slip at location.
+ * @param slipRate Current slip rate at location.
+ * @param normalTraction Normal traction at location.
+ *
+ * @returns Change in friction for a chance in slip (dT/dD).
+ */
+ PylithScalar calcFrictionSlope(const PylithScalar slip,
+ const PylithScalar slipRate,
+ const PylithScalar normalTraction);
+
/** Compute friction at vertex.
*
* @pre Must call retrievePropsAndVars for cell before calling
@@ -254,6 +269,8 @@
* @param numProperties Number of properties.
* @param stateVars State variables at location.
* @param numStateVars Number of state variables.
+ *
+ * @returns Friction (magnitude of shear traction) at vertex.
*/
virtual
PylithScalar _calcFriction(const PylithScalar t,
@@ -265,6 +282,27 @@
const PylithScalar* stateVars,
const int numStateVars) = 0;
+ /** Compute change in friction for a change in slip (Jacobian).
+ *
+ * @param slip Current slip at location.
+ * @param slipRate Current slip rate at location.
+ * @param normalTraction Normal traction at location.
+ * @param properties Properties at location.
+ * @param numProperties Number of properties.
+ * @param stateVars State variables at location.
+ * @param numStateVars Number of state variables.
+ *
+ * @returns Change in friction for a chance in slip (dT/dD).
+ */
+ virtual
+ PylithScalar _calcFrictionSlope(const PylithScalar slip,
+ const PylithScalar slipRate,
+ const PylithScalar normalTraction,
+ const PylithScalar* properties,
+ const int numProperties,
+ const PylithScalar* stateVars,
+ const int numStateVars) = 0;
+
/** Update state variables (for next time step).
*
* @param t Time in simulation.
Modified: short/3D/PyLith/trunk/libsrc/pylith/friction/RateStateAgeing.cc
===================================================================
--- short/3D/PyLith/trunk/libsrc/pylith/friction/RateStateAgeing.cc 2012-02-11 05:37:58 UTC (rev 19616)
+++ short/3D/PyLith/trunk/libsrc/pylith/friction/RateStateAgeing.cc 2012-02-11 05:38:06 UTC (rev 19617)
@@ -333,54 +333,79 @@
mu_f = a * asinh(sinhArg);
friction = -mu_f * normalTraction + properties[p_cohesion];
- std::cout << "slip: " << slip
- << ", slipRate: " << slipRate
- << ", stateVar: " << theta
- << ", bLnTermL: " << bLnTerm
- << ", expTerm: " << expTerm
- << ", sinhArg: " << sinhArg
- << ", mu_f: " << mu_f
- << std::endl;
#else
- const double slipRateLinear = _minSlipRate;
- const double slipRateFactor = 1.0e-3;
+ const PylithScalar slipRateLinear = _minSlipRate;
- const double f0 = properties[p_coef];
- const double a = properties[p_a];
- const double b = properties[p_b];
- const double L = properties[p_L];
- const double slipRate0 = properties[p_slipRate0];
- const double theta = stateVars[s_state];
+ const PylithScalar f0 = properties[p_coef];
+ const PylithScalar a = properties[p_a];
+ const PylithScalar b = properties[p_b];
+ const PylithScalar L = properties[p_L];
+ const PylithScalar slipRate0 = properties[p_slipRate0];
+ const PylithScalar theta = stateVars[s_state];
- if (slipRate > slipRateLinear) {
+ if (slipRate >= slipRateLinear) {
mu_f = f0 + a*log(slipRate / slipRate0) + b*log(slipRate0*theta/L);
- } else if (slipRate > slipRateFactor*slipRateLinear) {
- mu_f = f0 + a*log(slipRateLinear / slipRate0) + b*log(slipRate0*theta/L) -
- a*(1.0-slipRateFactor) *
- (1.0 - slipRate/slipRateLinear) / (1.0 - slipRateFactor);
} else {
mu_f = f0 + a*log(slipRateLinear / slipRate0) + b*log(slipRate0*theta/L) -
- a*(1.0-slipRateFactor);
+ a*(1.0 - slipRate/slipRateLinear);
} // else
friction = -mu_f * normalTraction + properties[p_cohesion];
+#if 0
std::cout << "slip: " << slip
<< ", slipRate: " << slipRate
<< ", stateVar: " << theta
<< ", mu_f: " << mu_f
<< std::endl;
+#endif
#endif
} // if
- PetscLogFlops(11);
+ PetscLogFlops(12);
return friction;
} // _calcFriction
// ----------------------------------------------------------------------
+// Compute change in friction for a change in slip (Jacobian).
+PylithScalar
+pylith::friction::RateStateAgeing::_calcFrictionSlope(const PylithScalar slip,
+ const PylithScalar slipRate,
+ const PylithScalar normalTraction,
+ const PylithScalar* properties,
+ const int numProperties,
+ const PylithScalar* stateVars,
+ const int numStateVars)
+{ // _calcFrictionSlope
+ assert(properties);
+ assert(_RateStateAgeing::numProperties == numProperties);
+ assert(numStateVars);
+ assert(_RateStateAgeing::numStateVars == numStateVars);
+
+ PylithScalar slope = 0.0;
+ if (normalTraction <= 0.0) {
+ // if fault is in compression
+
+ const PylithScalar a = properties[p_a];
+ const PylithScalar slipRate0 = properties[p_slipRate0];
+
+ const PylithScalar slipRateLinear = _minSlipRate;
+ const PylithScalar slipRateEff = std::max(slipRate, slipRateLinear);
+
+ //slope = -slipRateEff / (slipRate0 * normalTraction * a);
+ slope = -normalTraction * a; // log space
+ } // if
+
+ PetscLogFlops(5);
+
+ return slope;
+} // _calcFrictionSlope
+
+
+// ----------------------------------------------------------------------
// Update state variables (for next time step).
void
pylith::friction::RateStateAgeing::_updateStateVars(const PylithScalar t,
Modified: short/3D/PyLith/trunk/libsrc/pylith/friction/RateStateAgeing.hh
===================================================================
--- short/3D/PyLith/trunk/libsrc/pylith/friction/RateStateAgeing.hh 2012-02-11 05:37:58 UTC (rev 19616)
+++ short/3D/PyLith/trunk/libsrc/pylith/friction/RateStateAgeing.hh 2012-02-11 05:38:06 UTC (rev 19617)
@@ -143,6 +143,26 @@
const PylithScalar* stateVars,
const int numStateVars);
+ /** Compute change in friction for a change in slip (Jacobian).
+ *
+ * @param slip Current slip at location.
+ * @param slipRate Current slip rate at location.
+ * @param normalTraction Normal traction at location.
+ * @param properties Properties at location.
+ * @param numProperties Number of properties.
+ * @param stateVars State variables at location.
+ * @param numStateVars Number of state variables.
+ *
+ * @returns Change in friction for a chance in slip (dT/dD).
+ */
+ PylithScalar _calcFrictionSlope(const PylithScalar slip,
+ const PylithScalar slipRate,
+ const PylithScalar normalTraction,
+ const PylithScalar* properties,
+ const int numProperties,
+ const PylithScalar* stateVars,
+ const int numStateVars);
+
/** Update state variables (for next time step).
*
* @param t Time in simulation.
Modified: short/3D/PyLith/trunk/libsrc/pylith/friction/SlipWeakening.cc
===================================================================
--- short/3D/PyLith/trunk/libsrc/pylith/friction/SlipWeakening.cc 2012-02-11 05:37:58 UTC (rev 19616)
+++ short/3D/PyLith/trunk/libsrc/pylith/friction/SlipWeakening.cc 2012-02-11 05:38:06 UTC (rev 19617)
@@ -283,8 +283,8 @@
PylithScalar mu_f = 0.0;
if (normalTraction <= 0.0) {
// if fault is in compression
- const double slipPrev = stateVars[s_slipPrev];
- const double slipCum = stateVars[s_slipCum] + fabs(slip - slipPrev);
+ const PylithScalar slipPrev = stateVars[s_slipPrev];
+ const PylithScalar slipCum = stateVars[s_slipCum] + fabs(slip - slipPrev);
if (slipCum < properties[p_d0]) {
// if/else linear slip-weakening form of mu_f
@@ -294,17 +294,54 @@
} else {
mu_f = properties[p_coefD];
} // if/else
- friction = - mu_f * normalTraction + properties[p_cohesion];
+ friction = -mu_f * normalTraction + properties[p_cohesion];
} else { // else
friction = properties[p_cohesion];
} // if/else
- PetscLogFlops(6);
+ PetscLogFlops(10);
return friction;
} // _calcFriction
// ----------------------------------------------------------------------
+// Compute change in friction for a change in slip (Jacobian).
+PylithScalar
+pylith::friction::SlipWeakening::_calcFrictionSlope(const PylithScalar slip,
+ const PylithScalar slipRate,
+ const PylithScalar normalTraction,
+ const PylithScalar* properties,
+ const int numProperties,
+ const PylithScalar* stateVars,
+ const int numStateVars)
+{ // _calcFrictionSlope
+ assert(properties);
+ assert(_SlipWeakening::numProperties == numProperties);
+ assert(stateVars);
+ assert(_SlipWeakening::numStateVars == numStateVars);
+
+ PylithScalar slope = 0.0;
+ if (normalTraction <= 0.0) {
+ // if fault is in compression
+ const PylithScalar slipPrev = stateVars[s_slipPrev];
+ const PylithScalar slipCum = stateVars[s_slipCum] + fabs(slip - slipPrev);
+
+ if (slipCum < properties[p_d0]) {
+ // if/else linear slip-weakening form of mu_f
+ slope = -normalTraction * (properties[p_coefS] - properties[p_coefD])
+ / properties[p_d0];
+ } else {
+ slope = 0.0;
+ } // if/else
+ } // if
+
+ PetscLogFlops(7);
+
+ return slope;
+} // _calcFrictionSlope
+
+
+// ----------------------------------------------------------------------
// Update state variables (for next time step).
void
pylith::friction::SlipWeakening::_updateStateVars(const PylithScalar t,
@@ -332,6 +369,8 @@
stateVars[s_slipPrev] = slip;
stateVars[s_slipCum] = 0.0;
} // else
+
+ PetscLogFlops(3);
} // _updateStateVars
Modified: short/3D/PyLith/trunk/libsrc/pylith/friction/SlipWeakening.hh
===================================================================
--- short/3D/PyLith/trunk/libsrc/pylith/friction/SlipWeakening.hh 2012-02-11 05:37:58 UTC (rev 19616)
+++ short/3D/PyLith/trunk/libsrc/pylith/friction/SlipWeakening.hh 2012-02-11 05:38:06 UTC (rev 19617)
@@ -131,6 +131,26 @@
const PylithScalar* stateVars,
const int numStateVars);
+ /** Compute change in friction for a change in slip (Jacobian).
+ *
+ * @param slip Current slip at location.
+ * @param slipRate Current slip rate at location.
+ * @param normalTraction Normal traction at location.
+ * @param properties Properties at location.
+ * @param numProperties Number of properties.
+ * @param stateVars State variables at location.
+ * @param numStateVars Number of state variables.
+ *
+ * @returns Change in friction for a chance in slip (dT/dD).
+ */
+ PylithScalar _calcFrictionSlope(const PylithScalar slip,
+ const PylithScalar slipRate,
+ const PylithScalar normalTraction,
+ const PylithScalar* properties,
+ const int numProperties,
+ const PylithScalar* stateVars,
+ const int numStateVars);
+
/** Update state variables (for next time step).
*
* @param t Time in simulation.
Modified: short/3D/PyLith/trunk/libsrc/pylith/friction/SlipWeakeningTime.cc
===================================================================
--- short/3D/PyLith/trunk/libsrc/pylith/friction/SlipWeakeningTime.cc 2012-02-11 05:37:58 UTC (rev 19616)
+++ short/3D/PyLith/trunk/libsrc/pylith/friction/SlipWeakeningTime.cc 2012-02-11 05:38:06 UTC (rev 19617)
@@ -312,6 +312,19 @@
} // _calcFriction
// ----------------------------------------------------------------------
+// Compute change in friction for a change in slip (Jacobian).
+PylithScalar
+pylith::friction::SlipWeakeningTime::_calcFrictionSlope(const PylithScalar slip,
+ const PylithScalar slipRate,
+ const PylithScalar normalTraction,
+ const PylithScalar* properties,
+ const int numProperties,
+ const PylithScalar* stateVars,
+ const int numStateVars)
+{ // _calcFrictionSlope
+} // _calcFrictionSlope
+
+// ----------------------------------------------------------------------
// Update state variables (for next time step).
void
pylith::friction::SlipWeakeningTime::_updateStateVars(const PylithScalar t,
Modified: short/3D/PyLith/trunk/libsrc/pylith/friction/SlipWeakeningTime.hh
===================================================================
--- short/3D/PyLith/trunk/libsrc/pylith/friction/SlipWeakeningTime.hh 2012-02-11 05:37:58 UTC (rev 19616)
+++ short/3D/PyLith/trunk/libsrc/pylith/friction/SlipWeakeningTime.hh 2012-02-11 05:38:06 UTC (rev 19617)
@@ -125,6 +125,26 @@
const PylithScalar* stateVars,
const int numStateVars);
+ /** Compute change in friction for a change in slip (Jacobian).
+ *
+ * @param slip Current slip at location.
+ * @param slipRate Current slip rate at location.
+ * @param normalTraction Normal traction at location.
+ * @param properties Properties at location.
+ * @param numProperties Number of properties.
+ * @param stateVars State variables at location.
+ * @param numStateVars Number of state variables.
+ *
+ * @returns Change in friction for a chance in slip (dT/dD).
+ */
+ PylithScalar _calcFrictionSlope(const PylithScalar slip,
+ const PylithScalar slipRate,
+ const PylithScalar normalTraction,
+ const PylithScalar* properties,
+ const int numProperties,
+ const PylithScalar* stateVars,
+ const int numStateVars);
+
/** Update state variables (for next time step).
*
* @param t Time in simulation.
Modified: short/3D/PyLith/trunk/libsrc/pylith/friction/StaticFriction.cc
===================================================================
--- short/3D/PyLith/trunk/libsrc/pylith/friction/StaticFriction.cc 2012-02-11 05:37:58 UTC (rev 19616)
+++ short/3D/PyLith/trunk/libsrc/pylith/friction/StaticFriction.cc 2012-02-11 05:38:06 UTC (rev 19617)
@@ -167,4 +167,21 @@
} // _calcFriction
+// ----------------------------------------------------------------------
+// Compute change in friction for a change in slip (Jacobian).
+PylithScalar
+pylith::friction::StaticFriction::_calcFrictionSlope(const PylithScalar slip,
+ const PylithScalar slipRate,
+ const PylithScalar normalTraction,
+ const PylithScalar* properties,
+ const int numProperties,
+ const PylithScalar* stateVars,
+ const int numStateVars)
+{ // _calcFrictionSlope
+ const PylithScalar slope = 0;
+
+ return slope;
+} // _calcFrictionSlope
+
+
// End of file
Modified: short/3D/PyLith/trunk/libsrc/pylith/friction/StaticFriction.hh
===================================================================
--- short/3D/PyLith/trunk/libsrc/pylith/friction/StaticFriction.hh 2012-02-11 05:37:58 UTC (rev 19616)
+++ short/3D/PyLith/trunk/libsrc/pylith/friction/StaticFriction.hh 2012-02-11 05:38:06 UTC (rev 19617)
@@ -96,6 +96,26 @@
const PylithScalar* stateVars,
const int numStateVars);
+ /** Compute change in friction for a change in slip (Jacobian).
+ *
+ * @param slip Current slip at location.
+ * @param slipRate Current slip rate at location.
+ * @param normalTraction Normal traction at location.
+ * @param properties Properties at location.
+ * @param numProperties Number of properties.
+ * @param stateVars State variables at location.
+ * @param numStateVars Number of state variables.
+ *
+ * @returns Change in friction for a chance in slip (dT/dD).
+ */
+ PylithScalar _calcFrictionSlope(const PylithScalar slip,
+ const PylithScalar slipRate,
+ const PylithScalar normalTraction,
+ const PylithScalar* properties,
+ const int numProperties,
+ const PylithScalar* stateVars,
+ const int numStateVars);
+
// PRIVATE MEMBERS ////////////////////////////////////////////////////
private :
Modified: short/3D/PyLith/trunk/libsrc/pylith/friction/TimeWeakening.cc
===================================================================
--- short/3D/PyLith/trunk/libsrc/pylith/friction/TimeWeakening.cc 2012-02-11 05:37:58 UTC (rev 19616)
+++ short/3D/PyLith/trunk/libsrc/pylith/friction/TimeWeakening.cc 2012-02-11 05:38:06 UTC (rev 19617)
@@ -310,4 +310,21 @@
} // _updateStateVars
+// ----------------------------------------------------------------------
+// Compute change in friction for a change in slip (Jacobian).
+PylithScalar
+pylith::friction::TimeWeakening::_calcFrictionSlope(const PylithScalar slip,
+ const PylithScalar slipRate,
+ const PylithScalar normalTraction,
+ const PylithScalar* properties,
+ const int numProperties,
+ const PylithScalar* stateVars,
+ const int numStateVars)
+{ // _calcFrictionSlope
+ const PylithScalar slope = 0;
+
+ return slope;
+} // _calcFrictionSlope
+
+
// End of file
Modified: short/3D/PyLith/trunk/libsrc/pylith/friction/TimeWeakening.hh
===================================================================
--- short/3D/PyLith/trunk/libsrc/pylith/friction/TimeWeakening.hh 2012-02-11 05:37:58 UTC (rev 19616)
+++ short/3D/PyLith/trunk/libsrc/pylith/friction/TimeWeakening.hh 2012-02-11 05:38:06 UTC (rev 19617)
@@ -125,6 +125,26 @@
const PylithScalar* stateVars,
const int numStateVars);
+ /** Compute change in friction for a change in slip (Jacobian).
+ *
+ * @param slip Current slip at location.
+ * @param slipRate Current slip rate at location.
+ * @param normalTraction Normal traction at location.
+ * @param properties Properties at location.
+ * @param numProperties Number of properties.
+ * @param stateVars State variables at location.
+ * @param numStateVars Number of state variables.
+ *
+ * @returns Change in friction for a chance in slip (dT/dD).
+ */
+ PylithScalar _calcFrictionSlope(const PylithScalar slip,
+ const PylithScalar slipRate,
+ const PylithScalar normalTraction,
+ const PylithScalar* properties,
+ const int numProperties,
+ const PylithScalar* stateVars,
+ const int numStateVars);
+
/** Update state variables (for next time step).
*
* @param slip Current slip at location.
Modified: short/3D/PyLith/trunk/libsrc/pylith/materials/ElasticMaterial.cc
===================================================================
--- short/3D/PyLith/trunk/libsrc/pylith/materials/ElasticMaterial.cc 2012-02-11 05:37:58 UTC (rev 19616)
+++ short/3D/PyLith/trunk/libsrc/pylith/materials/ElasticMaterial.cc 2012-02-11 05:38:06 UTC (rev 19617)
@@ -304,7 +304,7 @@
return;
ALE::MemoryLogger& logger = ALE::MemoryLogger::singleton();
- logger.stagePush("Materials");
+ logger.stagePush("MaterialsFields");
assert(0 != _initialFields);
_initialFields->add("initial stress", "initial_stress");
@@ -450,7 +450,7 @@
return;
ALE::MemoryLogger& logger = ALE::MemoryLogger::singleton();
- logger.stagePush("Materials");
+ logger.stagePush("MaterialsFields");
assert(0 != _initialFields);
_initialFields->add("initial strain", "initial_strain");
Modified: short/3D/PyLith/trunk/libsrc/pylith/materials/Material.cc
===================================================================
--- short/3D/PyLith/trunk/libsrc/pylith/materials/Material.cc 2012-02-11 05:37:58 UTC (rev 19616)
+++ short/3D/PyLith/trunk/libsrc/pylith/materials/Material.cc 2012-02-11 05:38:06 UTC (rev 19617)
@@ -117,7 +117,7 @@
assert(0 != quadrature);
ALE::MemoryLogger& logger = ALE::MemoryLogger::singleton();
- logger.stagePush("Materials");
+ logger.stagePush("MaterialsFields");
// Get quadrature information
const int numQuadPts = quadrature->numQuadPts();
@@ -384,7 +384,7 @@
} // if
if (!useCurrentField) {
ALE::MemoryLogger& logger = ALE::MemoryLogger::singleton();
- logger.stagePush("Output");
+ logger.stagePush("OutputFields");
field->newSection(cells, totalFiberDim);
field->allocate();
logger.stagePop();
@@ -454,7 +454,7 @@
} // if
if (!useCurrentField) {
ALE::MemoryLogger& logger = ALE::MemoryLogger::singleton();
- logger.stagePush("Output");
+ logger.stagePush("OutputFields");
field->newSection(cells, totalFiberDim);
field->allocate();
logger.stagePop();
Modified: short/3D/PyLith/trunk/libsrc/pylith/meshio/CellFilterAvg.cc
===================================================================
--- short/3D/PyLith/trunk/libsrc/pylith/meshio/CellFilterAvg.cc 2012-02-11 05:37:58 UTC (rev 19616)
+++ short/3D/PyLith/trunk/libsrc/pylith/meshio/CellFilterAvg.cc 2012-02-11 05:38:06 UTC (rev 19617)
@@ -121,7 +121,7 @@
// Allocate field if necessary
ALE::MemoryLogger& logger = ALE::MemoryLogger::singleton();
- logger.stagePush("Output");
+ logger.stagePush("OutputFields");
if (0 == _fieldAvg) {
_fieldAvg = new field_type(fieldIn.mesh());
assert(0 != _fieldAvg);
Modified: short/3D/PyLith/trunk/libsrc/pylith/meshio/DataWriterHDF5.cc
===================================================================
--- short/3D/PyLith/trunk/libsrc/pylith/meshio/DataWriterHDF5.cc 2012-02-11 05:37:58 UTC (rev 19616)
+++ short/3D/PyLith/trunk/libsrc/pylith/meshio/DataWriterHDF5.cc 2012-02-11 05:38:06 UTC (rev 19617)
@@ -146,8 +146,9 @@
CHECK_PETSC_ERROR(err);
// Account for censored cells
- int cellDepth = (sieveMesh->depth() == -1) ? -1 : 1;
- err = MPI_Allreduce(&cellDepth, &cellDepth, 1, MPI_INT, MPI_MAX,
+ int cellDepthLocal = (sieveMesh->depth() == -1) ? -1 : 1;
+ int cellDepth = 0;
+ err = MPI_Allreduce(&cellDepthLocal, &cellDepth, 1, MPI_INT, MPI_MAX,
sieveMesh->comm());CHECK_PETSC_ERROR(err);
const int depth = (0 == label) ? cellDepth : labelId;
const std::string labelName = (0 == label) ?
@@ -379,8 +380,9 @@
const ALE::Obj<typename mesh_type::SieveMesh>& sieveMesh =
field.mesh().sieveMesh();
assert(!sieveMesh.isNull());
- int cellDepth = (sieveMesh->depth() == -1) ? -1 : 1;
- err = MPI_Allreduce(&cellDepth, &cellDepth, 1, MPI_INT, MPI_MAX,
+ int cellDepthLocal = (sieveMesh->depth() == -1) ? -1 : 1;
+ int cellDepth = 0;
+ err = MPI_Allreduce(&cellDepthLocal, &cellDepth, 1, MPI_INT, MPI_MAX,
sieveMesh->comm());CHECK_PETSC_ERROR(err);
const int depth = (0 == label) ? cellDepth : labelId;
const std::string labelName = (0 == label) ?
Modified: short/3D/PyLith/trunk/libsrc/pylith/meshio/DataWriterHDF5Ext.cc
===================================================================
--- short/3D/PyLith/trunk/libsrc/pylith/meshio/DataWriterHDF5Ext.cc 2012-02-11 05:37:58 UTC (rev 19616)
+++ short/3D/PyLith/trunk/libsrc/pylith/meshio/DataWriterHDF5Ext.cc 2012-02-11 05:38:06 UTC (rev 19617)
@@ -162,8 +162,9 @@
// Write cells
// Account for censored cells
- int cellDepth = (sieveMesh->depth() == -1) ? -1 : 1;
- err = MPI_Allreduce(&cellDepth, &cellDepth, 1, MPI_INT, MPI_MAX,
+ int cellDepthLocal = (sieveMesh->depth() == -1) ? -1 : 1;
+ int cellDepth = 0;
+ err = MPI_Allreduce(&cellDepthLocal, &cellDepth, 1, MPI_INT, MPI_MAX,
sieveMesh->comm());CHECK_PETSC_ERROR(err);
const int depth = (0 == label) ? cellDepth : labelId;
const std::string labelName = (0 == label) ?
@@ -446,8 +447,9 @@
const ALE::Obj<typename mesh_type::SieveMesh>& sieveMesh =
field.mesh().sieveMesh();
assert(!sieveMesh.isNull());
- int cellDepth = (sieveMesh->depth() == -1) ? -1 : 1;
- err = MPI_Allreduce(&cellDepth, &cellDepth, 1, MPI_INT, MPI_MAX,
+ int cellDepthLocal = (sieveMesh->depth() == -1) ? -1 : 1;
+ int cellDepth = 0;
+ err = MPI_Allreduce(&cellDepthLocal, &cellDepth, 1, MPI_INT, MPI_MAX,
sieveMesh->comm());CHECK_PETSC_ERROR(err);
const int depth = (0 == label) ? cellDepth : labelId;
const std::string labelName = (0 == label) ?
Modified: short/3D/PyLith/trunk/libsrc/pylith/meshio/DataWriterVTK.cc
===================================================================
--- short/3D/PyLith/trunk/libsrc/pylith/meshio/DataWriterVTK.cc 2012-02-11 05:37:58 UTC (rev 19616)
+++ short/3D/PyLith/trunk/libsrc/pylith/meshio/DataWriterVTK.cc 2012-02-11 05:38:06 UTC (rev 19617)
@@ -255,8 +255,9 @@
// Cannot just use mesh->depth() because boundaries report the wrong thing
const ALE::Obj<SieveMesh>& sieveMesh = field.mesh().sieveMesh();
assert(!sieveMesh.isNull());
- int cellDepth = (sieveMesh->depth() == -1) ? -1 : 1;
- err = MPI_Allreduce(&cellDepth, &cellDepth, 1, MPI_INT, MPI_MAX,
+ int cellDepthLocal = (sieveMesh->depth() == -1) ? -1 : 1;
+ int cellDepth = 0;
+ err = MPI_Allreduce(&cellDepthLocal, &cellDepth, 1, MPI_INT, MPI_MAX,
sieveMesh->comm());CHECK_PETSC_ERROR(err);
const int depth = (!label) ? cellDepth : labelId;
const std::string labelName = (!label) ?
Modified: short/3D/PyLith/trunk/libsrc/pylith/meshio/MeshBuilder.cc
===================================================================
--- short/3D/PyLith/trunk/libsrc/pylith/meshio/MeshBuilder.cc 2012-02-11 05:37:58 UTC (rev 19616)
+++ short/3D/PyLith/trunk/libsrc/pylith/meshio/MeshBuilder.cc 2012-02-11 05:38:06 UTC (rev 19617)
@@ -222,7 +222,7 @@
// Memory logging
ALE::MemoryLogger& logger = ALE::MemoryLogger::singleton();
//logger.setDebug(fault->debug()/2);
- logger.stagePush("FaultCreation");
+ logger.stagePush("Creation");
if (0 == rank) {
assert(coordinates.size() == numVertices*spaceDim);
@@ -256,7 +256,7 @@
faultBd = ALE::Selection<SieveFlexMesh>::boundary(tmpMesh);
logger.stagePop();
- logger.stagePush("FaultStratification");
+ logger.stagePush("Stratification");
fault->stratify();
logger.stagePop();
} else {
@@ -265,7 +265,7 @@
faultBd = ALE::Selection<SieveFlexMesh>::boundary(tmpMesh);
logger.stagePop();
- logger.stagePush("FaultStratification");
+ logger.stagePush("Stratification");
fault->getSieve()->setChart(SieveMesh::sieve_type::chart_type());
fault->getSieve()->allocate();
fault->stratify();
@@ -275,20 +275,20 @@
#if defined(ALE_MEM_LOGGING)
std::cout
<< std::endl
- << "FaultCreation " << logger.getNumAllocations("FaultCreation")
- << " allocations " << logger.getAllocationTotal("FaultCreation")
+ << "FaultCreation " << logger.getNumAllocations("Creation")
+ << " allocations " << logger.getAllocationTotal("Creation")
<< " bytes" << std::endl
- << "FaultCreation " << logger.getNumDeallocations("FaultCreation")
- << " deallocations " << logger.getDeallocationTotal("FaultCreation")
+ << "FaultCreation " << logger.getNumDeallocations("Creation")
+ << " deallocations " << logger.getDeallocationTotal("Creation")
<< " bytes" << std::endl
- << "FaultStratification " << logger.getNumAllocations("FaultStratification")
- << " allocations " << logger.getAllocationTotal("FaultStratification")
+ << "FaultStratification " << logger.getNumAllocations("Stratification")
+ << " allocations " << logger.getAllocationTotal("Stratification")
<< " bytes" << std::endl
- << "FaultStratification " << logger.getNumDeallocations("FaultStratification")
- << " deallocations " << logger.getDeallocationTotal("FaultStratification")
+ << "FaultStratification " << logger.getNumDeallocations("Stratification")
+ << " deallocations " << logger.getDeallocationTotal("Stratification")
<< " bytes" << std::endl << std::endl;
#endif
Modified: short/3D/PyLith/trunk/libsrc/pylith/meshio/VertexFilterVecNorm.cc
===================================================================
--- short/3D/PyLith/trunk/libsrc/pylith/meshio/VertexFilterVecNorm.cc 2012-02-11 05:37:58 UTC (rev 19616)
+++ short/3D/PyLith/trunk/libsrc/pylith/meshio/VertexFilterVecNorm.cc 2012-02-11 05:38:06 UTC (rev 19617)
@@ -91,7 +91,7 @@
// Allocation field if necessary
ALE::MemoryLogger& logger = ALE::MemoryLogger::singleton();
- logger.stagePush("Output");
+ logger.stagePush("OutputFields");
if (0 == _fieldVecNorm) {
_fieldVecNorm = new field_type(fieldIn.mesh());
_fieldVecNorm->label("vector norm");
Modified: short/3D/PyLith/trunk/libsrc/pylith/topology/MeshOrder.cc
===================================================================
--- short/3D/PyLith/trunk/libsrc/pylith/topology/MeshOrder.cc 2012-02-11 05:37:58 UTC (rev 19616)
+++ short/3D/PyLith/trunk/libsrc/pylith/topology/MeshOrder.cc 2012-02-11 05:38:06 UTC (rev 19617)
@@ -74,12 +74,11 @@
} else {
const int numCells = cells->size();
const int numVertices = vertices->size();
- const int numEntities = numCells + numVertices;
_cellsNormal = ALE::Interval<point_type>(0, numCells);
_verticesNormal = ALE::Interval<point_type>(numCells, numCells+numVertices);
- _cellsCensored = ALE::Interval<point_type>(numEntities, numEntities);
- _verticesCensored = ALE::Interval<point_type>(numEntities, numEntities);
+ _verticesCensored = ALE::Interval<point_type>(numCells+numVertices, numCells+numVertices);
+ _cellsCensored = ALE::Interval<point_type>(numCells+numVertices, numCells+numVertices);
} // if/else
} // initialize
Modified: short/3D/PyLith/trunk/libsrc/pylith/topology/SubMesh.hh
===================================================================
--- short/3D/PyLith/trunk/libsrc/pylith/topology/SubMesh.hh 2012-02-11 05:37:58 UTC (rev 19616)
+++ short/3D/PyLith/trunk/libsrc/pylith/topology/SubMesh.hh 2012-02-11 05:38:06 UTC (rev 19617)
@@ -126,6 +126,24 @@
*/
int dimension(void) const;
+ /** Get representative cone size for mesh.
+ *
+ * @returns Representative cone size for mesh.
+ */
+ int coneSize(void) const;
+
+ /** Get number of vertices in mesh.
+ *
+ * @returns Number of vertices in mesh.
+ */
+ int numVertices(void) const;
+
+ /** Get number of cells in mesh.
+ *
+ * @returns Number of cells in mesh.
+ */
+ int numCells(void) const;
+
/** Get MPI communicator associated with mesh.
*
* @returns MPI communicator.
Modified: short/3D/PyLith/trunk/libsrc/pylith/topology/SubMesh.icc
===================================================================
--- short/3D/PyLith/trunk/libsrc/pylith/topology/SubMesh.icc 2012-02-11 05:37:58 UTC (rev 19616)
+++ short/3D/PyLith/trunk/libsrc/pylith/topology/SubMesh.icc 2012-02-11 05:38:06 UTC (rev 19617)
@@ -64,6 +64,29 @@
return (!_mesh.isNull()) ? _mesh->getDimension() : 0;
}
+// Get representative cone size for mesh.
+inline
+int
+pylith::topology::SubMesh::coneSize(void) const {
+
+ return (!_mesh.isNull() && numCells() > 0) ?
+ _mesh->getSieve()->getConeSize(*_mesh->heightStratum(1)->begin()) : 0;
+}
+
+// Get number of vertices in mesh.
+inline
+int
+pylith::topology::SubMesh::numVertices(void) const {
+ return (!_mesh.isNull() && _mesh->depth() > 0) ? _mesh->depthStratum(0)->size() : 0;
+}
+
+// Get number of cells in mesh.
+inline
+int
+pylith::topology::SubMesh::numCells(void) const {
+ return (!_mesh.isNull() && _mesh->height() > 0) ? _mesh->heightStratum(1)->size() : 0;
+}
+
// Get MPI communicator associated with mesh.
inline
const MPI_Comm
Modified: short/3D/PyLith/trunk/pylith/bc/AbsorbingDampers.py
===================================================================
--- short/3D/PyLith/trunk/pylith/bc/AbsorbingDampers.py 2012-02-11 05:37:58 UTC (rev 19616)
+++ short/3D/PyLith/trunk/pylith/bc/AbsorbingDampers.py 2012-02-11 05:38:06 UTC (rev 19617)
@@ -119,6 +119,12 @@
logEvent = "%sinit" % self._loggingPrefix
self._eventLogger.eventBegin(logEvent)
+ from pylith.mpi.Communicator import mpi_comm_world
+ comm = mpi_comm_world()
+
+ if 0 == comm.rank:
+ self._info.log("Initializing absorbing boundary '%s'." % self.label())
+
Integrator.initialize(self, totalTime, numTimeSteps, normalizer)
BoundaryCondition.initialize(self, totalTime, numTimeSteps, normalizer)
Modified: short/3D/PyLith/trunk/pylith/bc/DirichletBC.py
===================================================================
--- short/3D/PyLith/trunk/pylith/bc/DirichletBC.py 2012-02-11 05:37:58 UTC (rev 19616)
+++ short/3D/PyLith/trunk/pylith/bc/DirichletBC.py 2012-02-11 05:38:06 UTC (rev 19617)
@@ -92,6 +92,12 @@
logEvent = "%sinit" % self._loggingPrefix
self._eventLogger.eventBegin(logEvent)
+ from pylith.mpi.Communicator import mpi_comm_world
+ comm = mpi_comm_world()
+
+ if 0 == comm.rank:
+ self._info.log("Initializing Dirichlet boundary '%s'." % self.label())
+
self.normalizer(normalizer)
BoundaryCondition.initialize(self, totalTime, numTimeSteps, normalizer)
Modified: short/3D/PyLith/trunk/pylith/bc/DirichletBoundary.py
===================================================================
--- short/3D/PyLith/trunk/pylith/bc/DirichletBoundary.py 2012-02-11 05:37:58 UTC (rev 19616)
+++ short/3D/PyLith/trunk/pylith/bc/DirichletBoundary.py 2012-02-11 05:38:06 UTC (rev 19617)
@@ -112,8 +112,13 @@
logEvent = "%sinit" % self._loggingPrefix
self._eventLogger.eventBegin(logEvent)
+ from pylith.mpi.Communicator import mpi_comm_world
+ comm = mpi_comm_world()
+
+ if 0 == comm.rank:
+ self._info.log("Initializing Dirichlet boundary '%s'." % self.label())
+
DirichletBC.initialize(self, totalTime, numTimeSteps, normalizer)
-
self.output.initialize(normalizer)
self.output.writeInfo()
Modified: short/3D/PyLith/trunk/pylith/bc/Neumann.py
===================================================================
--- short/3D/PyLith/trunk/pylith/bc/Neumann.py 2012-02-11 05:37:58 UTC (rev 19616)
+++ short/3D/PyLith/trunk/pylith/bc/Neumann.py 2012-02-11 05:38:06 UTC (rev 19617)
@@ -121,6 +121,12 @@
logEvent = "%sinit" % self._loggingPrefix
self._eventLogger.eventBegin(logEvent)
+ from pylith.mpi.Communicator import mpi_comm_world
+ comm = mpi_comm_world()
+
+ if 0 == comm.rank:
+ self._info.log("Initializing Neumann boundary '%s'." % self.label())
+
Integrator.initialize(self, totalTime, numTimeSteps, normalizer)
BoundaryCondition.initialize(self, totalTime, numTimeSteps, normalizer)
Modified: short/3D/PyLith/trunk/pylith/bc/PointForce.py
===================================================================
--- short/3D/PyLith/trunk/pylith/bc/PointForce.py 2012-02-11 05:37:58 UTC (rev 19616)
+++ short/3D/PyLith/trunk/pylith/bc/PointForce.py 2012-02-11 05:38:06 UTC (rev 19617)
@@ -80,6 +80,12 @@
logEvent = "%sinit" % self._loggingPrefix
self._eventLogger.eventBegin(logEvent)
+ from pylith.mpi.Communicator import mpi_comm_world
+ comm = mpi_comm_world()
+
+ if 0 == comm.rank:
+ self._info.log("Initializing point forces '%s'." % self.label())
+
Integrator.initialize(self, totalTime, numTimeSteps, normalizer)
BoundaryCondition.initialize(self, totalTime, numTimeSteps, normalizer)
Modified: short/3D/PyLith/trunk/pylith/faults/FaultCohesiveDyn.py
===================================================================
--- short/3D/PyLith/trunk/pylith/faults/FaultCohesiveDyn.py 2012-02-11 05:37:58 UTC (rev 19616)
+++ short/3D/PyLith/trunk/pylith/faults/FaultCohesiveDyn.py 2012-02-11 05:38:06 UTC (rev 19617)
@@ -91,7 +91,7 @@
"slip_rate",
"traction"]},
'cell': \
- {'info': ["distribution"],
+ {'info': ["partition"],
'data': []}}
return
@@ -154,8 +154,7 @@
if 0 == comm.rank:
self._info.log("Initializing fault '%s'." % self.label())
- Integrator.initialize(self, totalTime, numTimeSteps, normalizer)
-
+ Integrator.initialize(self, totalTime, numTimeSteps, normalizer)
FaultCohesive.initialize(self, totalTime, numTimeSteps, normalizer)
self._eventLogger.eventEnd(logEvent)
Modified: short/3D/PyLith/trunk/pylith/faults/FaultCohesiveKin.py
===================================================================
--- short/3D/PyLith/trunk/pylith/faults/FaultCohesiveKin.py 2012-02-11 05:37:58 UTC (rev 19616)
+++ short/3D/PyLith/trunk/pylith/faults/FaultCohesiveKin.py 2012-02-11 05:38:06 UTC (rev 19617)
@@ -92,7 +92,7 @@
'data': ["slip",
"traction_change"]},
'cell': \
- {'info': ["distribution"],
+ {'info': ["partition"],
'data': []}}
return
Modified: short/3D/PyLith/trunk/pylith/meshio/OutputFaultDyn.py
===================================================================
--- short/3D/PyLith/trunk/pylith/meshio/OutputFaultDyn.py 2012-02-11 05:37:58 UTC (rev 19616)
+++ short/3D/PyLith/trunk/pylith/meshio/OutputFaultDyn.py 2012-02-11 05:38:06 UTC (rev 19617)
@@ -59,8 +59,7 @@
"traction"])
vertexDataFields.meta['tip'] = "Names of vertex data fields to output."
- cellInfoFields = pyre.inventory.list("cell_info_fields",
- default=["distribution"])
+ cellInfoFields = pyre.inventory.list("cell_info_fields", default=[])
cellInfoFields.meta['tip'] = "Names of cell info fields to output."
Modified: short/3D/PyLith/trunk/pylith/meshio/OutputFaultKin.py
===================================================================
--- short/3D/PyLith/trunk/pylith/meshio/OutputFaultKin.py 2012-02-11 05:37:58 UTC (rev 19616)
+++ short/3D/PyLith/trunk/pylith/meshio/OutputFaultKin.py 2012-02-11 05:38:06 UTC (rev 19617)
@@ -61,8 +61,7 @@
"traction_change"])
vertexDataFields.meta['tip'] = "Names of vertex data fields to output."
- cellInfoFields = pyre.inventory.list("cell_info_fields",
- default=["distribution"])
+ cellInfoFields = pyre.inventory.list("cell_info_fields", default=[])
cellInfoFields.meta['tip'] = "Names of cell info fields to output."
Modified: short/3D/PyLith/trunk/pylith/perf/Fault.py
===================================================================
--- short/3D/PyLith/trunk/pylith/perf/Fault.py 2012-02-11 05:37:58 UTC (rev 19616)
+++ short/3D/PyLith/trunk/pylith/perf/Fault.py 2012-02-11 05:38:06 UTC (rev 19617)
@@ -26,19 +26,45 @@
"""
Fault object for holding mesh memory and performance information.
"""
- def __init__(self):
+ def __init__(self, dimension=0, maxConeSize=0, numVertices=0, numCells=0):
"""
Constructor.
"""
+ self.dimension = dimension
+ self.coneSize = maxConeSize
+ self.nvertices = numVertices
+ self.ncells = numCells
+ self.cellType = None
+ self.initialize()
return
+ def initialize(self):
+ """
+ Initialize application.
+ """
+ from Mesh import cellTypes
+ try:
+ if self.coneSize > 0:
+ self.cellType = cellTypes[(self.dimension,self.coneSize)]
+ except:
+ raise ValueError("Unknown cell type '%s' for dim %d and cone size %d." %\
+ (self.cellType,self.dimension,self.coneSize))
+ return
+
+
def tabulate(self, memDict):
"""
Tabulate memory use.
"""
- memDict['Creation'] = 0
- memDict['Stratification'] = 0
- memDict['Coordinates'] = 0
+ dim = self.dimension
+ ncells = self.ncells
+ nvertices = self.nvertices
+ coneSize = self.coneSize
+
+ memDict['Creation'] = self.sizeInt * (2 * (coneSize*ncells + nvertices + ncells) + coneSize*ncells)
+ memDict['Stratification'] = 2 * self.sizeArrow * (nvertices + ncells)
+ memDict['Coordinates'] = (self.sizeDouble * dim * nvertices) + (2 * self.sizeInt * nvertices) + (2 * self.sizeInt * nvertices)
+
return
Modified: short/3D/PyLith/trunk/pylith/perf/MemoryLogger.py
===================================================================
--- short/3D/PyLith/trunk/pylith/perf/MemoryLogger.py 2012-02-11 05:37:58 UTC (rev 19616)
+++ short/3D/PyLith/trunk/pylith/perf/MemoryLogger.py 2012-02-11 05:38:06 UTC (rev 19617)
@@ -145,10 +145,12 @@
import pylith.perf.Field
if not stage in self.memory: self.memory[stage] = {}
+ if not 'Fields' in self.memory[stage]:
+ self.memory[stage]['Fields'] = {}
if not field is None:
fieldModel = pylith.perf.Field.Field(field.label(), field.sectionSize(),
field.chartSize())
- fieldModel.tabulate(self.memory[stage])
+ fieldModel.tabulate(self.memory[stage]['Fields'])
return
@@ -183,7 +185,8 @@
import pylith.perf.Fault
if not stage in self.memory: self.memory[stage] = {}
- faultModel = pylith.perf.Fault.Fault()
+ faultModel = pylith.perf.Fault.Fault(mesh.dimension(), mesh.coneSize(),
+ mesh.numVertices(), mesh.numCells())
faultModel.tabulate(self.memory[stage])
return
Modified: short/3D/PyLith/trunk/pylith/perf/Mesh.py
===================================================================
--- short/3D/PyLith/trunk/pylith/perf/Mesh.py 2012-02-11 05:37:58 UTC (rev 19616)
+++ short/3D/PyLith/trunk/pylith/perf/Mesh.py 2012-02-11 05:38:06 UTC (rev 19617)
@@ -22,59 +22,65 @@
from Memory import Memory
+cellTypes = {(1,2): 'line2',
+ (1,3): 'line3',
+ (2,3): 'tri3',
+ (2,6): 'tri6',
+ (2,4): 'quad4',
+ (2,9): 'quad9',
+ (3,4): 'tet4',
+ (3,8): 'hex8',
+ (3,10): 'tet10',
+ (3,27): 'hex27'
+ }
+
+
class Mesh(Memory):
- cellTypes = {(1,2): 'line2',
- (1,3): 'line3',
- (2,3): 'tri3',
- (2,6): 'tri6',
- (2,4): 'quad4',
- (2,9): 'quad9',
- (3,4): 'tet4',
- (3,8): 'hex8',
- (3,10): 'tet10',
- (3,27): 'hex27'
- }
-
"""
Mesh object for holding mesh memory and performance information.
"""
- def __init__(self, dimension = 0, maxConeSize = 0, numVertices = 0, numCells = 0):
+
+ def __init__(self, dimension=0, maxConeSize=0, numVertices=0, numCells=0):
"""
Constructor.
"""
self.dimension = dimension
- self.coneSize = maxConeSize
+ self.coneSize = maxConeSize
self.nvertices = numVertices
- self.ncells = numCells
- self.cellType = None
+ self.ncells = numCells
+ self.cellType = None
self.initialize()
return
+
def cellTypeInfo(cls, cellType):
for k,cT in cls.cellTypes.iteritems():
if cT == cellType:
return k
raise ValueError("Unknown cell type '%s'." % cellType)
+
def initialize(self):
"""
Initialize application.
"""
try:
if self.coneSize > 0:
- self.cellType = self.cellTypes[(self.dimension,self.coneSize)]
+ self.cellType = cellTypes[(self.dimension,self.coneSize)]
except:
- raise ValueError("Unknown cell type '%s' for dim %d and cone size %d." % (self.cellType,self.dimension,self.coneSize))
+ raise ValueError("Unknown cell type '%s' for dim %d and cone size %d." %\
+ (self.cellType,self.dimension,self.coneSize))
return
+
def tabulate(self, memDict):
"""
Tabulate memory use.
"""
- memDict['Creation'] = self.sizeInt * (2 * (self.coneSize*self.ncells + self.nvertices + self.ncells) + self.coneSize*self.ncells)
+ memDict['Creation'] = self.sizeInt * (2 * (self.coneSize*self.ncells + self.nvertices + self.ncells) + self.coneSize*self.ncells)
memDict['Stratification'] = 2 * self.sizeArrow * (self.nvertices + self.ncells)
# Here we have data + atlas (could use uniform) + bc (could use Section)
- memDict['Coordinates'] = (self.sizeDouble * self.dimension * self.nvertices) + (2 * self.sizeInt * self.nvertices) + (2 * self.sizeInt * self.nvertices)
+ memDict['Coordinates'] = (self.sizeDouble * self.dimension * self.nvertices) + (2 * self.sizeInt * self.nvertices) + (2 * self.sizeInt * self.nvertices)
memDict['Overlap'] = 0 # Don't know overlap
memDict['RealSections'] = 0 # Real sections should be elsewhere
return
Modified: short/3D/PyLith/trunk/pylith/topology/Mesh.py
===================================================================
--- short/3D/PyLith/trunk/pylith/topology/Mesh.py 2012-02-11 05:37:58 UTC (rev 19616)
+++ short/3D/PyLith/trunk/pylith/topology/Mesh.py 2012-02-11 05:38:06 UTC (rev 19617)
@@ -74,33 +74,6 @@
return
- def dimension(self):
- """
- Return the topological mesh dimension
- """
- return ModuleMesh.dimension(self)
-
-
- def coneSize(self):
- """
- Return the representative cone size, or number of vertices on a cell
- """
- return ModuleMesh.coneSize(self)
-
-
- def numVertices(self):
- """
- Return the number of vertices.
- """
- return ModuleMesh.numVertices(self)
-
-
- def numCells(self):
- """
- Return the number of cells.
- """
- return ModuleMesh.numCells(self)
-
def groupSizes(self):
"""
Return the name and number of vertices for each group
@@ -111,4 +84,5 @@
groups.append((name,ModuleMesh.groupSize(self, name)))
return groups
+
# End of file
Modified: short/3D/PyLith/trunk/pylith/topology/MeshGenerator.py
===================================================================
--- short/3D/PyLith/trunk/pylith/topology/MeshGenerator.py 2012-02-11 05:37:58 UTC (rev 19616)
+++ short/3D/PyLith/trunk/pylith/topology/MeshGenerator.py 2012-02-11 05:38:06 UTC (rev 19617)
@@ -113,13 +113,13 @@
for interface in interfaces:
if 0 == comm.rank:
self._info.log("Counting vertices for fault '%s'." % interface.label())
- nvertices = interface.numVertices(mesh)
+ nvertices = interface.numVerticesNoMesh(mesh)
firstLagrangeVertex += nvertices
firstFaultCell += nvertices
if interface.useLagrangeConstraints():
firstFaultCell += nvertices
for interface in interfaces:
- nvertices = interface.numVertices(mesh)
+ nvertices = interface.numVerticesNoMesh(mesh)
if 0 == comm.rank:
self._info.log("Adjusting topology for fault '%s' with %d vertices." % \
(interface.label(), nvertices))
Modified: short/3D/PyLith/trunk/pylith/utils/profiling.py
===================================================================
--- short/3D/PyLith/trunk/pylith/utils/profiling.py 2012-02-11 05:37:58 UTC (rev 19616)
+++ short/3D/PyLith/trunk/pylith/utils/profiling.py 2012-02-11 05:38:06 UTC (rev 19617)
@@ -16,7 +16,7 @@
# ----------------------------------------------------------------------
#
-## @file pylith/utils/profile.py
+## @file pylith/utils/profiling.py
# ----------------------------------------------------------------------
def resourceUsage():
@@ -42,7 +42,11 @@
"""
Get CPU time and memory usage as a string.
"""
- return "CPU time: %s, Memory usage: %.2f MB" % resourceUsage()
+ from pylith.mpi.Communicator import mpi_comm_world
+ comm = mpi_comm_world()
+ (cputime, memory) = resourceUsage()
+ return "[%d] CPU time: %s, Memory usage: %.2f MB" % \
+ (comm.rank, cputime, memory)
# End of file
Modified: short/3D/PyLith/trunk/tests/topology/test_meshmem.py
===================================================================
--- short/3D/PyLith/trunk/tests/topology/test_meshmem.py 2012-02-11 05:37:58 UTC (rev 19616)
+++ short/3D/PyLith/trunk/tests/topology/test_meshmem.py 2012-02-11 05:38:06 UTC (rev 19617)
@@ -60,8 +60,9 @@
Run test.
"""
- filenameIn = "data/tet4.exo"
+ #filenameIn = "data/tet4.exo"
#filenameIn = "tri3_200m_gradient.exo"
+ filenameIn = "tet4_150m.exo"
from pylith.perf.MemoryLogger import MemoryLogger
self.logger = MemoryLogger()
@@ -135,22 +136,22 @@
mesh = dmesh
- # Refine mesh (if necessary)
- from pylith.topology.RefineUniform import RefineUniform
- refiner = RefineUniform()
- rmesh = refiner.refine(mesh)
- rmesh.memLoggingStage = "RefinedMesh"
+ if False: # Refine mesh (if necessary)
+ from pylith.topology.RefineUniform import RefineUniform
+ refiner = RefineUniform()
+ rmesh = refiner.refine(mesh)
+ rmesh.memLoggingStage = "RefinedMesh"
- print "Unrefined mesh logging stage",mesh.memLoggingStage
- self.logger.logMesh(mesh.memLoggingStage, mesh)
- material.ncells = MeshOps_numMaterialCells(mesh, material.id())
- self.logger.logMaterial(mesh.memLoggingStage, material)
+ print "Unrefined mesh logging stage",mesh.memLoggingStage
+ self.logger.logMesh(mesh.memLoggingStage, mesh)
+ material.ncells = MeshOps_numMaterialCells(mesh, material.id())
+ self.logger.logMaterial(mesh.memLoggingStage, material)
- self.logger.logMesh(rmesh.memLoggingStage, rmesh)
- material.ncells = MeshOps_numMaterialCells(rmesh, material.id())
- self.logger.logMaterial(rmesh.memLoggingStage, material)
+ self.logger.logMesh(rmesh.memLoggingStage, rmesh)
+ material.ncells = MeshOps_numMaterialCells(rmesh, material.id())
+ self.logger.logMaterial(rmesh.memLoggingStage, material)
- self._showStatus("After refining mesh")
+ self._showStatus("After refining mesh")
return
Modified: short/3D/PyLith/trunk/unittests/libtests/topology/TestRefineUniform.cc
===================================================================
--- short/3D/PyLith/trunk/unittests/libtests/topology/TestRefineUniform.cc 2012-02-11 05:37:58 UTC (rev 19616)
+++ short/3D/PyLith/trunk/unittests/libtests/topology/TestRefineUniform.cc 2012-02-11 05:38:06 UTC (rev 19617)
@@ -148,7 +148,7 @@
faultA.id(100);
if (0 != data.faultA) {
faultA.label(data.faultA);
- const int nvertices = faultA.numVertices(*mesh);
+ const int nvertices = faultA.numVerticesNoMesh(*mesh);
firstLagrangeVertex += nvertices;
firstFaultCell += 2*nvertices; // shadow + Lagrange vertices
} // if
@@ -157,7 +157,7 @@
faultB.id(101);
if (0 != data.faultB) {
faultA.label(data.faultB);
- const int nvertices = faultB.numVertices(*mesh);
+ const int nvertices = faultB.numVerticesNoMesh(*mesh);
firstLagrangeVertex += nvertices;
firstFaultCell += 2*nvertices; // shadow + Lagrange vertices
} // if
More information about the CIG-COMMITS
mailing list