[cig-commits] r15211 - in short/3D/PyLith/trunk: . libsrc/materials unittests/libtests/materials unittests/libtests/materials/data
brad at geodynamics.org
brad at geodynamics.org
Fri Jun 12 11:02:21 PDT 2009
Author: brad
Date: 2009-06-12 11:02:21 -0700 (Fri, 12 Jun 2009)
New Revision: 15211
Modified:
short/3D/PyLith/trunk/TODO
short/3D/PyLith/trunk/libsrc/materials/EffectiveStress.icc
short/3D/PyLith/trunk/libsrc/materials/PowerLaw3D.cc
short/3D/PyLith/trunk/unittests/libtests/materials/TestEffectiveStress.cc
short/3D/PyLith/trunk/unittests/libtests/materials/data/PowerLaw3DTimeDep.py
Log:
Fixed a few more power law bugs. Effective stress calculation gives strain increment, not total strain. Cleaned up effective stress calculation.
Modified: short/3D/PyLith/trunk/TODO
===================================================================
--- short/3D/PyLith/trunk/TODO 2009-06-12 15:54:47 UTC (rev 15210)
+++ short/3D/PyLith/trunk/TODO 2009-06-12 18:02:21 UTC (rev 15211)
@@ -33,12 +33,10 @@
ArbitratySlipFn (final slip, start time)
full-scale testing
2d/quad4
- axiialtract
+ axialtract
sheartract
dislocation2
2d/tri3
- axialdisp
- sheardisp
dislocation
dislocation2
cleanup
Modified: short/3D/PyLith/trunk/libsrc/materials/EffectiveStress.icc
===================================================================
--- short/3D/PyLith/trunk/libsrc/materials/EffectiveStress.icc 2009-06-12 15:54:47 UTC (rev 15210)
+++ short/3D/PyLith/trunk/libsrc/materials/EffectiveStress.icc 2009-06-12 18:02:21 UTC (rev 15211)
@@ -63,7 +63,7 @@
const int maxIterations = 50;
// Arbitrary factor by which to increase the brackets.
- const double bracketFactor = 1.6;
+ const double bracketFactor = 2;
double x1 = *px1;
double x2 = *px2;
@@ -110,10 +110,9 @@
const int maxIterations = 100;
// Desired accuracy for root. This is a bit arbitrary for now.
- const double accuracy = 1.0e-6;
+ const double accuracy = 1.0e-20;
- /// Determine if root has already been found, or if root is not bracketed.
- // Otherwise, organize search so that effStressFunc(xLow) is less than zero.
+ // Organize search so that effStressFunc(xLow) is less than zero.
double funcValueLow = material->effStressFunc(x1);
double funcValueHigh = material->effStressFunc(x2);
assert(funcValueLow * funcValueHigh <= 0.0);
@@ -123,21 +122,13 @@
double xHigh = 0.0;
bool converged = false;
- if (fabs(funcValueLow) < accuracy) {
- effStress = x1;
- converged = true;
- return effStress;
- } else if (fabs(funcValueHigh) < accuracy) {
- effStress = x2;
- converged = true;
- return effStress;
- } else if (funcValueLow < 0.0) {
+ if (funcValueLow < 0.0) {
xLow = x1;
xHigh = x2;
} else {
- xHigh = x1;
xLow = x2;
- }
+ xHigh = x1;
+ } // if/else
effStress = 0.5 * (x1 + x2);
double dxPrevious = fabs(x2 - x1);
@@ -152,6 +143,10 @@
while (iteration < maxIterations) {
funcXHigh = (effStress - xHigh) * funcDeriv - funcValue;
funcXLow = (effStress - xLow) * funcDeriv - funcValue;
+ if (fabs(funcXLow) < accuracy || fabs(funcXHigh) < accuracy) {
+ converged = true;
+ break;
+ } // if
// Use bisection if solution goes out of bounds or is not converging
// fast enough.
if ( (funcXHigh * funcXLow >= 0.0) ||
@@ -164,10 +159,6 @@
dx = funcValue / funcDeriv;
effStress = effStress - dx;
} // else
- if (fabs(dx) < accuracy) {
- converged = true;
- break;
- } // if
material->effStressFuncDerivFunc(&funcValue, &funcDeriv, effStress);
if (funcValue < 0.0) {
xLow = effStress;
Modified: short/3D/PyLith/trunk/libsrc/materials/PowerLaw3D.cc
===================================================================
--- short/3D/PyLith/trunk/libsrc/materials/PowerLaw3D.cc 2009-06-12 15:54:47 UTC (rev 15210)
+++ short/3D/PyLith/trunk/libsrc/materials/PowerLaw3D.cc 2009-06-12 18:02:21 UTC (rev 15211)
@@ -658,7 +658,9 @@
const double a = ae + alpha * dt * gammaTau;
const double y = a * a * effStressTpdt * effStressTpdt - b +
c * gammaTau - d * d * gammaTau * gammaTau;
+
PetscLogFlops(21);
+
return y;
} // effStressFunc
@@ -1360,13 +1362,13 @@
stateVars[s_stress+iComp] = devStressTpdt + diag[iComp] *
(meanStressTpdt + meanStressInitial);
devStressTau = (1.0 - alpha) * devStressT[iComp] + alpha * devStressTpdt;
- stateVars[s_viscousStrain+iComp] = _dt * gammaTau * devStressTau;
+ stateVars[s_viscousStrain+iComp] += _dt * gammaTau * devStressTau;
} // for
_needNewJacobian = true;
PetscLogFlops(14 + _tensorSize * 14);
-} // _updatePropertiesViscoelastic
+} // _updateStateVarsViscoelastic
// ----------------------------------------------------------------------
// Compute scalar product of two tensors.
Modified: short/3D/PyLith/trunk/unittests/libtests/materials/TestEffectiveStress.cc
===================================================================
--- short/3D/PyLith/trunk/unittests/libtests/materials/TestEffectiveStress.cc 2009-06-12 15:54:47 UTC (rev 15210)
+++ short/3D/PyLith/trunk/unittests/libtests/materials/TestEffectiveStress.cc 2009-06-12 18:02:21 UTC (rev 15211)
@@ -55,10 +55,10 @@
Quadratic(void) {};
~Quadratic(void) {};
double effStressFunc(const double x) {
- return 1.0e-2 - pow(x - 1.0e-3, 2);
+ return 1.0e+5 - 1.0/9.0e+3 * pow(x + 2.0e+4, 2);
};
double effStressDerivFunc(const double x) {
- return -2*(x-1.0e-03);
+ return -2*1.0/9.0e+3*(x+2.0e+4);
};
double effStressFuncDerivFunc(double* f,
double* df,
@@ -122,7 +122,7 @@
void
pylith::materials::TestEffectiveStress::testCalculateQuadratic(void)
{ // testCalculateQuadratic
- const double valueE = 0.101;
+ const double valueE = 1.0e+04;
_EffectiveStress::Quadratic material;
Modified: short/3D/PyLith/trunk/unittests/libtests/materials/data/PowerLaw3DTimeDep.py
===================================================================
--- short/3D/PyLith/trunk/unittests/libtests/materials/data/PowerLaw3DTimeDep.py 2009-06-12 15:54:47 UTC (rev 15210)
+++ short/3D/PyLith/trunk/unittests/libtests/materials/data/PowerLaw3DTimeDep.py 2009-06-12 18:02:21 UTC (rev 15211)
@@ -377,7 +377,7 @@
self.dt, effStressT, powerLawExpV, viscosityCoeffV)
# Find the root using Brent's method (from scipy)
- rootTolerance = 1.0e-12
+ rootTolerance = 1.0e-14
effStressTpdt = scipy.optimize.brentq(self._effStressFunc, x1, x2,
args=(ae, b, c, d, self.alpha,
self.dt, effStressT,
More information about the CIG-COMMITS
mailing list