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

brad at geodynamics.org brad at geodynamics.org
Fri Jun 12 14:04:23 PDT 2009


Author: brad
Date: 2009-06-12 14:04:23 -0700 (Fri, 12 Jun 2009)
New Revision: 15220

Modified:
   short/3D/PyLith/trunk/libsrc/materials/EffectiveStress.icc
Log:
Adjusted convergence tolerance. Eliminated change in algorithm for slow convergence.

Modified: short/3D/PyLith/trunk/libsrc/materials/EffectiveStress.icc
===================================================================
--- short/3D/PyLith/trunk/libsrc/materials/EffectiveStress.icc	2009-06-12 20:37:10 UTC (rev 15219)
+++ short/3D/PyLith/trunk/libsrc/materials/EffectiveStress.icc	2009-06-12 21:04:23 UTC (rev 15220)
@@ -98,7 +98,6 @@
     throw std::runtime_error("Unable to bracket effective stress.");
 } // _bracket
 
-#include <iostream> // TEMPORARY
 // ----------------------------------------------------------------------
 // Find root using Newton's method with bisection.
 template<typename material_type>
@@ -111,7 +110,7 @@
   const int maxIterations = 100;
 
   // Desired accuracy for root. This is a bit arbitrary for now.
-  const double accuracy = 1.0e-18;
+  const double accuracy = 1.0e-10;
 
   // Organize search so that effStressFunc(xLow) is less than zero.
   double funcValueLow = material->effStressFunc(x1);
@@ -144,19 +143,12 @@
   while (iteration < maxIterations) {
     funcXHigh = (effStress - xHigh) * funcDeriv - funcValue;
     funcXLow = (effStress - xLow) * funcDeriv - funcValue;
-    std::cout << "low: " << funcXLow
-	      << ", high: " << funcXHigh
-	      << ", func: " << funcValue
-	      << std::endl;
     if (fabs(funcValue) < accuracy) {
       converged = true;
       break;
     } // if
-    // Use bisection if solution goes out of bounds or is not converging
-    // fast enough.
-    if ( (funcXHigh * funcXLow >= 0.0) ||
-	 (fabs(2.0 * funcValue) > fabs(dxPrevious * funcDeriv))) {
-      dxPrevious = dx;
+    // Use bisection if solution goes out of bounds.
+    if (funcXHigh * funcXLow >= 0.0) {
       dx = 0.5 * (xHigh - xLow);
       effStress = xLow + dx;
     } else {



More information about the CIG-COMMITS mailing list