[cig-commits] [commit] knepley/upgrade-petsc-interface: Fixed stable time step for power-law. (3f1dca7)

cig_noreply at geodynamics.org cig_noreply at geodynamics.org
Wed Jan 15 18:58:23 PST 2014


Repository : ssh://geoshell/pylith

On branch  : knepley/upgrade-petsc-interface
Link       : https://github.com/geodynamics/pylith/compare/780dd6347510363e3138470cbcacc3db5d61142b...3f1dca7cdf8548b151adb1185932fd7966545e67

>---------------------------------------------------------------

commit 3f1dca7cdf8548b151adb1185932fd7966545e67
Author: Charles Williams <C.Williams at gns.cri.nz>
Date:   Thu Jan 16 15:58:10 2014 +1300

    Fixed stable time step for power-law.


>---------------------------------------------------------------

3f1dca7cdf8548b151adb1185932fd7966545e67
 libsrc/pylith/materials/PowerLaw3D.cc                              | 7 +++----
 libsrc/pylith/materials/PowerLawPlaneStrain.cc                     | 7 +++----
 unittests/libtests/materials/data/PowerLaw3DElastic.py             | 6 +++---
 unittests/libtests/materials/data/PowerLaw3DElasticData.cc         | 2 +-
 unittests/libtests/materials/data/PowerLaw3DTimeDep.py             | 6 +++---
 unittests/libtests/materials/data/PowerLaw3DTimeDepData.cc         | 2 +-
 unittests/libtests/materials/data/PowerLawPlaneStrainElastic.py    | 6 +++---
 .../libtests/materials/data/PowerLawPlaneStrainElasticData.cc      | 2 +-
 unittests/libtests/materials/data/PowerLawPlaneStrainTimeDep.py    | 6 +++---
 .../libtests/materials/data/PowerLawPlaneStrainTimeDepData.cc      | 2 +-
 10 files changed, 22 insertions(+), 24 deletions(-)

diff --git a/libsrc/pylith/materials/PowerLaw3D.cc b/libsrc/pylith/materials/PowerLaw3D.cc
index abb05b1..b2f3ac8 100644
--- a/libsrc/pylith/materials/PowerLaw3D.cc
+++ b/libsrc/pylith/materials/PowerLaw3D.cc
@@ -421,9 +421,8 @@ pylith::materials::PowerLaw3D::_stableTimeStepImplicit(
   if (effStress <= 0.0) {
     dtTest = pylith::PYLITH_MAXSCALAR;
   } else {
-    dtTest = 0.05 *
-    pow((referenceStress/effStress), (powerLawExp - 1.0)) *
-    (referenceStress/mu)/referenceStrainRate;
+    dtTest = pow((referenceStress/effStress), (powerLawExp - 1.0)) *
+      (referenceStress/mu)/(referenceStrainRate * 30.0);
   } //else
   const PylithScalar dtStable = dtTest;
 
@@ -431,7 +430,7 @@ pylith::materials::PowerLaw3D::_stableTimeStepImplicit(
   PylithScalar maxwellTime = 10.0 * dtStable;
   std::cout << "Maxwell time:  " << maxwellTime << std::endl;
 #endif
-  PetscLogFlops(21);
+  PetscLogFlops(27);
   return dtStable;
 } // _stableTimeStepImplicit
 
diff --git a/libsrc/pylith/materials/PowerLawPlaneStrain.cc b/libsrc/pylith/materials/PowerLawPlaneStrain.cc
index b17e6e1..e0e008c 100644
--- a/libsrc/pylith/materials/PowerLawPlaneStrain.cc
+++ b/libsrc/pylith/materials/PowerLawPlaneStrain.cc
@@ -421,9 +421,8 @@ pylith::materials::PowerLawPlaneStrain::_stableTimeStepImplicit(
   if (effStress <= 0.0) {
     dtTest = pylith::PYLITH_MAXSCALAR;
   } else {
-    dtTest = 0.05 *
-    pow((referenceStress/effStress), (powerLawExp - 1.0)) *
-    (referenceStress/mu)/referenceStrainRate;
+    dtTest = pow((referenceStress/effStress), (powerLawExp - 1.0)) *
+      (referenceStress/mu)/(referenceStrainRate * 30.0);
   } //else
   const PylithScalar dtStable = dtTest;
 
@@ -431,7 +430,7 @@ pylith::materials::PowerLawPlaneStrain::_stableTimeStepImplicit(
   PylithScalar maxwellTime = 10.0 * dtStable;
   std::cout << "Maxwell time:  " << maxwellTime << std::endl;
 #endif
-  PetscLogFlops(22);
+  PetscLogFlops(23);
   return dtStable;
 } // _stableTimeStepImplicit
 
diff --git a/unittests/libtests/materials/data/PowerLaw3DElastic.py b/unittests/libtests/materials/data/PowerLaw3DElastic.py
index 504565a..6af5f55 100644
--- a/unittests/libtests/materials/data/PowerLaw3DElastic.py
+++ b/unittests/libtests/materials/data/PowerLaw3DElastic.py
@@ -192,7 +192,7 @@ class PowerLaw3DElastic(ElasticMaterialApp):
                                           stressUpdated[1,:]],
                                          dtype=numpy.float64)
 
-    self.dtStableImplicit = 0.1*min(maxwellTimeA, maxwellTimeB)
+    self.dtStableImplicit = 0.2 * min(maxwellTimeA, maxwellTimeB)
     self.dtStableExplicit = 1000.0 / vpA
 
     return
@@ -228,8 +228,8 @@ class PowerLaw3DElastic(ElasticMaterialApp):
     effStress = (0.5 * devStressProd)**0.5
     maxwellTime = 1.0
     if (effStress != 0.0):
-      maxwellTime = 0.5 * (refStress/effStress)**(powerLawExponent - 1.0) * \
-                    (refStress/mu)/refStrainRate
+      maxwellTime = (refStress/effStress)**(powerLawExponent - 1.0) * \
+                    (refStress/mu)/(refStrainRate * 6.0)
 
     return maxwellTime
 
diff --git a/unittests/libtests/materials/data/PowerLaw3DElasticData.cc b/unittests/libtests/materials/data/PowerLaw3DElasticData.cc
index b893477..e136ef7 100644
--- a/unittests/libtests/materials/data/PowerLaw3DElasticData.cc
+++ b/unittests/libtests/materials/data/PowerLaw3DElasticData.cc
@@ -45,7 +45,7 @@ const PylithScalar pylith::materials::PowerLaw3DElasticData::_pressureScale =
 
 const PylithScalar pylith::materials::PowerLaw3DElasticData::_densityScale =   2.25000000e+04;
 
-const PylithScalar pylith::materials::PowerLaw3DElasticData::_dtStableImplicit =   4.44444444e+06;
+const PylithScalar pylith::materials::PowerLaw3DElasticData::_dtStableImplicit =   2.96296296e+06;
 
 const PylithScalar pylith::materials::PowerLaw3DElasticData::_dtStableExplicit =   1.92450090e-01;
 
diff --git a/unittests/libtests/materials/data/PowerLaw3DTimeDep.py b/unittests/libtests/materials/data/PowerLaw3DTimeDep.py
index 56edff0..ab12c12 100644
--- a/unittests/libtests/materials/data/PowerLaw3DTimeDep.py
+++ b/unittests/libtests/materials/data/PowerLaw3DTimeDep.py
@@ -210,7 +210,7 @@ class PowerLaw3DTimeDep(ElasticMaterialApp):
     maxwellTimeB = self._getMaxwellTime(muB, refStrainRateB, refStressB, \
                                         powerLawExponentB, stressB)
 
-    self.dtStableImplicit = 0.1*min(maxwellTimeA, maxwellTimeB)
+    self.dtStableImplicit = 0.2 * min(maxwellTimeA, maxwellTimeB)
     self.dtStableExplicit = 1000.0 / vpA
 
     return
@@ -290,8 +290,8 @@ class PowerLaw3DTimeDep(ElasticMaterialApp):
     effStress = (0.5 * devStressProd)**0.5
     maxwellTime = 1.0
     if (effStress != 0.0):
-      maxwellTime = 0.5 * (refStress/effStress)**(powerLawExponent - 1.0) * \
-                    (refStress/mu)/refStrainRate
+      maxwellTime = (refStress/effStress)**(powerLawExponent - 1.0) * \
+                    (refStress/mu)/(refStrainRate * 6.0)
 
     return maxwellTime
 
diff --git a/unittests/libtests/materials/data/PowerLaw3DTimeDepData.cc b/unittests/libtests/materials/data/PowerLaw3DTimeDepData.cc
index ff80d7b..9c12034 100644
--- a/unittests/libtests/materials/data/PowerLaw3DTimeDepData.cc
+++ b/unittests/libtests/materials/data/PowerLaw3DTimeDepData.cc
@@ -45,7 +45,7 @@ const PylithScalar pylith::materials::PowerLaw3DTimeDepData::_pressureScale =
 
 const PylithScalar pylith::materials::PowerLaw3DTimeDepData::_densityScale =   2.25000000e+04;
 
-const PylithScalar pylith::materials::PowerLaw3DTimeDepData::_dtStableImplicit =   4.44444444e+06;
+const PylithScalar pylith::materials::PowerLaw3DTimeDepData::_dtStableImplicit =   2.96296296e+06;
 
 const PylithScalar pylith::materials::PowerLaw3DTimeDepData::_dtStableExplicit =   1.92450090e-01;
 
diff --git a/unittests/libtests/materials/data/PowerLawPlaneStrainElastic.py b/unittests/libtests/materials/data/PowerLawPlaneStrainElastic.py
index aed8764..e025424 100644
--- a/unittests/libtests/materials/data/PowerLawPlaneStrainElastic.py
+++ b/unittests/libtests/materials/data/PowerLawPlaneStrainElastic.py
@@ -213,7 +213,7 @@ class PowerLawPlaneStrainElastic(ElasticMaterialApp):
                                         stressUpdated[1,:])
 
 
-    self.dtStableImplicit = 0.1*min(maxwellTimeA, maxwellTimeB)
+    self.dtStableImplicit = 0.2 * min(maxwellTimeA, maxwellTimeB)
     self.dtStableExplicit = 1000.0 / vpA
 
     return
@@ -247,8 +247,8 @@ class PowerLawPlaneStrainElastic(ElasticMaterialApp):
     effStress = (0.5 * devStressProd)**0.5
     maxwellTime = 1.0
     if (effStress != 0.0):
-      maxwellTime = 0.5 * (refStress/effStress)**(powerLawExponent - 1.0) * \
-                    (refStress/mu)/refStrainRate
+      maxwellTime = (refStress/effStress)**(powerLawExponent - 1.0) * \
+                    (refStress/mu)/(refStrainRate * 6.0)
 
     return maxwellTime
 
diff --git a/unittests/libtests/materials/data/PowerLawPlaneStrainElasticData.cc b/unittests/libtests/materials/data/PowerLawPlaneStrainElasticData.cc
index 8954679..eccecfb 100644
--- a/unittests/libtests/materials/data/PowerLawPlaneStrainElasticData.cc
+++ b/unittests/libtests/materials/data/PowerLawPlaneStrainElasticData.cc
@@ -45,7 +45,7 @@ const PylithScalar pylith::materials::PowerLawPlaneStrainElasticData::_pressureS
 
 const PylithScalar pylith::materials::PowerLawPlaneStrainElasticData::_densityScale =   2.25000000e+04;
 
-const PylithScalar pylith::materials::PowerLawPlaneStrainElasticData::_dtStableImplicit =   4.09893495e+06;
+const PylithScalar pylith::materials::PowerLawPlaneStrainElasticData::_dtStableImplicit =   2.73262330e+06;
 
 const PylithScalar pylith::materials::PowerLawPlaneStrainElasticData::_dtStableExplicit =   1.92450090e-01;
 
diff --git a/unittests/libtests/materials/data/PowerLawPlaneStrainTimeDep.py b/unittests/libtests/materials/data/PowerLawPlaneStrainTimeDep.py
index be6dee1..d4e8011 100644
--- a/unittests/libtests/materials/data/PowerLawPlaneStrainTimeDep.py
+++ b/unittests/libtests/materials/data/PowerLawPlaneStrainTimeDep.py
@@ -222,7 +222,7 @@ class PowerLawPlaneStrainTimeDep(ElasticMaterialApp):
                                         powerLawExponentB,
                                         self.stateVarsUpdated[1,5:])
 
-    self.dtStableImplicit = 0.1*min(maxwellTimeA, maxwellTimeB)
+    self.dtStableImplicit = 0.2 * min(maxwellTimeA, maxwellTimeB)
     self.dtStableExplicit = 1000.0 / vpA
 
     return
@@ -300,8 +300,8 @@ class PowerLawPlaneStrainTimeDep(ElasticMaterialApp):
     effStress = (0.5 * devStressProd)**0.5
     maxwellTime = 1.0
     if (effStress != 0.0):
-      maxwellTime = 0.5 * (refStress/effStress)**(powerLawExponent - 1.0) * \
-                    (refStress/mu)/refStrainRate
+      maxwellTime = (refStress/effStress)**(powerLawExponent - 1.0) * \
+                    (refStress/mu)/(refStrainRate * 6.0)
 
     return maxwellTime
 
diff --git a/unittests/libtests/materials/data/PowerLawPlaneStrainTimeDepData.cc b/unittests/libtests/materials/data/PowerLawPlaneStrainTimeDepData.cc
index 4aa88ab..a9da4fb 100644
--- a/unittests/libtests/materials/data/PowerLawPlaneStrainTimeDepData.cc
+++ b/unittests/libtests/materials/data/PowerLawPlaneStrainTimeDepData.cc
@@ -45,7 +45,7 @@ const PylithScalar pylith::materials::PowerLawPlaneStrainTimeDepData::_pressureS
 
 const PylithScalar pylith::materials::PowerLawPlaneStrainTimeDepData::_densityScale =   2.25000000e+04;
 
-const PylithScalar pylith::materials::PowerLawPlaneStrainTimeDepData::_dtStableImplicit =   4.44444444e+06;
+const PylithScalar pylith::materials::PowerLawPlaneStrainTimeDepData::_dtStableImplicit =   2.96296296e+06;
 
 const PylithScalar pylith::materials::PowerLawPlaneStrainTimeDepData::_dtStableExplicit =   1.92450090e-01;
 



More information about the CIG-COMMITS mailing list