[cig-commits] r19088 - in short/3D/PyLith/branches/v1.6-stable: libsrc/pylith/faults libsrc/pylith/friction modulesrc/faults modulesrc/friction pylith/faults pylith/friction tests/2d/frictionslide unittests/libtests/faults unittests/libtests/faults/data unittests/libtests/friction
brad at geodynamics.org
brad at geodynamics.org
Tue Oct 18 11:10:27 PDT 2011
Author: brad
Date: 2011-10-18 11:10:26 -0700 (Tue, 18 Oct 2011)
New Revision: 19088
Modified:
short/3D/PyLith/branches/v1.6-stable/libsrc/pylith/faults/FaultCohesiveDyn.cc
short/3D/PyLith/branches/v1.6-stable/libsrc/pylith/faults/FaultCohesiveDyn.hh
short/3D/PyLith/branches/v1.6-stable/libsrc/pylith/friction/RateStateAgeing.cc
short/3D/PyLith/branches/v1.6-stable/libsrc/pylith/friction/RateStateAgeing.hh
short/3D/PyLith/branches/v1.6-stable/modulesrc/faults/FaultCohesiveDyn.i
short/3D/PyLith/branches/v1.6-stable/modulesrc/friction/RateStateAgeing.i
short/3D/PyLith/branches/v1.6-stable/pylith/faults/FaultCohesive.py
short/3D/PyLith/branches/v1.6-stable/pylith/faults/FaultCohesiveDyn.py
short/3D/PyLith/branches/v1.6-stable/pylith/friction/FrictionModel.py
short/3D/PyLith/branches/v1.6-stable/pylith/friction/RateStateAgeing.py
short/3D/PyLith/branches/v1.6-stable/tests/2d/frictionslide/plot_friction.py
short/3D/PyLith/branches/v1.6-stable/tests/2d/frictionslide/pylithapp.cfg
short/3D/PyLith/branches/v1.6-stable/tests/2d/frictionslide/ratestate_weak.cfg
short/3D/PyLith/branches/v1.6-stable/unittests/libtests/faults/TestFaultCohesiveDyn.cc
short/3D/PyLith/branches/v1.6-stable/unittests/libtests/faults/TestFaultCohesiveDyn.hh
short/3D/PyLith/branches/v1.6-stable/unittests/libtests/faults/data/cohesivedyn.py
short/3D/PyLith/branches/v1.6-stable/unittests/libtests/friction/TestRateStateAgeing.cc
short/3D/PyLith/branches/v1.6-stable/unittests/libtests/friction/TestRateStateAgeing.hh
Log:
Added ability to change zero tolerances for constraining solution space and rate-state friction. Fixed small bug in constraining solution space (need to use zero tolerance to test for compressive normal stress due to roundoff errors).
Modified: short/3D/PyLith/branches/v1.6-stable/libsrc/pylith/faults/FaultCohesiveDyn.cc
===================================================================
--- short/3D/PyLith/branches/v1.6-stable/libsrc/pylith/faults/FaultCohesiveDyn.cc 2011-10-18 00:14:10 UTC (rev 19087)
+++ short/3D/PyLith/branches/v1.6-stable/libsrc/pylith/faults/FaultCohesiveDyn.cc 2011-10-18 18:10:26 UTC (rev 19088)
@@ -60,11 +60,9 @@
typedef ALE::ISieveVisitor::IndicesVisitor<RealSection,SieveSubMesh::order_type,PetscInt> IndicesVisitor;
// ----------------------------------------------------------------------
-const double pylith::faults::FaultCohesiveDyn::_zeroTolerance = 1.0e-10;
-
-// ----------------------------------------------------------------------
// Default constructor.
pylith::faults::FaultCohesiveDyn::FaultCohesiveDyn(void) :
+ _zeroTolerance(1.0e-10),
_dbInitialTract(0),
_friction(0),
_jacobian(0),
@@ -112,6 +110,21 @@
} // frictionModel
// ----------------------------------------------------------------------
+// Nondimensional tolerance for detecting near zero values.
+void
+pylith::faults::FaultCohesiveDyn::zeroTolerance(const double value)
+{ // zeroTolerance
+ if (value < 0.0) {
+ std::ostringstream msg;
+ msg << "Tolerance (" << value << ") for detecting values near zero for "
+ "fault " << label() << " must be nonnegative.";
+ throw std::runtime_error(msg.str());
+ } // if
+
+ _zeroTolerance = value;
+} // zeroTolerance
+
+// ----------------------------------------------------------------------
// Initialize fault. Determine orientation and setup boundary
void
pylith::faults::FaultCohesiveDyn::initialize(const topology::Mesh& mesh,
@@ -722,7 +735,7 @@
// Do not allow fault interpenetration and set fault opening to
// zero if fault is under compression.
- if (tractionNormal < 0.0 ||
+ if (tractionNormal < -_zeroTolerance ||
slipVertex[indexN] + dSlipVertex[indexN] < 0.0) {
dSlipVertex[indexN] = -slipVertex[indexN];
} // if
Modified: short/3D/PyLith/branches/v1.6-stable/libsrc/pylith/faults/FaultCohesiveDyn.hh
===================================================================
--- short/3D/PyLith/branches/v1.6-stable/libsrc/pylith/faults/FaultCohesiveDyn.hh 2011-10-18 00:14:10 UTC (rev 19087)
+++ short/3D/PyLith/branches/v1.6-stable/libsrc/pylith/faults/FaultCohesiveDyn.hh 2011-10-18 18:10:26 UTC (rev 19088)
@@ -67,6 +67,12 @@
*/
void frictionModel(friction::FrictionModel* const model);
+ /** Nondimensional tolerance for detecting near zero values.
+ *
+ * @param value Nondimensional tolerance
+ */
+ void zeroTolerance(const double value);
+
/** Initialize fault. Determine orientation and setup boundary
* condition parameters.
*
@@ -268,6 +274,9 @@
// PRIVATE MEMBERS ////////////////////////////////////////////////////
private :
+ /// Minimum resolvable value accounting for roundoff errors
+ double _zeroTolerance;
+
/// Database for initial tractions.
spatialdata::spatialdb::SpatialDB* _dbInitialTract;
@@ -279,9 +288,6 @@
PetscKSP _ksp; ///< PETSc KSP linear solver for sensitivity problem.
- /// Minimum resolvable value accounting for roundoff errors
- static const double _zeroTolerance;
-
// NOT IMPLEMENTED ////////////////////////////////////////////////////
private :
Modified: short/3D/PyLith/branches/v1.6-stable/libsrc/pylith/friction/RateStateAgeing.cc
===================================================================
--- short/3D/PyLith/branches/v1.6-stable/libsrc/pylith/friction/RateStateAgeing.cc 2011-10-18 00:14:10 UTC (rev 19087)
+++ short/3D/PyLith/branches/v1.6-stable/libsrc/pylith/friction/RateStateAgeing.cc 2011-10-18 18:10:26 UTC (rev 19088)
@@ -121,7 +121,8 @@
_RateStateAgeing::stateVars,
_RateStateAgeing::numStateVars,
_RateStateAgeing::dbStateVars,
- _RateStateAgeing::numDBStateVars))
+ _RateStateAgeing::numDBStateVars)),
+ _minSlipRate(1.0e-12)
{ // constructor
} // constructor
@@ -132,6 +133,21 @@
} // destructor
// ----------------------------------------------------------------------
+// Set floor for slip rate used in computing friction. Used to
+void
+pylith::friction::RateStateAgeing::minSlipRate(const double value)
+{ // minSlipRate
+ if (value < 0.0) {
+ std::ostringstream msg;
+ msg << "Minimum slip rate (" << value << ") for rate state friction model "
+ << label() << " must be nonnegative.";
+ throw std::runtime_error(msg.str());
+ } // if
+
+ _minSlipRate = value;
+} // minSlipRate
+
+// ----------------------------------------------------------------------
// Compute properties from values in spatial database.
void
pylith::friction::RateStateAgeing::_dbToProperties(
@@ -300,7 +316,7 @@
// Since regulatized friction -> 0 as slipRate -> 0, limit slip
// rate to some minimum value
- const double slipRateEff = std::max(1.0e-12, slipRate);
+ const double slipRateEff = std::max(1.0e-14, slipRate);
const double slipRate0 = properties[p_slipRate0];
const double a = properties[p_a];
@@ -358,7 +374,7 @@
// Since regulatized friction -> 0 as slipRate -> 0, limit slip
// rate to some minimum value
- const double slipRateEff = std::max(1.0e-12, slipRate);
+ const double slipRateEff = std::max(1.0e-14, slipRate);
const double dt = _dt;
const double thetaTVertex = stateVars[s_state];
Modified: short/3D/PyLith/branches/v1.6-stable/libsrc/pylith/friction/RateStateAgeing.hh
===================================================================
--- short/3D/PyLith/branches/v1.6-stable/libsrc/pylith/friction/RateStateAgeing.hh 2011-10-18 00:14:10 UTC (rev 19087)
+++ short/3D/PyLith/branches/v1.6-stable/libsrc/pylith/friction/RateStateAgeing.hh 2011-10-18 18:10:26 UTC (rev 19088)
@@ -58,6 +58,13 @@
/// Destructor.
~RateStateAgeing(void);
+ /** Set floor for slip rate used in computing friction. Used to
+ * avoid zero friction at zero slip rate.
+ *
+ * @param value Floor for slip rate.
+ */
+ void minSlipRate(const double value);
+
// PROTECTED METHODS //////////////////////////////////////////////////
protected :
@@ -155,6 +162,9 @@
// PRIVATE MEMBERS ////////////////////////////////////////////////////
private :
+ /// Floor for slip rate used in friction calculation.
+ double _minSlipRate;
+
/// Indices for properties in section and spatial database.
static const int p_coef;
static const int p_slipRate0;
Modified: short/3D/PyLith/branches/v1.6-stable/modulesrc/faults/FaultCohesiveDyn.i
===================================================================
--- short/3D/PyLith/branches/v1.6-stable/modulesrc/faults/FaultCohesiveDyn.i 2011-10-18 00:14:10 UTC (rev 19087)
+++ short/3D/PyLith/branches/v1.6-stable/modulesrc/faults/FaultCohesiveDyn.i 2011-10-18 18:10:26 UTC (rev 19088)
@@ -53,6 +53,12 @@
*/
void frictionModel(pylith::friction::FrictionModel* const model);
+ /** Nondimensional tolerance for detecting near zero values.
+ *
+ * @param value Nondimensional tolerance
+ */
+ void zeroTolerance(const double value);
+
/** Initialize fault. Determine orientation and setup boundary
* condition parameters.
*
Modified: short/3D/PyLith/branches/v1.6-stable/modulesrc/friction/RateStateAgeing.i
===================================================================
--- short/3D/PyLith/branches/v1.6-stable/modulesrc/friction/RateStateAgeing.i 2011-10-18 00:14:10 UTC (rev 19087)
+++ short/3D/PyLith/branches/v1.6-stable/modulesrc/friction/RateStateAgeing.i 2011-10-18 18:10:26 UTC (rev 19088)
@@ -36,6 +36,13 @@
/// Destructor.
~RateStateAgeing(void);
+ /** Set floor for slip rate used in computing friction. Used to
+ * avoid zero friction at zero slip rate.
+ *
+ * @param value Floor for slip rate.
+ */
+ void minSlipRate(const double value);
+
// PROTECTED METHODS //////////////////////////////////////////////
protected :
Modified: short/3D/PyLith/branches/v1.6-stable/pylith/faults/FaultCohesive.py
===================================================================
--- short/3D/PyLith/branches/v1.6-stable/pylith/faults/FaultCohesive.py 2011-10-18 00:14:10 UTC (rev 19087)
+++ short/3D/PyLith/branches/v1.6-stable/pylith/faults/FaultCohesive.py 2011-10-18 18:10:26 UTC (rev 19088)
@@ -94,6 +94,8 @@
self.inventory.meshFilename)
# Hardwire collocated quadrature
+ self.faultQuadrature.inventory.cell._configure()
+ self.faultQuadrature._configure()
self.faultQuadrature.cell.collocateQuad = True
self.faultQuadrature.cell.order = self.faultQuadrature.cell.degree
return
Modified: short/3D/PyLith/branches/v1.6-stable/pylith/faults/FaultCohesiveDyn.py
===================================================================
--- short/3D/PyLith/branches/v1.6-stable/pylith/faults/FaultCohesiveDyn.py 2011-10-18 00:14:10 UTC (rev 19087)
+++ short/3D/PyLith/branches/v1.6-stable/pylith/faults/FaultCohesiveDyn.py 2011-10-18 18:10:26 UTC (rev 19088)
@@ -42,7 +42,7 @@
Python object for managing FaultCohesiveDyn facilities and properties.
\b Properties
- @li None
+ @li \b zero_tolerance Tolerance for detecting zero values.
\b Facilities
@li \b db_initial_tractions Spatial database for initial tractions.
@@ -56,6 +56,10 @@
import pyre.inventory
+ zeroTolerance = pyre.inventory.float("zero_tolerance", default=1.0e-12,
+ validator=pyre.inventory.greaterEqual(0.0))
+ zeroTolerance.meta['tip'] = "Tolerance for detecting zero values."
+
db = pyre.inventory.facility("db_initial_tractions", family="spatial_database",
factory=NullComponent)
db.meta['tip'] = "Spatial database for initial tractions."
@@ -192,6 +196,7 @@
if not isinstance(self.inventory.db, NullComponent):
ModuleFaultCohesiveDyn.dbInitialTract(self, self.inventory.db)
ModuleFaultCohesiveDyn.frictionModel(self, self.inventory.friction)
+ ModuleFaultCohesiveDyn.zeroTolerance(self, self.inventory.zeroTolerance)
self.output = self.inventory.output
return
Modified: short/3D/PyLith/branches/v1.6-stable/pylith/friction/FrictionModel.py
===================================================================
--- short/3D/PyLith/branches/v1.6-stable/pylith/friction/FrictionModel.py 2011-10-18 00:14:10 UTC (rev 19087)
+++ short/3D/PyLith/branches/v1.6-stable/pylith/friction/FrictionModel.py 2011-10-18 18:10:26 UTC (rev 19088)
@@ -38,7 +38,7 @@
Validate descriptive label.
"""
if 0 == len(value):
- raise ValueError("Discriptive label for friction model not specified.")
+ raise ValueError("Descriptive label for friction model not specified.")
return value
Modified: short/3D/PyLith/branches/v1.6-stable/pylith/friction/RateStateAgeing.py
===================================================================
--- short/3D/PyLith/branches/v1.6-stable/pylith/friction/RateStateAgeing.py 2011-10-18 00:14:10 UTC (rev 19087)
+++ short/3D/PyLith/branches/v1.6-stable/pylith/friction/RateStateAgeing.py 2011-10-18 18:10:26 UTC (rev 19088)
@@ -33,6 +33,30 @@
Factory: friction_model.
"""
+ # INVENTORY //////////////////////////////////////////////////////////
+
+ class Inventory(FrictionModel.Inventory):
+ """
+ Python object for managing RateStateAgeing facilities and properties.
+ """
+
+ ## @class Inventory
+ ## Python object for managing RateStateAgeing facilities and properties.
+ ##
+ ## \b Properties
+ ## @li \b min_slip_rate Floor for nondimensional slip rate used in
+ ## friction calculation.
+ ##
+ ## \b Facilities
+ ## @li None
+
+ import pyre.inventory
+
+ minSlipRate = pyre.inventory.float("min_slip_rate", default=1.0e-12,
+ validator=pyre.inventory.greaterEqual(0.0))
+ minSlipRate.meta['tip'] = "Floor for nondimensional slip rate used in "\
+ "friction calculation."
+
# PUBLIC METHODS /////////////////////////////////////////////////////
def __init__(self, name="ratestateageing"):
@@ -58,6 +82,20 @@
# PRIVATE METHODS ////////////////////////////////////////////////////
+ def _configure(self):
+ """
+ Setup members using inventory.
+ """
+ try:
+ FrictionModel._configure(self)
+ ModuleRateStateAgeing.minSlipRate(self, self.inventory.minSlipRate)
+ except ValueError, err:
+ aliases = ", ".join(self.aliases)
+ raise ValueError("Error while configuring friction model "
+ "(%s):\n%s" % (aliases, err.message))
+ return
+
+
def _createModuleObj(self):
"""
Call constructor for module object for access to C++ object.
Modified: short/3D/PyLith/branches/v1.6-stable/tests/2d/frictionslide/plot_friction.py
===================================================================
--- short/3D/PyLith/branches/v1.6-stable/tests/2d/frictionslide/plot_friction.py 2011-10-18 00:14:10 UTC (rev 19087)
+++ short/3D/PyLith/branches/v1.6-stable/tests/2d/frictionslide/plot_friction.py 2011-10-18 18:10:26 UTC (rev 19088)
@@ -1,6 +1,6 @@
#!/usr/bin/env python
-sim = "ratestate_weak"
+sim = "ratestate_stable"
# ======================================================================
import tables
Modified: short/3D/PyLith/branches/v1.6-stable/tests/2d/frictionslide/pylithapp.cfg
===================================================================
--- short/3D/PyLith/branches/v1.6-stable/tests/2d/frictionslide/pylithapp.cfg 2011-10-18 00:14:10 UTC (rev 19087)
+++ short/3D/PyLith/branches/v1.6-stable/tests/2d/frictionslide/pylithapp.cfg 2011-10-18 18:10:26 UTC (rev 19088)
@@ -78,6 +78,7 @@
[pylithapp.timedependent.interfaces]
fault = pylith.faults.FaultCohesiveDyn
+fault.zero_tolerance = 1.0e-14
[pylithapp.timedependent.interfaces.fault]
id = 100
@@ -96,8 +97,8 @@
sub_pc_factor_shift_type = nonzero
# KSP
-ksp_rtol = 1.0e-8
-ksp_atol = 1.0e-12
+ksp_rtol = 1.0e-14
+ksp_atol = 1.0e-18
ksp_max_it = 500
ksp_gmres_restart = 20
@@ -106,8 +107,8 @@
ksp_converged_reason = true
# SNES
-snes_rtol = 1.0e-8
-snes_atol = 1.0e-10
+snes_rtol = 1.0e-12
+snes_atol = 1.0e-16
snes_max_it = 500
snes_monitor = true
Modified: short/3D/PyLith/branches/v1.6-stable/tests/2d/frictionslide/ratestate_weak.cfg
===================================================================
--- short/3D/PyLith/branches/v1.6-stable/tests/2d/frictionslide/ratestate_weak.cfg 2011-10-18 00:14:10 UTC (rev 19087)
+++ short/3D/PyLith/branches/v1.6-stable/tests/2d/frictionslide/ratestate_weak.cfg 2011-10-18 18:10:26 UTC (rev 19088)
@@ -13,6 +13,7 @@
[pylithapp.timedependent.interfaces.fault]
friction = pylith.friction.RateStateAgeing
friction.label = Rate state ageing
+friction.min_slip_rate = 1.0e-14
# Set the friction model parameters.
# reference coefficient of friction: 0.6
Modified: short/3D/PyLith/branches/v1.6-stable/unittests/libtests/faults/TestFaultCohesiveDyn.cc
===================================================================
--- short/3D/PyLith/branches/v1.6-stable/unittests/libtests/faults/TestFaultCohesiveDyn.cc 2011-10-18 00:14:10 UTC (rev 19087)
+++ short/3D/PyLith/branches/v1.6-stable/unittests/libtests/faults/TestFaultCohesiveDyn.cc 2011-10-18 18:10:26 UTC (rev 19088)
@@ -98,6 +98,20 @@
} // testDBInitialTract
// ----------------------------------------------------------------------
+// Test zeroTolerance().
+void
+pylith::faults::TestFaultCohesiveDyn::testZeroTolerance(void)
+{ // testZeroTolerance
+ FaultCohesiveDyn fault;
+
+ CPPUNIT_ASSERT_EQUAL(1.0e-10, fault._zeroTolerance); // default
+
+ const double value = 1.0e-20;
+ fault.zeroTolerance(value);
+ CPPUNIT_ASSERT_EQUAL(value, fault._zeroTolerance);
+ } // zeroTolerance
+
+// ----------------------------------------------------------------------
// Test initialize().
void
pylith::faults::TestFaultCohesiveDyn::testInitialize(void)
Modified: short/3D/PyLith/branches/v1.6-stable/unittests/libtests/faults/TestFaultCohesiveDyn.hh
===================================================================
--- short/3D/PyLith/branches/v1.6-stable/unittests/libtests/faults/TestFaultCohesiveDyn.hh 2011-10-18 00:14:10 UTC (rev 19087)
+++ short/3D/PyLith/branches/v1.6-stable/unittests/libtests/faults/TestFaultCohesiveDyn.hh 2011-10-18 18:10:26 UTC (rev 19088)
@@ -53,6 +53,7 @@
CPPUNIT_TEST( testConstructor );
CPPUNIT_TEST( testDBInitialTract );
+ CPPUNIT_TEST( testZeroTolerance );
// Tests in derived classes:
// testInitialize()
@@ -89,6 +90,9 @@
/// Test dbInitialTract().
void testDBInitialTract(void);
+ /// Test zeroTolerance().
+ void testZeroTolerance(void);
+
/// Test initialize().
void testInitialize(void);
Modified: short/3D/PyLith/branches/v1.6-stable/unittests/libtests/faults/data/cohesivedyn.py
===================================================================
--- short/3D/PyLith/branches/v1.6-stable/unittests/libtests/faults/data/cohesivedyn.py 2011-10-18 00:14:10 UTC (rev 19087)
+++ short/3D/PyLith/branches/v1.6-stable/unittests/libtests/faults/data/cohesivedyn.py 2011-10-18 18:10:26 UTC (rev 19088)
@@ -1,5 +1,5 @@
-cell = "hex8"
-dim = "3d"
+cell = "tri3d"
+dim = "2d"
testCase = "open"
import numpy
@@ -310,6 +310,7 @@
slipVertex[:,1] = 0
mask = slipVertex[:,1] < 0.0
slipVertex[mask,1] = 0
+ print "slip",slipVertex
slipVertex = faultToGlobal(slipVertex, C)
slipVertex = numpy.reshape(slipVertex, m)
Modified: short/3D/PyLith/branches/v1.6-stable/unittests/libtests/friction/TestRateStateAgeing.cc
===================================================================
--- short/3D/PyLith/branches/v1.6-stable/unittests/libtests/friction/TestRateStateAgeing.cc 2011-10-18 00:14:10 UTC (rev 19087)
+++ short/3D/PyLith/branches/v1.6-stable/unittests/libtests/friction/TestRateStateAgeing.cc 2011-10-18 18:10:26 UTC (rev 19088)
@@ -38,6 +38,20 @@
} // setUp
// ----------------------------------------------------------------------
+// Test floor for minimum slip rate.
+void
+pylith::friction::TestRateStateAgeing::testMinSlipRate(void)
+{ // testMinSlipRate
+ RateStateAgeing model;
+
+ CPPUNIT_ASSERT_EQUAL(1.0e-12, model._minSlipRate); // default
+
+ const double value = 1.0e-20;
+ model.minSlipRate(value);
+ CPPUNIT_ASSERT_EQUAL(value, model._minSlipRate);
+} // testMinSlipRate
+
+// ----------------------------------------------------------------------
// Test properties metadata.
void
pylith::friction::TestRateStateAgeing::testPropertiesMetadata(void)
Modified: short/3D/PyLith/branches/v1.6-stable/unittests/libtests/friction/TestRateStateAgeing.hh
===================================================================
--- short/3D/PyLith/branches/v1.6-stable/unittests/libtests/friction/TestRateStateAgeing.hh 2011-10-18 00:14:10 UTC (rev 19087)
+++ short/3D/PyLith/branches/v1.6-stable/unittests/libtests/friction/TestRateStateAgeing.hh 2011-10-18 18:10:26 UTC (rev 19088)
@@ -43,6 +43,7 @@
// CPPUNIT TEST SUITE /////////////////////////////////////////////////
CPPUNIT_TEST_SUITE( TestRateStateAgeing );
+ CPPUNIT_TEST( testMinSlipRate );
CPPUNIT_TEST( testPropertiesMetadata );
CPPUNIT_TEST( testStateVarsMetadata );
CPPUNIT_TEST( testDBToProperties );
@@ -63,6 +64,9 @@
/// Setup testing data.
void setUp(void);
+ /// Test floor for minimum slip rate.
+ void testMinSlipRate(void);
+
/// Test properties metadata.
void testPropertiesMetadata(void);
More information about the CIG-COMMITS
mailing list