[cig-commits] r15065 - in short/3D/PyLith/trunk/libsrc: . materials

willic3 at geodynamics.org willic3 at geodynamics.org
Tue May 26 17:29:59 PDT 2009


Author: willic3
Date: 2009-05-26 17:29:59 -0700 (Tue, 26 May 2009)
New Revision: 15065

Modified:
   short/3D/PyLith/trunk/libsrc/Makefile.am
   short/3D/PyLith/trunk/libsrc/materials/EffectiveStress.cc
   short/3D/PyLith/trunk/libsrc/materials/EffectiveStress.hh
   short/3D/PyLith/trunk/libsrc/materials/Makefile.am
   short/3D/PyLith/trunk/libsrc/materials/PowerLaw3D.cc
   short/3D/PyLith/trunk/libsrc/materials/PowerLaw3D.hh
Log:
Finished debugging compile problems for PowerLaw3D and EffectiveStress.



Modified: short/3D/PyLith/trunk/libsrc/Makefile.am
===================================================================
--- short/3D/PyLith/trunk/libsrc/Makefile.am	2009-05-26 22:59:53 UTC (rev 15064)
+++ short/3D/PyLith/trunk/libsrc/Makefile.am	2009-05-27 00:29:59 UTC (rev 15065)
@@ -76,6 +76,8 @@
 	materials/ElasticIsotropic3D.cc \
 	materials/ViscoelasticMaxwell.cc \
 	materials/MaxwellIsotropic3D.cc \
+	materials/PowerLaw3D.cc \
+	materials/EffectiveStress.cc \
 	meshio/BinaryIO.cc \
 	meshio/GMVFile.cc \
 	meshio/GMVFileAscii.cc \

Modified: short/3D/PyLith/trunk/libsrc/materials/EffectiveStress.cc
===================================================================
--- short/3D/PyLith/trunk/libsrc/materials/EffectiveStress.cc	2009-05-26 22:59:53 UTC (rev 15064)
+++ short/3D/PyLith/trunk/libsrc/materials/EffectiveStress.cc	2009-05-27 00:29:59 UTC (rev 15065)
@@ -17,6 +17,7 @@
 #include "petsc.h" // USES PetscLogFlops
 
 #include <stdexcept> // USES std::runtime_error
+#include <cmath> // USES fabs
 
 // ----------------------------------------------------------------------
 // Compute effective stress, given an initial guess, a vector of parameters,
@@ -28,13 +29,12 @@
 pylith::materials::EffectiveStress::getEffStress(
 				 const double effStressInitialGuess,
 				 const double stressScale,
-				 EffStressStruct effStressParams,
+				 const EffStressStruct& effStressParams,
 				 effStressFunc_fn_type effStressFunc,
 				 effStressFuncDFunc_fn_type effStressFuncDFunc)
 { // getEffStress
   // Check parameters
-  if (effStressInitialGuess < 0.0)
-    throw std::runtime_error("Effective stress initial guess must be >= 0.");
+  assert(effStressInitialGuess >= 0.0);
 
   // Bracket the root.
   double x1 = 0.0;
@@ -48,13 +48,11 @@
   } // else
 
   PetscLogFlops(4);
-  _bracketEffStress(x1, x2, effStressParams, effStressFunc);
+  _bracketEffStress(&x1, &x2, effStressParams, effStressFunc);
 
   // Find effective stress using Newton's method with bisection.
-  const double xx1 = x1;
-  const double xx2 = x2;
-  const double effStress = _findEffStress(xx1,
-					  xx2,
+  const double effStress = _findEffStress(x1,
+					  x2,
 					  effStressParams,
 					  effStressFunc,
 					  effStressFuncDFunc);
@@ -69,7 +67,7 @@
 pylith::materials::EffectiveStress::_bracketEffStress(
 				     double* px1,
 				     double* px2,
-				     EffStressStruct& effStressParams,
+				     const EffStressStruct& effStressParams,
 				     effStressFunc_fn_type effStressFunc)
 { // _bracketEffStress
   // Arbitrary number of iterations to bracket the root
@@ -86,19 +84,19 @@
   int iteration = 0;
   bool bracketed = false;
 
-  while ((iteration < maxIterations) && (bracketed == false)) {
+  while (iteration < maxIterations) {
     if ((funcValue1 * funcValue2) < 0.0) {
       bracketed = true;
+      break;
+    }
+    if (fabs(funcValue1) < fabs(funcValue2)) {
+      x1 += bracketFactor * (x1 - x2);
+      x1 = std::max(x1, 0.0);
+      funcValue1 = effStressFunc(x1, effStressParams);
     } else {
-      if (std::abs(f1) < std::abs(f2)) {
-	x1 += bracketFactor * (x1 - x2);
-	x1 = std::max(x1, 0.0);
-	funcValue1 = effStressFunc(x1, effStressParams);
-      } else {
-	x2 += bracketFactor * (x1 - x2);
-	x2 = std::max(x2, 0.0);
-	funcValue2 = effStressFunc(x2, effStressParams);
-      } // else
+      x2 += bracketFactor * (x1 - x2);
+      x2 = std::max(x2, 0.0);
+      funcValue2 = effStressFunc(x2, effStressParams);
     } // else
     ++iteration;
   } // while
@@ -114,65 +112,64 @@
 
 // ----------------------------------------------------------------------
 // Find root using Newton's method with bisection.
-void
+double
 pylith::materials::EffectiveStress::_findEffStress(
-				     const double xx1,
-				     const double xx2,
-				     EffStressStruct effStressParams,
-				     effStressFuncType effStressFunc,
-				     effStressFuncDFuncType effStressFuncDFunc)
+				 const double x1,
+				 const double x2,
+				 const EffStressStruct& effStressParams,
+				 effStressFunc_fn_type effStressFunc,
+				 effStressFuncDFunc_fn_type effStressFuncDFunc)
 { // _findEffStress
   // Arbitrary number of iterations to find the root
   const int maxIterations = 100;
 
   // Desired accuracy for root. This is a bit arbitrary for now.
-  const double accuracy = 1.0e-10;
+  const double accuracy = 1.0e-6;
 
   /// Determine if root has already been found, or if root is not bracketed.
   // Otherwise, organize search so that effStressFunc(xLow) is less than zero.
-  double funcValueLow = effStressFunc(xx1, effStressParams);
-  double funcValueHigh = effStressFunc(xx2, effStressParams);
-  if (funcValueLow * funcValueHigh > 0.0)
-    throw std::runtime_error("Effective stress is not bracketed.");
+  double funcValueLow = effStressFunc(x1, effStressParams);
+  double funcValueHigh = effStressFunc(x2, effStressParams);
+  assert(funcValueLow * funcValueHigh <= 0.0);
 
   double effStress = 0.0;
   double xLow = 0.0;
   double xHigh = 0.0;
   bool converged = false;
 
-  if (std::abs(funcValueLow) < accuracy) {
-    effStress = xx1;
+  if (fabs(funcValueLow) < accuracy) {
+    effStress = x1;
     converged = true;
     return effStress;
-  } else if (std::abs(funcValueHigh) < accuracy) {
-    effStress = xx2;
+  } else if (fabs(funcValueHigh) < accuracy) {
+    effStress = x2;
     converged = true;
     return effStress;
   } else if (funcValueLow < 0.0) {
-    xLow = xx1;
-    xHigh = xx2;
+    xLow = x1;
+    xHigh = x2;
   } else {
-    xHigh = xx1;
-    xLow = xx2;
+    xHigh = x1;
+    xLow = x2;
   }
 
-  effStress = 0.5 * (xx1 + xx2);
-  double dxPrevious = std::abs(xx2 - xx1);
+  effStress = 0.5 * (x1 + x2);
+  double dxPrevious = fabs(x2 - x1);
   double dx = dxPrevious;
   double funcValue = 0.0;
   double funcDeriv = 0.0;
   double funcXHigh = 0.0;
   double funcXLow = 0.0;
-  effStressFuncDFunc(effStress, effStressParams, funcValue, funcDeriv);
+  effStressFuncDFunc(effStress, effStressParams, &funcValue, &funcDeriv);
   int iteration = 0;
 
-  while ((iteration < maxIterations) && (converged == false)) {
+  while (iteration < maxIterations) {
     funcXHigh = (effStress - xHigh) * funcDeriv - funcValue;
     funcXLow = (effStress - xLow) * funcDeriv - funcValue;
     // Use bisection if solution goes out of bounds or is not converging
     // fast enough.
     if ( (funcXHigh * funcXLow >= 0.0) ||
-	 (std::abs(2.0 * funcValue) > std::abs(dxPrevious * funcDeriv))) {
+	 (fabs(2.0 * funcValue) > fabs(dxPrevious * funcDeriv))) {
       dxPrevious = dx;
       dx = 0.5 * (xHigh - xLow);
       effStress = xLow + dx;
@@ -181,11 +178,11 @@
       dx = funcValue/funcDeriv;
       effStress = effStress - dx;
     } // else
-    if (std::abs(dx) < accuracy) {
+    if (fabs(dx) < accuracy) {
       converged = true;
-      return effStress;
+      break;
     } // if
-    effStressFuncDFunc(effStress, effStressParams, funcValue, funcDeriv);
+    effStressFuncDFunc(effStress, effStressParams, &funcValue, &funcDeriv);
     if (funcValue < 0.0) {
       xLow = effStress;
     } else {
@@ -198,6 +195,7 @@
     throw std::runtime_error("Cannot find root of effective stress function.");
 
   PetscLogFlops(5 + 15 * iteration);
+  return effStress;
 
 } // _findEffStress
 

Modified: short/3D/PyLith/trunk/libsrc/materials/EffectiveStress.hh
===================================================================
--- short/3D/PyLith/trunk/libsrc/materials/EffectiveStress.hh	2009-05-26 22:59:53 UTC (rev 15064)
+++ short/3D/PyLith/trunk/libsrc/materials/EffectiveStress.hh	2009-05-27 00:29:59 UTC (rev 15065)
@@ -54,17 +54,17 @@
   /// Member prototype for effStressFunc()
   typedef static double (*effStressFunc_fn_type)
   (const double,
-   EffStressStruct);
+   const EffStressStruct&);
 
   /// Member prototype for effStressDFunc()
   typedef static double (*effStressDFunc_fn_type)
   (const double,
-   EffStressStruct);
+   const EffStressStruct&);
   
   /// Member prototype for effStressFuncDFunc()
   typedef static void (*effStressFuncDFunc_fn_type)
   (const double,
-   EffStressStruct,
+   const EffStressStruct&,
    double*,
    double*);
 
@@ -96,7 +96,7 @@
   static
   double getEffStress(const double effStressInitialGuess,
 		      const double stressScale,
-		      EffStressStruct effStressParams,
+		      const EffStressStruct& effStressParams,
 		      effStressFunc_fn_type effStressFunc,
 		      effStressFuncDFunc_fn_type effStressFuncDFunc);
 
@@ -111,9 +111,10 @@
    * @param effStressFunc Function to compute effective stress only.
    *
    */
+  static
   void _bracketEffStress(double* px1,
 			 double* px2,
-			 EffStressStruct effStressParams,
+			 const EffStressStruct& effStressParams,
 			 effStressFunc_fn_type effStressFunc);
 
   /** Solve for effective stress using Newton's method with bisection.
@@ -126,9 +127,10 @@
    *
    * @returns Computed effective stress.
    */
+  static
   double _findEffStress(double xx1,
 			double xx2,
-			EffStressStruct effStressParams,
+			const EffStressStruct& effStressParams,
 			effStressFunc_fn_type effStressFunc,
 			effStressFuncDFunc_fn_type effStressFuncDFunc);
 

Modified: short/3D/PyLith/trunk/libsrc/materials/Makefile.am
===================================================================
--- short/3D/PyLith/trunk/libsrc/materials/Makefile.am	2009-05-26 22:59:53 UTC (rev 15064)
+++ short/3D/PyLith/trunk/libsrc/materials/Makefile.am	2009-05-27 00:29:59 UTC (rev 15065)
@@ -25,11 +25,14 @@
 	GenMaxwellIsotropic3D.icc \
 	MaxwellIsotropic3D.hh \
 	MaxwellIsotropic3D.icc \
+	PowerLaw3D.hh \
+	PowerLaw3D.icc \
 	Metadata.hh \
 	Metadata.icc \
 	Material.hh \
 	Material.icc \
 	ViscoelasticMaxwell.hh \
+	EffectiveStress.hh \
 	materialsfwd.hh
 
 

Modified: short/3D/PyLith/trunk/libsrc/materials/PowerLaw3D.cc
===================================================================
--- short/3D/PyLith/trunk/libsrc/materials/PowerLaw3D.cc	2009-05-26 22:59:53 UTC (rev 15064)
+++ short/3D/PyLith/trunk/libsrc/materials/PowerLaw3D.cc	2009-05-27 00:29:59 UTC (rev 15065)
@@ -638,11 +638,10 @@
 // ----------------------------------------------------------------------
 // Effective stress function that computes effective stress function only
 // (no derivative).
-static
 double
 pylith::materials::PowerLaw3D::effStressFunc(
      const double effStressTpdt,
-     pylith::materials::EffectiveStress::EffStressStruct& effStressParams)
+     const EffectiveStress::EffStressStruct& effStressParams)
 { // effStressFunc
   double ae = effStressParams.ae;
   double b = effStressParams.b;
@@ -667,11 +666,10 @@
 // ----------------------------------------------------------------------
 // Effective stress function that computes effective stress function
 // derivative only (no function value).
-static
 double
 pylith::materials::PowerLaw3D::effStressDFunc(
      const double effStressTpdt,
-     pylith::materials::EffectiveStress::EffStressStruct& effStressParams)
+     const EffectiveStress::EffStressStruct& effStressParams)
 { // effStressDFunc
   double ae = effStressParams.ae;
   double c = effStressParams.c;
@@ -699,11 +697,10 @@
 // ----------------------------------------------------------------------
 // Effective stress function that computes effective stress function
 // and derivative.
-static
 void
 pylith::materials::PowerLaw3D::effStressFuncDFunc(
      const double effStressTpdt,
-     pylith::materials::EffectiveStress::EffStressStruct& effStressParams,
+     const EffectiveStress::EffStressStruct& effStressParams,
      double* py,
      double* pdy)
 { // effStressFuncDFunc

Modified: short/3D/PyLith/trunk/libsrc/materials/PowerLaw3D.hh
===================================================================
--- short/3D/PyLith/trunk/libsrc/materials/PowerLaw3D.hh	2009-05-26 22:59:53 UTC (rev 15064)
+++ short/3D/PyLith/trunk/libsrc/materials/PowerLaw3D.hh	2009-05-27 00:29:59 UTC (rev 15065)
@@ -64,7 +64,7 @@
    */
   static double effStressFunc(
     const double effStressTpdt,
-    pylith::materials::EffectiveStress::EffStressStruct& effStressParams);
+    const EffectiveStress::EffStressStruct& effStressParams);
 
   /** Compute effective stress function derivative.
    *
@@ -75,7 +75,7 @@
    */
   static double effStressDFunc(
     const double effStressTpdt,
-    pylith::materials::EffectiveStress::EffStressStruct& effStressParams);
+    const EffectiveStress::EffStressStruct& effStressParams);
 
   /** Compute effective stress function and derivative.
    *
@@ -87,7 +87,7 @@
    */
   static void effStressFuncDFunc(
     const double effStressTpdt,
-    pylith::materials::EffectiveStress::EffStressStruct& effStressParams,
+    const EffectiveStress::EffStressStruct& effStressParams,
     double* py,
     double* pdy);
 



More information about the CIG-COMMITS mailing list