[cig-commits] r16553 - in short/3D/PyLith/trunk: libsrc/faults libsrc/friction playpen/friction/bar_shearwave/quad4 pylith/faults
brad at geodynamics.org
brad at geodynamics.org
Wed Apr 14 09:55:47 PDT 2010
Author: brad
Date: 2010-04-14 09:55:46 -0700 (Wed, 14 Apr 2010)
New Revision: 16553
Modified:
short/3D/PyLith/trunk/libsrc/faults/FaultCohesiveDyn.cc
short/3D/PyLith/trunk/libsrc/friction/SlipWeakening.cc
short/3D/PyLith/trunk/playpen/friction/bar_shearwave/quad4/pylithapp.cfg
short/3D/PyLith/trunk/playpen/friction/bar_shearwave/quad4/shearwave-slipweakening.cfg
short/3D/PyLith/trunk/pylith/faults/FaultCohesiveDyn.py
Log:
Fixed bug where output of friction information was not hooked up to fault output. Removed some debugging output.
Modified: short/3D/PyLith/trunk/libsrc/faults/FaultCohesiveDyn.cc
===================================================================
--- short/3D/PyLith/trunk/libsrc/faults/FaultCohesiveDyn.cc 2010-04-14 08:01:02 UTC (rev 16552)
+++ short/3D/PyLith/trunk/libsrc/faults/FaultCohesiveDyn.cc 2010-04-14 16:55:46 UTC (rev 16553)
@@ -120,6 +120,7 @@
assert(0 != _friction);
assert(0 != _faultMesh);
assert(0 != _fields);
+ _friction->normalizer(*_normalizer);
_friction->initialize(*_faultMesh, _quadrature, _fields->get("area"));
const spatialdata::geocoords::CoordSys* cs = mesh.coordsys();
@@ -1030,6 +1031,14 @@
_calcTractions(&buffer, dispT);
return buffer;
+ } else if (_friction->hasProperty(name) || _friction->hasStateVar(name)) {
+ assert(0 != _fields);
+ if (!_fields->hasField("buffer (other)"))
+ _fields->add("buffer (other)", "buffer");
+ topology::Field<topology::SubMesh>& buffer = _fields->get("buffer (other)");
+ _friction->getField(&buffer, name);
+ return buffer;
+
} else {
std::ostringstream msg;
msg << "Request for unknown vertex field '" << name << "' for fault '"
@@ -1799,8 +1808,6 @@
const double dlp = dLagrangeTpdt[0];
(*slip)[0] -= Spp * dlp;
- std::cout << "Slip: (" << (*slip)[0] << ")" << std::endl;
-
PetscLogFlops(2);
} // _sensitivitySolveLumped1D
@@ -1833,8 +1840,6 @@
(*slip)[0] -= Spp * dlp;
(*slip)[1] -= Sqq * dlq;
- std::cout << "Slip: (" << (*slip)[0] << ", " << (*slip)[1] << ")" << std::endl;
-
PetscLogFlops(7);
} // _sensitivitySolveLumped2D
@@ -1873,8 +1878,6 @@
(*slip)[1] -= Sqq * dlq;
(*slip)[2] -= Srr * dlr;
- std::cout << "Slip: (" << (*slip)[0] << ", " << (*slip)[1] << ", " << (*slip)[2] << ")" << std::endl;
-
PetscLogFlops(9);
} // _sensitivitySolveLumped3D
@@ -1912,8 +1915,6 @@
{ // _constrainSolnSpace2D
assert(0 != dLagrangeTpdt);
- std::cout << "Normal traction:" << tractionTpdt[1] << std::endl;
-
const double slipMag = fabs(slip[0]);
const double slipRateMag = fabs(slipRate[0]);
@@ -1922,14 +1923,11 @@
if (tractionNormal < 0 && 0.0 == slip[1]) {
// if in compression and no opening
- std::cout << "FAULT IN COMPRESSION" << std::endl;
const double frictionStress = _friction->calcFriction(slipMag, slipRateMag,
tractionNormal);
- std::cout << "frictionStress: " << frictionStress << std::endl;
if (tractionShearMag > frictionStress ||
(tractionShearMag < frictionStress && slipMag > 0.0)) {
// traction is limited by friction, so have sliding
- std::cout << "LIMIT TRACTION, HAVE SLIDING" << std::endl;
// Update slip based on value required to stick versus friction
const double dlp = -(tractionShearMag - frictionStress) * area *
@@ -1938,13 +1936,10 @@
(*dLagrangeTpdt)[1] = 0.0;
} else {
// else friction exceeds value necessary, so stick
- std::cout << "STICK" << std::endl;
// no changes to solution
} // if/else
} else {
// if in tension, then traction is zero.
- std::cout << "FAULT IN TENSION" << std::endl;
-
(*dLagrangeTpdt)[0] = -tractionTpdt[0] * area;
(*dLagrangeTpdt)[1] = -tractionTpdt[1] * area;
} // else
@@ -1963,8 +1958,6 @@
{ // _constrainSolnSpace3D
assert(0 != dLagrangeTpdt);
- std::cout << "Normal traction:" << tractionTpdt[2] << std::endl;
-
const double slipShearMag = sqrt(slip[0] * slip[0] +
slip[1] * slip[1]);
double slipRateMag = sqrt(slipRate[0]*slipRate[0] +
@@ -1977,15 +1970,11 @@
if (tractionNormal < 0.0 && 0.0 == slip[2]) {
// if in compression and no opening
- std::cout << "FAULT IN COMPRESSION" << std::endl;
const double frictionStress =
_friction->calcFriction(slipShearMag, slipRateMag, tractionNormal);
- std::cout << "frictionStress: " << frictionStress << std::endl;
if (tractionShearMag > frictionStress ||
(tractionShearMag < frictionStress && slipShearMag > 0.0)) {
// traction is limited by friction, so have sliding
- std::cout << "LIMIT TRACTION, HAVE SLIDING" << std::endl;
-
// Update slip based on value required to stick versus friction
const double dlp = -(tractionShearMag - frictionStress) * area *
tractionTpdt[0] / tractionShearMag;
@@ -1998,13 +1987,10 @@
} else {
// else friction exceeds value necessary, so stick
- std::cout << "STICK" << std::endl;
// no changes to solution
} // if/else
} else {
// if in tension, then traction is zero.
- std::cout << "FAULT IN TENSION" << std::endl;
-
(*dLagrangeTpdt)[0] = -tractionTpdt[0] * area;
(*dLagrangeTpdt)[1] = -tractionTpdt[1] * area;
(*dLagrangeTpdt)[2] = -tractionTpdt[2] * area;
Modified: short/3D/PyLith/trunk/libsrc/friction/SlipWeakening.cc
===================================================================
--- short/3D/PyLith/trunk/libsrc/friction/SlipWeakening.cc 2010-04-14 08:01:02 UTC (rev 16552)
+++ short/3D/PyLith/trunk/libsrc/friction/SlipWeakening.cc 2010-04-14 16:55:46 UTC (rev 16553)
@@ -220,7 +220,7 @@
// Nondimensionalize state variables.
void
pylith::friction::SlipWeakening::_nondimStateVars(double* const values,
- const int nvalues) const
+ const int nvalues) const
{ // _nondimStateVars
assert(0 != _normalizer);
assert(0 != values);
@@ -236,7 +236,7 @@
// Dimensionalize state variables.
void
pylith::friction::SlipWeakening::_dimStateVars(double* const values,
- const int nvalues) const
+ const int nvalues) const
{ // _dimStateVars
assert(0 != _normalizer);
assert(0 != values);
@@ -302,7 +302,6 @@
stateVars[s_slipPrev] = stateVars[s_slipCum];
stateVars[s_slipCum] += fabs(slip - slipPrev);
- std::cout << "SLIP RATE: " << slipRate << std::endl;
} // _updateStateVars
Modified: short/3D/PyLith/trunk/playpen/friction/bar_shearwave/quad4/pylithapp.cfg
===================================================================
--- short/3D/PyLith/trunk/playpen/friction/bar_shearwave/quad4/pylithapp.cfg 2010-04-14 08:01:02 UTC (rev 16552)
+++ short/3D/PyLith/trunk/playpen/friction/bar_shearwave/quad4/pylithapp.cfg 2010-04-14 16:55:46 UTC (rev 16553)
@@ -42,6 +42,10 @@
# Set materials to an array with 1 material 'elastic'.
materials = [elastic]
+# Nondimensionalize problem using wave propagation parameters.
+normalizer = spatialdata.units.NondimElasticDynamic
+normalizer.shear_wave_speed = 1.0*km/s
+
[pylithapp.timedependent.formulation.time_step]
total_time = 0.0*s
dt = 0.05*s
Modified: short/3D/PyLith/trunk/playpen/friction/bar_shearwave/quad4/shearwave-slipweakening.cfg
===================================================================
--- short/3D/PyLith/trunk/playpen/friction/bar_shearwave/quad4/shearwave-slipweakening.cfg 2010-04-14 08:01:02 UTC (rev 16552)
+++ short/3D/PyLith/trunk/playpen/friction/bar_shearwave/quad4/shearwave-slipweakening.cfg 2010-04-14 16:55:46 UTC (rev 16553)
@@ -24,7 +24,7 @@
bc.x_neg = pylith.bc.AbsorbingDampers
[pylithapp.timedependent.formulation.time_step]
-total_time = 12.0*s
+total_time = 1.0*s
dt = 0.05*s
# ----------------------------------------------------------------------
@@ -115,6 +115,8 @@
# Give basename for VTK fault output.
[pylithapp.timedependent.interfaces.fault.output]
+vertex_info_fields = [strike_dir,normal_dir,initial_traction,static_coefficient,dynamic_coefficient,slip_weakening_parameter,cohesion]
+vertex_data_fields = [slip,traction,cumulative_slip,previous_slip]
writer.filename = output/shearwave-slipweakening-fault.vtk
# Give basename for VTK output of state variables.
Modified: short/3D/PyLith/trunk/pylith/faults/FaultCohesiveDyn.py
===================================================================
--- short/3D/PyLith/trunk/pylith/faults/FaultCohesiveDyn.py 2010-04-14 08:01:02 UTC (rev 16552)
+++ short/3D/PyLith/trunk/pylith/faults/FaultCohesiveDyn.py 2010-04-14 16:55:46 UTC (rev 16553)
@@ -103,6 +103,11 @@
if not isinstance(self.inventory.db, NullComponent):
self.availableFields['vertex']['info'] += ["initial_traction"]
+
+ self.availableFields['vertex']['info'] += \
+ self.friction.availableFields['vertex']['info']
+ self.availableFields['vertex']['data'] += \
+ self.friction.availableFields['vertex']['data']
return
More information about the CIG-COMMITS
mailing list