[cig-commits] r16317 - short/3D/PyLith/trunk/libsrc/materials

[willic3 at geodynamics.org]

Author: willic3
Date: 2010-02-22 20:02:28 -0800 (Mon, 22 Feb 2010)
New Revision: 16317

Modified:
   short/3D/PyLith/trunk/libsrc/materials/EffectiveStress.icc
   short/3D/PyLith/trunk/libsrc/materials/PowerLaw3D.cc
Log:
Changed the way dtStable is computed when stresses are zero.
Changed bracketing method to deal with cases where the effective stress
should be nearly zero.



Modified: short/3D/PyLith/trunk/libsrc/materials/EffectiveStress.icc
===================================================================
--- short/3D/PyLith/trunk/libsrc/materials/EffectiveStress.icc	2010-02-23 04:00:59 UTC (rev 16316)
+++ short/3D/PyLith/trunk/libsrc/materials/EffectiveStress.icc	2010-02-23 04:02:28 UTC (rev 16317)
@@ -68,6 +68,7 @@
   const double xMin = -1.0e-10;
   double x1 = *px1;
   double x2 = *px2;
+  double xLow = 1.0e-10;
 
   double funcValue1 = material->effStressFunc(x1);
   double funcValue2 = material->effStressFunc(x2);
@@ -82,13 +83,15 @@
 
     if (fabs(funcValue1) < fabs(funcValue2)) {
       x1 += bracketFactor * (x1 - x2);
-      x1 = std::max(x1, xMin);
+      x1 = std::max(x1, xLow);
       funcValue1 = material->effStressFunc(x1);
     } else {
       x2 += bracketFactor * (x1 - x2);
-      x2 = std::max(x2, xMin);
+      x2 = std::max(x2, xLow);
       funcValue2 = material->effStressFunc(x2);
     } // else
+    if (x1 == xLow || x2 == xLow)
+      xLow *= 0.5;
     ++iteration;
   } // while
 

Modified: short/3D/PyLith/trunk/libsrc/materials/PowerLaw3D.cc
===================================================================
--- short/3D/PyLith/trunk/libsrc/materials/PowerLaw3D.cc	2010-02-23 04:00:59 UTC (rev 16316)
+++ short/3D/PyLith/trunk/libsrc/materials/PowerLaw3D.cc	2010-02-23 04:02:28 UTC (rev 16317)
@@ -27,6 +27,7 @@
 #include <cassert> // USES assert()
 #include <cstring> // USES memcpy()
 #include <sstream> // USES std::ostringstream
+#include <iostream> // USES std::cout
 #include <stdexcept> // USES std::runtime_error
 
 // ----------------------------------------------------------------------
@@ -408,14 +409,12 @@
 			      stress[4],
 			      stress[5] };
   const double devStressProd = _scalarProduct(devStress, devStress);
-  const double effStress = sqrt(0.5 * devStressProd);
-  double dtTest = 1.0;
-  if (effStress != 0.0)
-    dtTest = 0.05 *
-      pow((referenceStress/effStress), (powerLawExp - 1.0)) *
-      (referenceStress/mu)/referenceStrainRate;
+  const double effStress = (devStressProd <= 0.0) ? referenceStress :
+    sqrt(0.5 * devStressProd);
+  const double dtStable = 0.05 *
+    pow((referenceStress/effStress), (powerLawExp - 1.0)) *
+    (referenceStress/mu)/referenceStrainRate;
 
-  const double dtStable = dtTest;
 #if 0 // DEBUGGING
   double maxwellTime = 10.0 * dtStable;
   std::cout << "Maxwell time:  " << maxwellTime << std::endl;
@@ -630,7 +629,7 @@
       
       const double effStressInitialGuess = effStressT;
 
-      double effStressTpdt =
+      effStressTpdt =
 	EffectiveStress::calculate<PowerLaw3D>(effStressInitialGuess,
 					       stressScale, this);
     } // if
@@ -1252,7 +1251,7 @@
 
     const double effStressInitialGuess = effStressT;
 
-    double effStressTpdt =
+    effStressTpdt =
       EffectiveStress::calculate<PowerLaw3D>(effStressInitialGuess,
 					     stressScale, this);