[cig-commits] r14975 - short/3D/PyLith/trunk/libsrc/materials
willic3 at geodynamics.org
willic3 at geodynamics.org
Mon May 11 06:25:16 PDT 2009
Author: willic3
Date: 2009-05-11 06:25:15 -0700 (Mon, 11 May 2009)
New Revision: 14975
Added:
short/3D/PyLith/trunk/libsrc/materials/EffectiveStress.hh
Modified:
short/3D/PyLith/trunk/libsrc/materials/EffectiveStress.cc
short/3D/PyLith/trunk/libsrc/materials/PowerLaw3D.cc
short/3D/PyLith/trunk/libsrc/materials/PowerLaw3D.hh
Log:
More work on PowerLaw3D. The functions are pretty much finished.
Just need to fix how arguments are passed, etc.
Modified: short/3D/PyLith/trunk/libsrc/materials/EffectiveStress.cc
===================================================================
--- short/3D/PyLith/trunk/libsrc/materials/EffectiveStress.cc 2009-05-11 04:52:20 UTC (rev 14974)
+++ short/3D/PyLith/trunk/libsrc/materials/EffectiveStress.cc 2009-05-11 13:25:15 UTC (rev 14975)
@@ -26,10 +26,10 @@
// actual initial guess is zero.
double
pylith::materials::EffectiveStress::getEffStress(
- const double effStressInitialGuess,
- const double* effStressParams,
- ***effStressFunc,
- ***effStressFuncDFunc)
+ const double effStressInitialGuess,
+ const EffStressStruct& effStressParams,
+ effStressFuncType* effStressFunc,
+ effStressFuncDFuncType* effStressFuncDFunc)
{ // getEffStress
// Check parameters
if (effStressInitialGuess < 0.0)
@@ -42,16 +42,16 @@
x1 = effStressInitialGuess - 0.5 * effStressInitialGuess;
x2 = effStressInitialGuess + 0.5 * effStressInitialGuess;
} else {
- x1 = effStressParams[0] - 0.5 * effStressParams[0];
- x2 = effStressParams[0] + 0.5 * effStressParams[0];
+ x1 = effStressParams.stressScale - 0.5 * effStressParams.stressScale;
+ x2 = effStressParams.stressScale + 0.5 * effStressParams.stressScale;
} // else
PetscLogFlops(4);
_bracketEffStress(x1, x2, effStressParams, effStressFunc);
// Find effective stress using Newton's method with bisection.
- const double effStress = _findEffStress(x1,
- x2,
+ const double effStress = _findEffStress(&x1,
+ &x2,
effStressParams,
effStressFunc,
effStressFuncDFunc);
@@ -64,16 +64,18 @@
// Bracket effective stress.
void
pylith::materials::EffectiveStress::_bracketEffStress(
- double* const x1,
- double* const x2,
- const double* effStressParams,
- ??? effStressFunc)
+ double* px1,
+ double* px2,
+ const EffStressStruct& effStressParams,
+ effStressFuncType* effStressFunc)
{ // _bracketEffStress
// Arbitrary number of iterations to bracket the root
const int maxIterations = 50;
// Arbitrary factor by which to increase the brackets.
const double bracketFactor = 1.6;
+ double x1 = *px1;
+ double x2 = *px2;
double funcValue1 = effStressFunc(x1, effStressParams);
double funcValue2 = effStressFunc(x2, effStressParams);
@@ -98,6 +100,9 @@
++iteration;
} // while
+ *px1 = x1;
+ *px2 = x2;
+
PetscLogFlops(5 * iteration);
if (bracketed == false)
throw std::runtime_error("Unable to bracket effective stress.");
@@ -108,11 +113,11 @@
// Find root using Newton's method with bisection.
void
pylith::materials::EffectiveStress::_findEffStress(
- double* const x1,
- double* const x2,
- const double* effStressParams,
- ??? effStressFunc,
- ??? effStressFuncDFunc)
+ const double x1,
+ const double x2,
+ const EffStressStruct& effStressParams,
+ effStressFuncType* effStressFunc,
+ effStressFuncDFuncType* effStressFuncDFunc)
{ // _findEffStress
// Arbitrary number of iterations to find the root
const int maxIterations = 100;
@@ -130,12 +135,15 @@
double effStress = 0.0;
double xLow = 0.0;
double xHigh = 0.0;
+ bool converged = false;
if (std::abs(funcValueLow) < accuracy) {
effStress = x1;
+ converged = true;
return effStress;
} else if (std::abs(funcValueHigh) < accuracy) {
effStress = x2;
+ converged = true;
return effStress;
} else if (funcValueLow < 0.0) {
xLow = x1;
@@ -150,31 +158,45 @@
double dx = dxPrevious;
double funcValue = 0.0;
double funcDeriv = 0.0;
+ double funcXHigh = 0.0;
+ double funcXLow = 0.0;
effStressFuncDFunc(effStress, effStressParams, funcValue, funcDeriv);
int iteration = 0;
- bool converged = false;
- // ******* finish fixing through here ***********
while ((iteration < maxIterations) && (converged == false)) {
- if ((funcValue1 * funcValue2) < 0.0) {
- bracketed = true;
+ 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))) {
+ dxPrevious = dx;
+ dx = 0.5 * (xHigh - xLow);
+ effStress = xLow + dx;
} 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
+ dxPrevious = dx;
+ dx = funcValue/funcDeriv;
+ effStress = effStress - dx;
} // else
+ if (std::abs(dx) < accuracy) {
+ converged = true;
+ return effStress;
+ } // if
+ effStressFuncDFunc(effStress, effStressParams, funcValue, funcDeriv);
+ if (funcValue < 0.0) {
+ xLow = effStress;
+ } else {
+ xHigh = effStress;
+ } // else
+ ++iteration;
} // while
- if (bracketed == false)
- throw std::runtime_error("Unable to bracket effective stress.");
+ if (converged == false)
+ throw std::runtime_error("Cannot find root of effective stress function.");
-} // _bracketEffStress
+ PetscLogFlops(5 + 15 * iteration);
+} // _findEffStress
+
// End of file
Added: short/3D/PyLith/trunk/libsrc/materials/EffectiveStress.hh
===================================================================
--- short/3D/PyLith/trunk/libsrc/materials/EffectiveStress.hh (rev 0)
+++ short/3D/PyLith/trunk/libsrc/materials/EffectiveStress.hh 2009-05-11 13:25:15 UTC (rev 14975)
@@ -0,0 +1,102 @@
+// -*- C++ -*-
+//
+// ----------------------------------------------------------------------
+//
+// Brad T. Aagaard
+// U.S. Geological Survey
+//
+// {LicenseText}
+//
+// ----------------------------------------------------------------------
+//
+
+/** @file libsrc/materials/EffectiveStress.hh
+ *
+ * @brief C++ EffectiveStress object.
+ *
+ * This class contains bracketing and root-finding functions for materials that
+ * use an effective stress formulation.
+ */
+
+#if !defined(pylith_materials_effectivestress_hh)
+#define pylith_materials_effectivestress_hh
+
+/// Namespace for pylith package
+namespace pylith {
+ namespace materials {
+ class EffectiveStress;
+ } // materials
+
+} // pylith
+
+/// C++ abstract base class for Material object.
+class pylith::materials::EffectiveStress
+{ // class EffectiveStress
+
+ // PUBLIC STRUCTS /////////////////////////////////////////////////////
+public :
+
+ struct EffStressStruct {
+ double stressScale;
+ double ae;
+ double b;
+ double c;
+ double d;
+ double alpha;
+ double dt;
+ double effStressT;
+ double powerLawExp;
+ double viscosityCoeff;
+ };
+
+ // PUBLIC METHODS /////////////////////////////////////////////////////
+public :
+
+ /** Get effective stress from initial guess.
+ *
+ * @param effStressInitialGuess Initial guess for effective stress.
+ * @param effStressParams Parameters used in computing effective stress.
+ * @param effStressFunc Function to compute effective stress only.
+ * @param effStressFuncDFunc Function to compute effective stress and derivative.
+ *
+ * @returns Computed effective stress.
+ */
+ static double getEffStress(const double effStressInitialGuess,
+ const EffStressStruct& effStressParams,
+ effStressFuncType* effStressFunc,
+ effStressFuncDFuncType* effStressFuncDFunc);
+
+ // PRIVATE METHODS /////////////////////////////////////////////////////
+private :
+
+ /** Bracket effective stress.
+ *
+ * @param x1 Initial guess for first bracket.
+ * @param x2 Initial guess for second bracket.
+ * @param effStressParams Parameters used in computing effective stress.
+ * @param effStressFunc Function to compute effective stress only.
+ *
+ */
+ void bracketEffStress(double x1,
+ double x2,
+ const double effStressParams,
+ static double &effStressFunc);
+
+ /** Solve for effective stress using Newton's method with bisection.
+ *
+ * @param x1 Initial guess for first bracket.
+ * @param x2 Initial guess for second bracket.
+ * @param effStressParams Parameters used in computing effective stress.
+ * @param effStressFunc Function to compute effective stress only.
+ *
+ */
+ void bracketEffStress(double x1,
+ double x2,
+ const double effStressParams,
+ static double &effStressFunc);
+}; // class EffectiveStress
+
+#endif // pylith_materials_effectivestress_hh
+
+
+// End of file
Modified: short/3D/PyLith/trunk/libsrc/materials/PowerLaw3D.cc
===================================================================
--- short/3D/PyLith/trunk/libsrc/materials/PowerLaw3D.cc 2009-05-11 04:52:20 UTC (rev 14974)
+++ short/3D/PyLith/trunk/libsrc/materials/PowerLaw3D.cc 2009-05-11 13:25:15 UTC (rev 14975)
@@ -568,27 +568,26 @@
const double stressScale = mu;
PetscLogFlops(92);
- // Put parameters into a vector and call root-finding algorithm.
- // This could also be a struct.
- const double effStressParams[] = {stressScale,
- ae,
- b,
- c,
- d,
- alpha,
- _dt,
- effectiveStressT,
- powerLawExp,
- viscosityCoeff};
- // I think the PETSc root-finding procedure is too involved for what we want
- // here. I would like the function to work something like:
+ // Put parameters into a struct and call root-finding algorithm.
+ _effStressParams.stressScale = stressScale;
+ _effStressParams.ae = ae;
+ _effStressParams.b = b;
+ _effStressParams.c = c;
+ _effStressParams.d = d;
+ _effStressParams.alpha = alpha;
+ _effStressParams.dt = _dt;
+ _effStressParams.effStressT = effStressT;
+ _effStressParams.powerLawExp = powerLawExp;
+ _effStressParams.viscosityCoeff = viscosityCoeff;
+
const double effStressInitialGuess = effStressT;
+
double effStressTpdt =
EffectiveStress::getEffStress(
effStressInitialGuess,
- effStressParams,
- pylith::materials::PowerLaw3D::_effStressFunc,
- pylith::materials::PowerLaw3D::_effStressFuncDFunc);
+ &_effStressParams,
+ *pylith::materials::PowerLaw3D::effStressFunc,
+ *pylith::materials::PowerLaw3D::effStressFuncDFunc);
// Compute stresses from effective stress.
const double effStressTau = (1.0 - alpha) * effStressT +
@@ -620,18 +619,19 @@
// Effective stress function that computes effective stress function only
// (no derivative).
double
-pylith::materials::PowerLaw3D::effStressFunc(double effStressTpdt,
- double *params)
+pylith::materials::PowerLaw3D::effStressFunc(
+ const double effStressTpdt,
+ const EffStressStruct& effStressParams)
{ // effStressFunc
- double ae = params[0];
- double b = params[1];
- double c = params[2];
- double d = params[3];
- double alpha = params[4];
- double dt = params[5];
- double effStressT = params[6];
- double powerLawExp = params[7];
- double viscosityCoeff = params[8];
+ double ae = effStressParams.ae;
+ double b = effStressParams.b;
+ double c = effStressParams.c;
+ double d = effStressParams.d;
+ double alpha = effStressParams.alpha;
+ double dt = effStressParams.dt;
+ double effStressT = effStressParams.effStressT;
+ double powerLawExp = effStressParams.powerLawExp;
+ double viscosityCoeff = effStressParams.viscosityCoeff;
double factor1 = 1.0-alpha;
double effStressTau = factor1 * effStressT + alpha * effStressTpdt;
double gammaTau = 0.5 * ((effStressTau/viscosityCoeff)^
@@ -647,17 +647,18 @@
// Effective stress function that computes effective stress function
// derivative only (no function value).
double
-pylith::materials::PowerLaw3D::effStressDFunc(double effStressTpdt,
- double *params)
+pylith::materials::PowerLaw3D::effStressDFunc(
+ const double effStressTpdt,
+ const EffStressStruct& effStressParams)
{ // effStressDFunc
- double ae = params[0];
- double c = params[2];
- double d = params[3];
- double alpha = params[4];
- double dt = params[5];
- double effStressT = params[6];
- double powerLawExp = params[7];
- double viscosityCoeff = params[8];
+ double ae = effStressParams.ae;
+ double c = effStressParams.c;
+ double d = effStressParams.d;
+ double alpha = effStressParams.alpha;
+ double dt = effStressParams.dt;
+ double effStressT = effStressParams.effStressT;
+ double powerLawExp = effStressParams.powerLawExp;
+ double viscosityCoeff = effStressParams.viscosityCoeff;
double factor1 = 1.0-alpha;
double effStressTau = factor1 * effStressT + alpha * effStressTpdt;
double gammaTau = 0.5 * ((effStressTau/viscosityCoeff)^
@@ -677,20 +678,24 @@
// Effective stress function that computes effective stress function
// and derivative.
void
-pylith::materials::PowerLaw3D::effStressFuncDFunc(double effStressTpdt,
- double *params,
- double *y,
- double *dy)
-{ // effStressFunc
- double ae = params[0];
- double b = params[1];
- double c = params[2];
- double d = params[3];
- double alpha = params[4];
- double dt = params[5];
- double effStressT = params[6];
- double powerLawExp = params[7];
- double viscosityCoeff = params[8];
+pylith::materials::PowerLaw3D::effStressFuncDFunc(
+ const double effStressTpdt,
+ const EffStressStruct& effStressParams,
+ double* py,
+ double* pdy)
+{ // effStressFuncDFunc
+ double y = *py;
+ double dy = *pdy;
+
+ double ae = effStressParams.ae;
+ double b = effStressParams.b;
+ double c = effStressParams.c;
+ double d = effStressParams.d;
+ double alpha = effStressParams.alpha;
+ double dt = effStressParams.dt;
+ double effStressT = effStressParams.effStressT;
+ double powerLawExp = effStressParams.powerLawExp;
+ double viscosityCoeff = effStressParams.viscosityCoeff;
double factor1 = 1.0-alpha;
double effStressTau = factor1 * effStressT + alpha * effStressTpdt;
double gammaTau = 0.5 * ((effStressTau/viscosityCoeff)^
@@ -699,11 +704,18 @@
((effStressTau/viscosityCoeff) ^ (powerLawExp - 2.0))/
(viscosityCoeff * viscosityCoeff);
double a = ae + alpha * dt * gammaTau;
- double *y = a * a * effStressTpdt * effStressTpdt - b +
- c * gammaTau - d * d * gammaTau * gammaTau;
- double *dy = 2.0 * a * a * effStressTpdt + dGammaTau *
+ y = a * a * effStressTpdt * effStressTpdt -
+ b +
+ c * gammaTau -
+ d * d * gammaTau * gammaTau;
+ dy = 2.0 * a * a * effStressTpdt +
+ dGammaTau *
(2.0 * a * alpha * dt * effStressTpdt * effStressTpdt +
c - 2.0 * d * d * gammaTau);
+
+ *py = y;
+ *pdy = dy;
+
PetscLogFlops(46);
} // effStressFuncDFunc
@@ -962,27 +974,26 @@
const double stressScale = mu;
PetscLogFlops(92);
- // Put parameters into a vector and call root-finding algorithm.
- // This could also be a struct.
- const double effStressParams[] = {stressScale,
- ae,
- b,
- c,
- d,
- alpha,
- _dt,
- effectiveStressT,
- powerLawExp,
- viscosityCoeff};
- // I think the PETSc root-finding procedure is too involved for what we want
- // here. I would like the function to work something like:
+ // Put parameters into a struct and call root-finding algorithm.
+ _effStressParams.stressScale = stressScale;
+ _effStressParams.ae = ae;
+ _effStressParams.b = b;
+ _effStressParams.c = c;
+ _effStressParams.d = d;
+ _effStressParams.alpha = alpha;
+ _effStressParams.dt = _dt;
+ _effStressParams.effStressT = effStressT;
+ _effStressParams.powerLawExp = powerLawExp;
+ _effStressParams.viscosityCoeff = viscosityCoeff;
+
const double effStressInitialGuess = effStressT;
+
const double effStressTpdt =
EffectiveStress::getEffStress(
effStressInitialGuess,
- effStressParams,
- pylith::materials::PowerLaw3D::_effStressFunc,
- pylith::materials::PowerLaw3D::_effStressFuncDFunc);
+ &_effStressParams,
+ *pylith::materials::PowerLaw3D::effStressFunc,
+ *pylith::materials::PowerLaw3D::effStressFuncDFunc);
// Compute quantities at intermediate time tau used to compute values at
// end of time step.
@@ -1275,27 +1286,26 @@
const double stressScale = mu;
PetscLogFlops(92);
- // Put parameters into a vector and call root-finding algorithm.
- // This could also be a struct.
- const double effStressParams[] = {stressScale,
- ae,
- b,
- c,
- d,
- alpha,
- _dt,
- effectiveStressT,
- powerLawExp,
- viscosityCoeff};
- // I think the PETSc root-finding procedure is too involved for what we want
- // here. I would like the function to work something like:
+ // Put parameters into a struct and call root-finding algorithm.
+ _effStressParams.stressScale = stressScale;
+ _effStressParams.ae = ae;
+ _effStressParams.b = b;
+ _effStressParams.c = c;
+ _effStressParams.d = d;
+ _effStressParams.alpha = alpha;
+ _effStressParams.dt = _dt;
+ _effStressParams.effStressT = effStressT;
+ _effStressParams.powerLawExp = powerLawExp;
+ _effStressParams.viscosityCoeff = viscosityCoeff;
+
const double effStressInitialGuess = effStressT;
+
double effStressTpdt =
EffectiveStress::getEffStress(
effStressInitialGuess,
- effStressParams,
- pylith::materials::PowerLaw3D::_effStressFunc,
- pylith::materials::PowerLaw3D::_effStressFuncDFunc);
+ &_effStressParams,
+ *pylith::materials::PowerLaw3D::effStressFunc,
+ *pylith::materials::PowerLaw3D::effStressFuncDFunc);
// Compute stress and viscous strain and update appropriate state variables.
const double effStressTau = (1.0 - alpha) * effStressT +
Modified: short/3D/PyLith/trunk/libsrc/materials/PowerLaw3D.hh
===================================================================
--- short/3D/PyLith/trunk/libsrc/materials/PowerLaw3D.hh 2009-05-11 04:52:20 UTC (rev 14974)
+++ short/3D/PyLith/trunk/libsrc/materials/PowerLaw3D.hh 2009-05-11 13:25:15 UTC (rev 14975)
@@ -26,16 +26,9 @@
#if !defined(pylith_materials_powerlaw3d_hh)
#define pylith_materials_powerlaw3d_hh
-#include "ElasticMaterial.hh"
+#include "ElasticMaterial.hh" // ISA ElasticMaterial
+#include "EffectiveStress.hh" // Uses EffectiveStress
-/// Namespace for pylith package
-namespace pylith {
- namespace materials {
- class PowerLaw3D;
- class TestPowerLaw3D; // unit testing
- } // materials
-} // pylith
-
/// 3-D, isotropic, linear Maxwell viscoelastic material.
class pylith::materials::PowerLaw3D : public ElasticMaterial
{ // class PowerLaw3D
@@ -444,6 +437,9 @@
// PRIVATE MEMBERS ////////////////////////////////////////////////////
private :
+ /// Structure to hold parameters for effective stress computation.
+ pylith::materials::EffectiveStress::EffStressStruct _effStressParams;
+
/// Method to use for _calcElasticConsts().
calcElasticConsts_fn_type _calcElasticConstsFn;
More information about the CIG-COMMITS
mailing list