[cig-commits] r15557 - in short/3D/PyLith/trunk: doc/userguide doc/userguide/materials examples/3d/hex8 examples/twocells/twohex8 libsrc/feassemble libsrc/materials playpen/faultfaces/test pylith/materials unittests/libtests/materials unittests/libtests/materials/data
willic3 at geodynamics.org
willic3 at geodynamics.org
Tue Aug 18 16:06:30 PDT 2009
Author: willic3
Date: 2009-08-18 16:06:29 -0700 (Tue, 18 Aug 2009)
New Revision: 15557
Modified:
short/3D/PyLith/trunk/doc/userguide/materials/materials.lyx
short/3D/PyLith/trunk/doc/userguide/userguide.lyx
short/3D/PyLith/trunk/examples/3d/hex8/gravity_istress.cfg
short/3D/PyLith/trunk/examples/3d/hex8/pylithapp.cfg
short/3D/PyLith/trunk/examples/twocells/twohex8/pylithapp.cfg
short/3D/PyLith/trunk/examples/twocells/twohex8/twohex8.mesh
short/3D/PyLith/trunk/libsrc/feassemble/ElasticityImplicit.cc
short/3D/PyLith/trunk/libsrc/materials/PowerLaw3D.cc
short/3D/PyLith/trunk/libsrc/materials/PowerLaw3D.hh
short/3D/PyLith/trunk/playpen/faultfaces/test/pylithapp.cfg
short/3D/PyLith/trunk/pylith/materials/PowerLaw3D.py
short/3D/PyLith/trunk/unittests/libtests/materials/TestPowerLaw3D.cc
short/3D/PyLith/trunk/unittests/libtests/materials/data/PowerLaw3DElastic.py
short/3D/PyLith/trunk/unittests/libtests/materials/data/PowerLaw3DElasticData.cc
short/3D/PyLith/trunk/unittests/libtests/materials/data/PowerLaw3DTimeDep.py
short/3D/PyLith/trunk/unittests/libtests/materials/data/PowerLaw3DTimeDepData.cc
Log:
Updated documentation and code for revised power-law, including unit tests.
Modified: short/3D/PyLith/trunk/doc/userguide/materials/materials.lyx
===================================================================
--- short/3D/PyLith/trunk/doc/userguide/materials/materials.lyx 2009-08-18 23:00:18 UTC (rev 15556)
+++ short/3D/PyLith/trunk/doc/userguide/materials/materials.lyx 2009-08-18 23:06:29 UTC (rev 15557)
@@ -1,4 +1,4 @@
-#LyX 1.6.0 created this file. For more info see http://www.lyx.org/
+#LyX 1.6.3 created this file. For more info see http://www.lyx.org/
\lyxformat 345
\begin_document
\begin_header
@@ -391,6 +391,9 @@
settings) will cause two files to be created for each material group: an
info file, which describes the material properties used in the model, and
a state variables file, which contains the state variable information.
+ Note that the material properties described by the info file are the properties
+ used internally by PyLith.
+ They do not necessarily correspond to properties specified by the user.
If the problem has more than one time step, a state variable output file
will be created for each requested time step.
We have requested that the values be averaged over each cell.
@@ -443,7 +446,7 @@
\end_layout
\begin_layout Standard
-The properties and state variable available for output in each material
+The properties and state variables available for output in each material
model are listed in Table
\begin_inset CommandInset ref
LatexCommand ref
@@ -1674,18 +1677,23 @@
\end_inset
, the second deviatoric stress invariant,
-\begin_inset Formula $\sqrt{J_{2}^{\prime}}$
+\begin_inset Formula $J_{2}^{\prime}$
\end_inset
-, and the effective deviatoric strain,
+, the effective deviatoric strain,
\begin_inset Formula $\overline{e}$
\end_inset
+, and the second deviatoric strain invariant,
+\begin_inset Formula $L_{2}^{\prime}$
+\end_inset
+
, as
\begin_inset Formula \begin{gather}
-\overline{\sigma}=\sqrt{\frac{3}{2}\underline{S}\cdot\underline{S}}\,\,\label{eq:17}\\
-\sqrt{J_{2}^{\prime}}=\sqrt{\frac{1}{2}\underline{S}\cdot\underline{S}}\,.\nonumber \\
-\overline{e}=\sqrt{\frac{2}{3}\underline{e^{\prime}}\cdot\underline{e^{\prime}}}\,\,\nonumber \end{gather}
+\overline{\sigma}=\sqrt{\frac{3}{2}\underline{S}\cdot\underline{S}}\,\,\nonumber \\
+J_{2}^{\prime}=\frac{1}{2}\underline{S}\cdot\underline{S}\,.\label{eq:17}\\
+\overline{e}=\sqrt{\frac{2}{3}\underline{e}\cdot\underline{e}}\,\,\nonumber \\
+L_{2}^{\prime}=\frac{1}{2}\underline{e}\cdot\underline{e}\,\,\nonumber \end{gather}
\end_inset
@@ -2247,13 +2255,13 @@
into deviatoric and volumetric parts:
\begin_inset Formula \begin{gather}
-^{t+\Delta t}\underline{S}=\frac{E}{1+\nu}\left(^{t+\Delta t}\underline{e}^{\prime}-^{t+\Delta t}\underline{e}^{C}-\underline{e}^{I}\right)+\underline{S}^{I}\label{eq:42}\\
+^{t+\Delta t}\underline{S}=\frac{E}{1+\nu}\left(^{t+\Delta t}\underline{e}-^{t+\Delta t}\underline{e}^{C}-\underline{e}^{I}\right)+\underline{S}^{I}\label{eq:42}\\
^{t+\Delta t}P=\frac{E}{1-2\nu}\left(^{t+\Delta t}\theta-\theta^{I}\right)+P^{I}\:,\nonumber \end{gather}
\end_inset
where
-\begin_inset Formula $^{t+\Delta t}\underline{e}^{\prime}$
+\begin_inset Formula $^{t+\Delta t}\underline{e}$
\end_inset
is the total deviatoric strain,
@@ -2299,13 +2307,13 @@
may also be written as
\begin_inset Formula \begin{gather}
-^{t+\Delta t}\underline{S}=\frac{E}{1+\mathrm{\nu}}(^{t+\Delta t}\underline{e}^{\prime\prime}-\underline{\Delta e}^{C})+\underline{S}^{I}\,,\label{eq:43}\end{gather}
+^{t+\Delta t}\underline{S}=\frac{E}{1+\mathrm{\nu}}(^{t+\Delta t}\underline{e}^{\prime}-\underline{\Delta e}^{C})+\underline{S}^{I}\,,\label{eq:43}\end{gather}
\end_inset
where
\begin_inset Formula \begin{gather}
-^{t+\Delta t}\underline{e}^{\prime\prime}=^{t+\Delta t}\underline{e}^{\prime}-^{t}\underline{e}^{C}-\underline{e}^{I}\,\,,\,\,\,\underline{\Delta e}^{C}=^{t+\Delta t}\underline{e}^{C}-^{t}\underline{e}^{C}\,.\label{eq:44}\end{gather}
+^{t+\Delta t}\underline{e}^{\prime}=^{t+\Delta t}\underline{e}-^{t}\underline{e}^{C}-\underline{e}^{I}\,\,,\,\,\,\underline{\Delta e}^{C}=^{t+\Delta t}\underline{e}^{C}-^{t}\underline{e}^{C}\,.\label{eq:44}\end{gather}
\end_inset
@@ -2407,7 +2415,7 @@
, we obtain
\begin_inset Formula \begin{gather}
-^{t+\Delta t}\underline{S}=\frac{E}{1+\nu}\left\{ ^{t+\Delta t}\underline{e}^{\prime\prime}-\frac{\Delta t}{2\eta}\left[(1-\alpha)^{t}\underline{S}+\alpha^{t+\Delta t}\underline{S}\right]\right\} +\underline{S}^{I}\,\,.\label{eq:52}\end{gather}
+^{t+\Delta t}\underline{S}=\frac{E}{1+\nu}\left\{ ^{t+\Delta t}\underline{e}^{\prime}-\frac{\Delta t}{2\eta}\left[(1-\alpha)^{t}\underline{S}+\alpha^{t+\Delta t}\underline{S}\right]\right\} +\underline{S}^{I}\,\,.\label{eq:52}\end{gather}
\end_inset
@@ -2417,7 +2425,7 @@
,
\begin_inset Formula \begin{gather}
-^{t+\Delta t}\underline{S}=\frac{1}{\frac{1+\mathrm{\nu}}{E}+\frac{\alpha\Delta t}{2\eta}}\left[^{t+\Delta t}\underline{e}^{\prime\prime}-\frac{\Delta t}{2\eta}(1-\alpha)^{t}\underline{S}+\frac{1+\mathrm{\nu}}{E}\underline{S}^{I}\right]\,\,.\label{eq:53}\end{gather}
+^{t+\Delta t}\underline{S}=\frac{1}{\frac{1+\mathrm{\nu}}{E}+\frac{\alpha\Delta t}{2\eta}}\left[^{t+\Delta t}\underline{e}^{\prime}-\frac{\Delta t}{2\eta}(1-\alpha)^{t}\underline{S}+\frac{1+\mathrm{\nu}}{E}\underline{S}^{I}\right]\,\,.\label{eq:53}\end{gather}
\end_inset
@@ -2458,7 +2466,7 @@
Using the chain rule,
\begin_inset Formula \begin{gather}
-\frac{\partial^{t+\Delta t}S_{i}}{\partial^{t+\Delta t}\epsilon_{j}}=\frac{\partial^{t+\Delta t}S_{i}}{\partial^{t+\Delta t}e_{k}^{\prime\prime}}\frac{\partial^{t+\Delta t}e_{k}^{\prime\prime}}{\partial^{t+\Delta t}e_{l}^{\prime}}\frac{\partial^{t+\Delta t}e_{l}^{\prime}}{\partial^{t+\Delta t}\epsilon_{j}}\,\,.\label{eq:58}\end{gather}
+\frac{\partial^{t+\Delta t}S_{i}}{\partial^{t+\Delta t}\epsilon_{j}}=\frac{\partial^{t+\Delta t}S_{i}}{\partial^{t+\Delta t}e_{k}^{\prime}}\frac{\partial^{t+\Delta t}e_{k}^{\prime}}{\partial^{t+\Delta t}e_{l}}\frac{\partial^{t+\Delta t}e_{l}}{\partial^{t+\Delta t}\epsilon_{j}}\,\,.\label{eq:58}\end{gather}
\end_inset
@@ -2471,7 +2479,7 @@
, we obtain
\begin_inset Formula \begin{gather}
-\frac{\partial^{t+\Delta t}e_{k}^{\prime\prime}}{\partial^{t+\Delta t}e_{l}^{\prime}}=\delta_{kl}\,\,,\label{eq:59}\end{gather}
+\frac{\partial^{t+\Delta t}e_{k}^{\prime}}{\partial^{t+\Delta t}e_{l}}=\delta_{kl}\,\,,\label{eq:59}\end{gather}
\end_inset
@@ -2484,11 +2492,11 @@
:
\begin_inset Formula \begin{gather}
-\frac{\partial^{t+\Delta t}e_{l}^{\prime}}{\partial^{t+\Delta t}e_{j}}=\frac{1}{3}\left[\begin{array}{ccc}
+\frac{\partial^{t+\Delta t}e_{l}}{\partial^{t+\Delta t}\epsilon_{j}}=\frac{1}{3}\left[\begin{array}{ccc}
2 & -1 & -1\\
-1 & 2 & -1\\
-1 & -1 & 2\end{array}\right];\,\,1\leq l,j\leq3\label{eq:60}\\
-\frac{\partial^{t+\Delta t}e_{l}^{\prime}}{\partial^{t+\Delta t}\epsilon_{j}}=\delta_{lj}\,\,;\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\textrm{otherwise.}\nonumber \end{gather}
+\frac{\partial^{t+\Delta t}e_{l}}{\partial^{t+\Delta t}\epsilon_{j}}=\delta_{lj}\,\,;\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\textrm{otherwise.}\nonumber \end{gather}
\end_inset
@@ -2501,7 +2509,7 @@
, we have
\begin_inset Formula \begin{gather}
-\frac{\partial^{t+\Delta t}S_{i}}{\partial^{t+\Delta t}e_{k}^{\prime\prime}}=\frac{\delta_{ik}}{\frac{1+\nu}{E}+\frac{\alpha\Delta t}{2\eta}}\,\,.\label{eq:61}\end{gather}
+\frac{\partial^{t+\Delta t}S_{i}}{\partial^{t+\Delta t}e_{k}^{\prime}}=\frac{\delta_{ik}}{\frac{1+\nu}{E}+\frac{\alpha\Delta t}{2\eta}}\,\,.\label{eq:61}\end{gather}
\end_inset
@@ -2569,41 +2577,282 @@
\end_layout
\begin_layout Standard
-A power-law Maxwell viscoelastic material may be parameterized with the
- elastic parameters
-\begin_inset Formula $E$
+Laboratory results on rock rheology are typically performed using a triaxial
+ experiment, and the creep data are fit to a power-law equation of the form
+ (e.g.,
+\begin_inset CommandInset citation
+LatexCommand cite
+key "Kirby:Kronenberg:1987"
+
\end_inset
+):
+\begin_inset Formula \begin{equation}
+\dot{\epsilon}_{11}^{C}=A_{E}\exp\left(\frac{-Q}{RT}\right)\left(\sigma_{1}-\sigma_{3}\right)^{n}=A_{E}\exp\left(\frac{-Q}{RT}\right)\sigma_{d}^{n}\:,\label{eq:64}\end{equation}
+
+\end_inset
+
+where
+\begin_inset Formula $\dot{\epsilon}_{11}^{C}$
+\end_inset
+
+ is the strain rate in the direction of the maximum principal stress
+\begin_inset Formula $\left(\sigma_{1}\right)$
+\end_inset
+
+,
+\begin_inset Formula $A_{E}$
+\end_inset
+
+ is the experimentally-derived pre-exponential constant,
+\begin_inset Formula $Q$
+\end_inset
+
+ is the activation enthalpy,
+\begin_inset Formula $R$
+\end_inset
+
+ is the universal gas constant,
+\begin_inset Formula $T$
+\end_inset
+
+ is the absolute temperature,
+\begin_inset Formula $n$
+\end_inset
+
+ is the power-law exponent,
+\begin_inset Formula $\sigma_{3}\:\left(=\sigma_{2}\right)$
+\end_inset
+
+ is equal to the confining pressure, and
+\begin_inset Formula $\sigma_{d}$
+\end_inset
+
+ is the differential stress.
+ To properly formulate the flow law, it must be generalized so that the
+ results are not influenced by the experiment type or the choice of coordinate
+ systems (e.g.,
+\begin_inset CommandInset citation
+LatexCommand cite
+key "Paterson:1994"
+
+\end_inset
+
+).
+ The flow law may then be generalized in terms of the deviatoric stress
+ and strain rate invariants:
+\begin_inset Formula \begin{equation}
+\sqrt{\dot{L}_{2}^{\prime C}}=A_{M}\exp\left(\frac{-Q}{RT}\right)\sqrt{J_{2}^{\prime}}^{n}\:,\label{eq:65}\end{equation}
+
+\end_inset
+
+where
+\begin_inset Formula $A_{M}$
+\end_inset
+
+ is now a pre-exponential constant used in the formulation for modeling.
+ In practice, it is necessary to compute each strain rate component using
+ the flow law.
+ This is accomplished using:
+\begin_inset Formula \begin{equation}
+\dot{e}_{ij}^{C}=A_{M}\exp\left(\frac{-Q}{RT}\right)\sqrt{J_{2}^{\prime}}^{n-1}S_{ij}\:.\label{eq:66}\end{equation}
+
+\end_inset
+
+Note that equations
+\begin_inset CommandInset ref
+LatexCommand ref
+reference "eq:65"
+
+\end_inset
+
and
-\begin_inset Formula $\nu$
+\begin_inset CommandInset ref
+LatexCommand ref
+reference "eq:66"
+
\end_inset
-, the viscosity coefficient
-\begin_inset Formula $\eta$
+ are consistent, since equation
+\begin_inset CommandInset ref
+LatexCommand ref
+reference "eq:65"
+
\end_inset
-, and the power-law exponent
+ may be obtained from equation
+\begin_inset CommandInset ref
+LatexCommand ref
+reference "eq:66"
+
+\end_inset
+
+ by taking the scalar inner product of both sides, multiplying by 1/2, and
+ taking the square root.
+\end_layout
+
+\begin_layout Standard
+In a triaxial experiment with confining pressure
+\begin_inset Formula $P_{c}$
+\end_inset
+
+, we have
+\begin_inset Formula \begin{gather}
+\sigma_{2}=\sigma_{3}=P_{c}\nonumber \\
+\sigma_{1}=\sigma_{1}^{app}\label{eq:67}\\
+P=\frac{\sigma_{1}+2P_{c}}{3}\:,\nonumber \end{gather}
+
+\end_inset
+
+where
+\begin_inset Formula $\sigma_{1}^{app}$
+\end_inset
+
+ is the applied load.
+ The deviatoric stresses are then:
+\begin_inset Formula \begin{gather}
+S_{1}=\frac{2}{3}\left(\sigma_{1}-P_{c}\right)\nonumber \\
+S_{2}=S_{3}=-\frac{1}{3}\left(\sigma_{1}-P_{c}\right)\:.\label{eq:68}\end{gather}
+
+\end_inset
+
+This gives
+\begin_inset Formula \begin{gather}
+S_{1}=\frac{2}{3}\left(\sigma_{1}-\sigma_{3}\right)=\frac{2}{3}\sigma_{d}\nonumber \\
+S_{2}=S_{3}=-\frac{1}{3}\left(\sigma_{1}-\sigma_{3}\right)=-\frac{1}{3}\sigma_{d}\:.\label{eq:69}\end{gather}
+
+\end_inset
+
+In terms of the second deviatoric stress invariant, we then have
+\begin_inset Formula \begin{equation}
+\sqrt{J_{2}^{\prime}}=\frac{\sigma_{d}}{\sqrt{3}}\:.\label{eq:70}\end{equation}
+
+\end_inset
+
+
+\end_layout
+
+\begin_layout Standard
+Under the assumption that the creep measured in the laboratory experiments
+ is incompressible, we have
+\begin_inset Formula \begin{gather}
+\dot{e}_{11}^{C}=\dot{\epsilon}_{11}\nonumber \\
+\dot{e}_{22}^{C}=\dot{e}_{33}^{C}=-\frac{1}{2}\dot{\epsilon}_{11}\:.\label{eq:71}\end{gather}
+
+\end_inset
+
+In terms of the second deviatoric strain rate invariant we then have
+\begin_inset Formula \begin{equation}
+\sqrt{\dot{L}_{2}^{\prime C}}=\frac{\sqrt{3}}{2}\dot{\epsilon}_{11}\:.\label{eq:72}\end{equation}
+
+\end_inset
+
+Substituting equations
+\begin_inset CommandInset ref
+LatexCommand ref
+reference "eq:70"
+
+\end_inset
+
+ and
+\begin_inset CommandInset ref
+LatexCommand ref
+reference "eq:72"
+
+\end_inset
+
+ into equation
+\begin_inset CommandInset ref
+LatexCommand ref
+reference "eq:64"
+
+\end_inset
+
+, we obtain
+\begin_inset Formula \begin{equation}
+\sqrt{\dot{L}_{2}^{\prime C}}=A_{E}\frac{\sqrt{3}^{n+1}}{2}\exp\left(\frac{-Q}{RT}\right)\sqrt{J_{2}^{\prime}}^{n}\:,\label{eq:73}\end{equation}
+
+\end_inset
+
+and therefore,
+\begin_inset Formula \begin{equation}
+A_{M}=\frac{\sqrt{3}^{n+1}}{2}A_{E}\:.\label{eq:74}\end{equation}
+
+\end_inset
+
+
+\end_layout
+
+\begin_layout Standard
+Since users will generally want to make direct use of laboratory results,
+ PyLith converts the coefficients internally, and users are asked to specify
+ the power-law exponent,
\begin_inset Formula $n$
\end_inset
-.
- The parameters defining this model are shown in Table
+, as well as a power-law coefficient:
+\begin_inset Formula \begin{equation}
+A_{T}=A_{E}\exp\left(\frac{-Q}{RT}\right)\:.\label{eq:75}\end{equation}
+
+\end_inset
+
+In Table
\begin_inset CommandInset ref
LatexCommand ref
reference "tab:powerLaw"
\end_inset
+, the properties
+\family typewriter
+power-law-coefficient
+\family default
+ and
+\family typewriter
+power-law-exponent
+\family default
+ refer to
+\begin_inset Formula $A_{T}$
+\end_inset
+
+ and
+\begin_inset Formula $n$
+\end_inset
+
+, respectively.
+ To avoid scaling issues, PyLith internally makes use of a viscosity coefficient
+,
+\begin_inset Formula $\eta$
+\end_inset
+
+, rather than the power-law coefficient:
+\begin_inset Formula \begin{gather}
+\dot{e}_{ij}^{C}=\frac{\sqrt{J_{2}^{\prime}}^{n-1}S_{ij}}{2\eta^{n}}\label{eq:76}\\
+\eta=\frac{\exp\left(\frac{Q}{nRT}\right)}{\left(2A_{M}\right)^{\frac{1}{n}}}\:.\nonumber \end{gather}
+
+\end_inset
+
+This viscosity coefficient is the value written to info files, as listed
+ in Table
+\begin_inset CommandInset ref
+LatexCommand ref
+reference "tab:material-model-output"
+
+\end_inset
+
.
- The creep strain increment is approximated as
+\end_layout
+
+\begin_layout Standard
+The creep strain increment is approximated as
\begin_inset Formula \begin{gather}
-\underline{\Delta e}^{C}=\frac{\Delta t\sqrt{^{\tau}J_{2}^{\prime}}^{n-1}\,^{\tau}\underline{S}}{2\eta^{n}}=\frac{\Delta t^{\tau}\overline{\sigma}^{n-1}\,^{\tau}\underline{S}}{2\sqrt{3}\eta^{n}}\,.\label{eq:64}\end{gather}
+\underline{\Delta e}^{C}=\frac{\Delta t\sqrt{^{\tau}J_{2}^{\prime}}^{n-1}\,^{\tau}\underline{S}}{2\eta^{n}}=\frac{\Delta t^{\tau}\overline{\sigma}^{n-1}\,^{\tau}\underline{S}}{2\sqrt{3}\eta^{n}}\,.\label{eq:77}\end{gather}
\end_inset
Therefore,
\begin_inset Formula \begin{gather}
-\Delta\bar{e}^{C}=\frac{\Delta t\sqrt{^{\tau}J_{2}^{\prime}}^{n}}{\sqrt{3}\eta^{n}}=\frac{\Delta t^{\tau}\overline{\sigma}^{n}}{\sqrt{3}^{n+1}\eta^{n}}\,,\,\textrm{and}\,^{\tau}\gamma=\frac{\sqrt{^{\tau}J_{2}^{\prime}}^{n-1}}{2\eta^{n}}\,.\label{eq:65}\end{gather}
+\Delta\bar{e}^{C}=\frac{\Delta t\sqrt{^{\tau}J_{2}^{\prime}}^{n}}{\sqrt{3}\eta^{n}}=\frac{\Delta t^{\tau}\overline{\sigma}^{n}}{\sqrt{3}^{n+1}\eta^{n}}\,,\,\textrm{and}\,^{\tau}\gamma=\frac{\sqrt{^{\tau}J_{2}^{\prime}}^{n-1}}{2\eta^{n}}\,.\label{eq:78}\end{gather}
\end_inset
@@ -2617,14 +2866,14 @@
,
\begin_inset CommandInset ref
LatexCommand ref
-reference "eq:64"
+reference "eq:77"
\end_inset
, and
\begin_inset CommandInset ref
LatexCommand ref
-reference "eq:65"
+reference "eq:78"
\end_inset
@@ -2637,27 +2886,27 @@
, we obtain:
\begin_inset Formula \begin{gather}
-^{t+\Delta t}\underline{S}=\frac{E}{1+\nu}\left\{ ^{t+\Delta t}\underline{e}^{\prime\prime}-\Delta t^{\tau}\gamma\left[\left(1-\alpha\right)^{t}\underline{S}+\alpha{}^{t+\Delta t}\underline{S}\right]\right\} +\underline{S}^{I}\,,\label{eq:66}\end{gather}
+^{t+\Delta t}\underline{S}=\frac{E}{1+\nu}\left\{ ^{t+\Delta t}\underline{e}^{\prime}-\Delta t^{\tau}\gamma\left[\left(1-\alpha\right)^{t}\underline{S}+\alpha{}^{t+\Delta t}\underline{S}\right]\right\} +\underline{S}^{I}\,,\label{eq:79}\end{gather}
\end_inset
which may be rewritten:
\begin_inset Formula \begin{gather}
-^{t+\Delta t}\underline{S}\left(\frac{1+\nu}{E}+\alpha\Delta t^{\tau}\gamma\right)={}^{t+\Delta t}\underline{e}^{\prime\prime}-\Delta t^{\tau}\gamma\left(1-\alpha\right)^{t}\underline{S}+\frac{1+\nu}{E}\underline{S}^{I}\,.\label{eq:67}\end{gather}
+^{t+\Delta t}\underline{S}\left(\frac{1+\nu}{E}+\alpha\Delta t^{\tau}\gamma\right)={}^{t+\Delta t}\underline{e}^{\prime}-\Delta t^{\tau}\gamma\left(1-\alpha\right)^{t}\underline{S}+\frac{1+\nu}{E}\underline{S}^{I}\,.\label{eq:80}\end{gather}
\end_inset
Taking the scalar inner product of both sides we obtain:
\begin_inset Formula \begin{gather}
-a^{2}\,\,{}^{t+\Delta t}J_{2}^{\prime}-b+c^{\tau}\gamma-d^{2}\,^{\tau}\gamma^{2}=F=0\,,\label{eq:68}\end{gather}
+a^{2}\,\,{}^{t+\Delta t}J_{2}^{\prime}-b+c^{\tau}\gamma-d^{2}\,^{\tau}\gamma^{2}=F=0\,,\label{eq:81}\end{gather}
\end_inset
where
\begin_inset Formula \begin{gather}
a=\frac{1+\nu}{E}+\alpha\Delta t^{\tau}\gamma\,\,\nonumber \\
-b=\frac{1}{2}{}^{t+\Delta t}\underline{e}^{\prime\prime}\cdot{}^{t+\Delta t}\underline{e}^{\prime\prime}+\frac{1+\nu}{E}{}^{t+\Delta t}\underline{e}^{\prime\prime}\cdot\underline{S}^{I}+\left(\frac{1+\nu}{E}\right)^{2}\,^{I}J_{2}^{\prime}\,.\label{eq:69}\\
-c=\Delta t\left(1-\alpha\right){}^{t+\Delta t}\underline{e}^{\prime\prime}\cdot^{t}\underline{S}+\Delta t\left(1-\alpha\right)\frac{1+\nu}{E}\,^{t}\underline{S}\cdot\underline{S}^{I}\,\,\nonumber \\
+b=\frac{1}{2}{}^{t+\Delta t}\underline{e}^{\prime}\cdot{}^{t+\Delta t}\underline{e}^{\prime}+\frac{1+\nu}{E}{}^{t+\Delta t}\underline{e}^{\prime}\cdot\underline{S}^{I}+\left(\frac{1+\nu}{E}\right)^{2}\,^{I}J_{2}^{\prime}\,.\label{eq:82}\\
+c=\Delta t\left(1-\alpha\right){}^{t+\Delta t}\underline{e}^{\prime}\cdot^{t}\underline{S}+\Delta t\left(1-\alpha\right)\frac{1+\nu}{E}\,^{t}\underline{S}\cdot\underline{S}^{I}\,\,\nonumber \\
d=\Delta t\left(1-\alpha\right)\sqrt{^{t}J_{2}^{\prime}}\,\,\nonumber \end{gather}
\end_inset
@@ -2665,7 +2914,7 @@
Equation
\begin_inset CommandInset ref
LatexCommand ref
-reference "eq:68"
+reference "eq:81"
\end_inset
@@ -2679,14 +2928,14 @@
time step may be found from Equations
\begin_inset CommandInset ref
LatexCommand ref
-reference "eq:64"
+reference "eq:78"
\end_inset
and
\begin_inset CommandInset ref
LatexCommand ref
-reference "eq:65"
+reference "eq:79"
\end_inset
@@ -2713,19 +2962,19 @@
We begin by rewriting Equation
\begin_inset CommandInset ref
LatexCommand ref
-reference "eq:67"
+reference "eq:80"
\end_inset
as
\begin_inset Formula \begin{gather}
-F=^{t+\Delta t}S_{i}\left(a_{E}+\alpha\Delta t^{\tau}\gamma\right)-^{t+\Delta t}e_{i}^{\prime\prime}+\Delta t^{\tau}\gamma\left(1-\alpha\right)^{t}S_{i}-a_{E}S_{i}^{I}=0\label{eq:70}\end{gather}
+F=^{t+\Delta t}S_{i}\left(a_{E}+\alpha\Delta t^{\tau}\gamma\right)-^{t+\Delta t}e_{i}^{\prime}+\Delta t^{\tau}\gamma\left(1-\alpha\right)^{t}S_{i}-a_{E}S_{i}^{I}=0\label{eq:83}\end{gather}
\end_inset
where
\begin_inset Formula \begin{gather}
-a_{E}=\frac{1+\nu}{E}\,.\label{eq:71}\end{gather}
+a_{E}=\frac{1+\nu}{E}\,.\label{eq:84}\end{gather}
\end_inset
@@ -2735,7 +2984,7 @@
is
\begin_inset Formula \begin{gather}
-\frac{\partial F}{\partial^{t+\Delta t}e_{k}^{\prime\prime}}=-\delta_{ik}\:,\label{eq:72}\end{gather}
+\frac{\partial F}{\partial^{t+\Delta t}e_{k}^{\prime}}=-\delta_{ik}\:,\label{eq:85}\end{gather}
\end_inset
@@ -2745,14 +2994,14 @@
is
\begin_inset Formula \begin{gather}
-\frac{\partial F}{\partial^{t+\Delta t}S_{i}}=a_{E}+\alpha\Delta t^{\tau}\gamma+\frac{\partial^{\tau}\gamma}{\partial^{t+\Delta t}S_{i}}\Delta t\left[\alpha^{t+\Delta t}S_{i}+\left(1-\alpha\right)^{t}S_{i}\right]\:.\label{eq:73}\end{gather}
+\frac{\partial F}{\partial^{t+\Delta t}S_{i}}=a_{E}+\alpha\Delta t^{\tau}\gamma+\frac{\partial^{\tau}\gamma}{\partial^{t+\Delta t}S_{i}}\Delta t\left[\alpha^{t+\Delta t}S_{i}+\left(1-\alpha\right)^{t}S_{i}\right]\:.\label{eq:86}\end{gather}
\end_inset
From Equation
\begin_inset CommandInset ref
LatexCommand ref
-reference "eq:65"
+reference "eq:78"
\end_inset
@@ -2765,20 +3014,20 @@
,
\begin_inset Formula \begin{gather}
-^{\tau}\gamma=\frac{1}{2\eta^{n}}\left[\alpha\sqrt{^{t+\Delta t}J_{2}^{\prime}}+\left(1-\alpha\right)\sqrt{^{t}J_{2}^{\prime}}\right]^{n-1}\:.\label{eq:74}\end{gather}
+^{\tau}\gamma=\frac{1}{2\eta^{n}}\left[\alpha\sqrt{^{t+\Delta t}J_{2}^{\prime}}+\left(1-\alpha\right)\sqrt{^{t}J_{2}^{\prime}}\right]^{n-1}\:.\label{eq:87}\end{gather}
\end_inset
Then
\begin_inset Formula \begin{gather}
-\frac{\partial^{\tau}\gamma}{\partial{}^{t+\Delta t}S_{i}}=\frac{\partial^{\tau}\gamma}{\partial\sqrt{^{t+\Delta t}J_{2}^{\prime}}}\frac{\partial\sqrt{^{t+\Delta t}J_{2}^{\prime}}}{\partial^{t+\Delta t}S_{l}}\label{eq:75}\\
+\frac{\partial^{\tau}\gamma}{\partial{}^{t+\Delta t}S_{i}}=\frac{\partial^{\tau}\gamma}{\partial\sqrt{^{t+\Delta t}J_{2}^{\prime}}}\frac{\partial\sqrt{^{t+\Delta t}J_{2}^{\prime}}}{\partial^{t+\Delta t}S_{l}}\label{eq:88}\\
=\frac{\alpha\left(n-1\right)\sqrt{^{\tau}J_{2}^{\prime}}^{n-2}{}^{t+\Delta t}T_{i}}{4\eta^{n}}\,,\nonumber \end{gather}
\end_inset
Where
\begin_inset Formula \begin{gather}
-^{t+\Delta t}T_{i}=^{t+\Delta t}S_{i}\:;\:\:1\leq i\leq3\label{eq:76}\\
+^{t+\Delta t}T_{i}=^{t+\Delta t}S_{i}\:;\:\:1\leq i\leq3\label{eq:89}\\
^{t+\Delta t}T_{i}=2^{t+\Delta t}S_{i}\:;\:\:\textrm{otherwise.}\nonumber \end{gather}
\end_inset
@@ -2786,27 +3035,27 @@
Then using Equations
\begin_inset CommandInset ref
LatexCommand ref
-reference "eq:72"
+reference "eq:85"
\end_inset
,
\begin_inset CommandInset ref
LatexCommand ref
-reference "eq:73"
+reference "eq:86"
\end_inset
,
\begin_inset CommandInset ref
LatexCommand ref
-reference "eq:75"
+reference "eq:88"
\end_inset
, and the quotient rule for derivatives of an implicit function,
\begin_inset Formula \begin{gather}
-\frac{\partial^{t+\Delta t}S_{i}}{\partial{}^{t+\Delta t}e_{k}^{\prime\prime}}=\frac{\delta_{ik}}{a_{E}+\alpha\Delta t\left[^{\tau}\gamma+\frac{^{\tau}S_{i}\left(n-1\right){}^{t+\Delta t}T_{i}\sqrt{^{\tau}J_{2}^{\prime}}^{n-2}}{4\sqrt{^{t+\Delta t}J_{2}^{\prime}}\eta^{n}}\right]}\,.\label{eq:77}\end{gather}
+\frac{\partial^{t+\Delta t}S_{i}}{\partial{}^{t+\Delta t}e_{k}^{\prime}}=\frac{\delta_{ik}}{a_{E}+\alpha\Delta t\left[^{\tau}\gamma+\frac{^{\tau}S_{i}\left(n-1\right){}^{t+\Delta t}T_{i}\sqrt{^{\tau}J_{2}^{\prime}}^{n-2}}{4\sqrt{^{t+\Delta t}J_{2}^{\prime}}\eta^{n}}\right]}\,.\label{eq:90}\end{gather}
\end_inset
@@ -2844,13 +3093,13 @@
1 & 1 & 1 & 0 & 0 & 0\\
0 & 0 & 0 & 0 & 0 & 0\\
0 & 0 & 0 & 0 & 0 & 0\\
-0 & 0 & 0 & 0 & 0 & 0\end{array}\right]+\frac{1}{3}\frac{\partial{}^{t+\Delta t}S_{i}}{\partial{}^{t+\Delta t}e_{k}^{\prime\prime}}\left[\begin{array}{cccccc}
+0 & 0 & 0 & 0 & 0 & 0\end{array}\right]+\frac{1}{3}\frac{\partial{}^{t+\Delta t}S_{i}}{\partial{}^{t+\Delta t}e_{k}^{\prime}}\left[\begin{array}{cccccc}
2 & -1 & -1 & 0 & 0 & 0\\
-1 & 2 & -1 & 0 & 0 & 0\\
-1 & -1 & 2 & 0 & 0 & 0\\
0 & 0 & 0 & 3 & 0 & 0\\
0 & 0 & 0 & 0 & 3 & 0\\
-0 & 0 & 0 & 0 & 0 & 3\end{array}\right]\,.\label{eq:78}\end{gather}
+0 & 0 & 0 & 0 & 0 & 3\end{array}\right]\,.\label{eq:91}\end{gather}
\end_inset
@@ -2862,14 +3111,14 @@
goes to infinity), Equations
\begin_inset CommandInset ref
LatexCommand ref
-reference "eq:77"
+reference "eq:90"
\end_inset
and
\begin_inset CommandInset ref
LatexCommand ref
-reference "eq:78"
+reference "eq:91"
\end_inset
@@ -2881,13 +3130,13 @@
we require the derivative of equation
\begin_inset CommandInset ref
LatexCommand ref
-reference "eq:68"
+reference "eq:81"
\end_inset
, which may be written:
\begin_inset Formula \begin{gather}
-\frac{\partial F}{\partial\sqrt{^{t+\Delta t}J_{2}^{\prime}}}=2a^{2}\sqrt{^{t+\Delta t}J_{2}^{\prime}}+\frac{\alpha\left(n-1\right)\sqrt{^{\tau}J_{2}^{\prime}}^{n-2}}{2\eta^{n}}\left(2a\alpha\Delta t{}^{t+\Delta t}J_{2}^{\prime}+c-2d^{2}\,^{\tau}\gamma\right)\,.\label{eq:79}\end{gather}
+\frac{\partial F}{\partial\sqrt{^{t+\Delta t}J_{2}^{\prime}}}=2a^{2}\sqrt{^{t+\Delta t}J_{2}^{\prime}}+\frac{\alpha\left(n-1\right)\sqrt{^{\tau}J_{2}^{\prime}}^{n-2}}{2\eta^{n}}\left(2a\alpha\Delta t{}^{t+\Delta t}J_{2}^{\prime}+c-2d^{2}\,^{\tau}\gamma\right)\,.\label{eq:92}\end{gather}
\end_inset
@@ -3503,7 +3752,7 @@
\begin_inset Text
\begin_layout Plain Layout
-\begin_inset Formula $\eta$
+\begin_inset Formula $A_{T}$
\end_inset
@@ -3517,7 +3766,7 @@
\begin_layout Plain Layout
\family typewriter
-viscosity-coefficient
+power-law-coefficient
\end_layout
\end_inset
Modified: short/3D/PyLith/trunk/doc/userguide/userguide.lyx
===================================================================
--- short/3D/PyLith/trunk/doc/userguide/userguide.lyx 2009-08-18 23:00:18 UTC (rev 15556)
+++ short/3D/PyLith/trunk/doc/userguide/userguide.lyx 2009-08-18 23:06:29 UTC (rev 15557)
@@ -1,4 +1,4 @@
-#LyX 1.6.2 created this file. For more info see http://www.lyx.org/
+#LyX 1.6.3 created this file. For more info see http://www.lyx.org/
\lyxformat 345
\begin_document
\begin_header
@@ -583,5 +583,42 @@
, 3369-3376.
\end_layout
+\begin_layout Bibliography
+\begin_inset CommandInset bibitem
+LatexCommand bibitem
+label "12"
+key "Kirby:Kronenberg:1987"
+
+\end_inset
+
+Kirby, S.
+ H.
+ and A.
+ K.
+ Kronenberg (1987), Rheology of the lithosphere: Selected topics,
+\shape italic
+Reviews of Geophysics, 25
+\shape default
+, 1219-1244.
+\end_layout
+
+\begin_layout Bibliography
+\begin_inset CommandInset bibitem
+LatexCommand bibitem
+label "13"
+key "Paterson:1994"
+
+\end_inset
+
+Paterson, W.
+ S.
+ B.
+ (1994),
+\shape italic
+The Physics of Glaciers, Third Edition
+\shape default
+, Elsevier Science Ltd., Oxford, 480 pp.
+\end_layout
+
\end_body
\end_document
Modified: short/3D/PyLith/trunk/examples/3d/hex8/gravity_istress.cfg
===================================================================
--- short/3D/PyLith/trunk/examples/3d/hex8/gravity_istress.cfg 2009-08-18 23:00:18 UTC (rev 15556)
+++ short/3D/PyLith/trunk/examples/3d/hex8/gravity_istress.cfg 2009-08-18 23:06:29 UTC (rev 15557)
@@ -101,6 +101,3 @@
cell_info_fields = [density,mu,lambda,maxwell_time]
cell_filter = pylith.meshio.CellFilterAvgMesh
writer.filename = gravity_istress-viscoelastic.vtk
-
-[pylithapp.petsc]
-ksp_atol = 1.0e-12
Modified: short/3D/PyLith/trunk/examples/3d/hex8/pylithapp.cfg
===================================================================
--- short/3D/PyLith/trunk/examples/3d/hex8/pylithapp.cfg 2009-08-18 23:00:18 UTC (rev 15556)
+++ short/3D/PyLith/trunk/examples/3d/hex8/pylithapp.cfg 2009-08-18 23:06:29 UTC (rev 15557)
@@ -6,16 +6,16 @@
# ----------------------------------------------------------------------
# Turn on some journals to show progress.
[pylithapp.journal.info]
-##timedependent = 1
-##implicit = 1
-##petsc = 1
-##solverlinear = 1
-##meshiocubit = 1
-##implicitelasticity = 1
-##faultcohesivekin = 1
-##fiatlagrange = 1
+timedependent = 1
+implicit = 1
+petsc = 1
+solverlinear = 1
+meshiocubit = 1
+implicitelasticity = 1
+faultcohesivekin = 1
+fiatlagrange = 1
quadrature3d = 1
-##pylithapp = 1
+pylithapp = 1
##problem = 1
##implicit = 1
@@ -84,7 +84,7 @@
log_summary = true
ksp_max_it = 100
ksp_gmres_restart = 50
-#start_in_debugger = true
+# start_in_debugger = true
snes_monitor = true
snes_view = true
Modified: short/3D/PyLith/trunk/examples/twocells/twohex8/pylithapp.cfg
===================================================================
--- short/3D/PyLith/trunk/examples/twocells/twohex8/pylithapp.cfg 2009-08-18 23:00:18 UTC (rev 15556)
+++ short/3D/PyLith/trunk/examples/twocells/twohex8/pylithapp.cfg 2009-08-18 23:06:29 UTC (rev 15557)
@@ -81,7 +81,7 @@
# using a spatial database file.
db_properties = spatialdata.spatialdb.UniformDB
db_properties.values = [vp,vs,density,viscosity]
-db_properties.data = [5773.502691896258*m/s,3333.333333333333*m/s,2700.0*kg/m**3,1.0e18*Pa*s]
+db_properties.data = [5291.502622129181*m/s,3000.0*m/s,2700.0*kg/m**3,7.10046e19*Pa*s]
# If we instead wanted to used the 'matprops.spatialdb' file to define
# material properties we would comment out the three 'db' lines above
@@ -99,7 +99,24 @@
# arguments may be found in the PETSc documentation.
# ----------------------------------------------------------------------
[pylithapp.petsc]
-pc_type = jacobi
+ksp_rtol = 1.0e-8
+pc_type = asm
-# start_in_debugger = true
-# debugger_timeout = 100
+# Change the preconditioner settings (must turn off
+# shift_positive_definite and turn on shift_nonzero).
+sub_pc_factor_shift_positive_definite = 0
+sub_pc_factor_shift_nonzero =
+
+ksp_monitor = true
+ksp_view = true
+log_summary = true
+ksp_max_it = 100
+ksp_gmres_restart = 50
+
+
+start_in_debugger = true
+
+snes_monitor = true
+snes_view = true
+ksp_converged_reason = true
+snes_converged_reason = true
Modified: short/3D/PyLith/trunk/examples/twocells/twohex8/twohex8.mesh
===================================================================
--- short/3D/PyLith/trunk/examples/twocells/twohex8/twohex8.mesh 2009-08-18 23:00:18 UTC (rev 15556)
+++ short/3D/PyLith/trunk/examples/twocells/twohex8/twohex8.mesh 2009-08-18 23:06:29 UTC (rev 15557)
@@ -114,6 +114,17 @@
}
// This group of vertices may be used to specify boundary conditions.
+ // There are 6 vertices corresponding to indices 1, 3, 5, 7, 9, and 11.
+ group = {
+ name = y_pos
+ type = vertices
+ count = 6
+ indices = {
+ 1 3 5 7 9 11
+ }
+ }
+
+ // This group of vertices may be used to specify boundary conditions.
// There are 6 vertices corresponding to indices 0, 1, 2, 3, 4, and 5.
group = {
name = z_neg
Modified: short/3D/PyLith/trunk/libsrc/feassemble/ElasticityImplicit.cc
===================================================================
--- short/3D/PyLith/trunk/libsrc/feassemble/ElasticityImplicit.cc 2009-08-18 23:00:18 UTC (rev 15556)
+++ short/3D/PyLith/trunk/libsrc/feassemble/ElasticityImplicit.cc 2009-08-18 23:06:29 UTC (rev 15557)
@@ -287,6 +287,7 @@
_logger->eventEnd(computeEvent);
} // if
+ // residualSection->view("After gravity contribution");
// Compute B(transpose) * sigma, first computing strains
_logger->eventBegin(stressEvent);
calcTotalStrainFn(&strainCell, basisDeriv, dispTpdtCell,
@@ -308,6 +309,7 @@
_logger->eventBegin(updateEvent);
residualVisitor.clear();
sieveMesh->updateAdd(*c_iter, residualVisitor);
+ // residualSection->view("After stress contribution");
_logger->eventEnd(updateEvent);
} // for
} // integrateResidual
Modified: short/3D/PyLith/trunk/libsrc/materials/PowerLaw3D.cc
===================================================================
--- short/3D/PyLith/trunk/libsrc/materials/PowerLaw3D.cc 2009-08-18 23:00:18 UTC (rev 15556)
+++ short/3D/PyLith/trunk/libsrc/materials/PowerLaw3D.cc 2009-08-18 23:06:29 UTC (rev 15557)
@@ -51,14 +51,14 @@
{ "density", 1, pylith::topology::FieldBase::SCALAR },
{ "mu", 1, pylith::topology::FieldBase::SCALAR },
{ "lambda", 1, pylith::topology::FieldBase::SCALAR },
- { "viscosity_coeff", 1, pylith::topology::FieldBase::SCALAR },
+ { "viscosity_coefficient", 1, pylith::topology::FieldBase::SCALAR },
{ "power_law_exponent", 1, pylith::topology::FieldBase::SCALAR }
};
// Values expected in properties spatial database
const int numDBProperties = 5;
const char* dbProperties[] = {"density", "vs", "vp" ,
- "viscosity-coeff",
+ "power-law-coefficient",
"power-law-exponent"};
/// Number of state variables.
@@ -99,11 +99,11 @@
const int pylith::materials::PowerLaw3D::p_lambda =
pylith::materials::PowerLaw3D::p_mu + 1;
-const int pylith::materials::PowerLaw3D::p_viscosityCoeff =
+const int pylith::materials::PowerLaw3D::p_viscosityCoefficient =
pylith::materials::PowerLaw3D::p_lambda + 1;
const int pylith::materials::PowerLaw3D::p_powerLawExponent =
- pylith::materials::PowerLaw3D::p_viscosityCoeff + 1;
+ pylith::materials::PowerLaw3D::p_viscosityCoefficient + 1;
// Indices of property database values (order must match dbProperties).
const int pylith::materials::PowerLaw3D::db_density = 0;
@@ -114,11 +114,11 @@
const int pylith::materials::PowerLaw3D::db_vp =
pylith::materials::PowerLaw3D::db_vs + 1;
-const int pylith::materials::PowerLaw3D::db_viscosityCoeff =
+const int pylith::materials::PowerLaw3D::db_powerLawCoefficient =
pylith::materials::PowerLaw3D::db_vp + 1;
const int pylith::materials::PowerLaw3D::db_powerLawExponent =
- pylith::materials::PowerLaw3D::db_viscosityCoeff + 1;
+ pylith::materials::PowerLaw3D::db_powerLawCoefficient + 1;
// Indices of state variables.
const int pylith::materials::PowerLaw3D::s_viscousStrain = 0;
@@ -198,10 +198,10 @@
const double density = dbValues[db_density];
const double vs = dbValues[db_vs];
const double vp = dbValues[db_vp];
- const double viscosityCoeff = dbValues[db_viscosityCoeff];
+ const double powerLawCoefficient = dbValues[db_powerLawCoefficient];
const double powerLawExponent = dbValues[db_powerLawExponent];
- if (density <= 0.0 || vs <= 0.0 || vp <= 0.0 || viscosityCoeff <= 0.0
+ if (density <= 0.0 || vs <= 0.0 || vp <= 0.0 || powerLawCoefficient <= 0.0
|| powerLawExponent < 1.0) {
std::ostringstream msg;
msg << "Spatial database returned illegal value for physical "
@@ -209,7 +209,7 @@
<< "density: " << density << "\n"
<< "vp: " << vp << "\n"
<< "vs: " << vs << "\n"
- << "viscosityCoeff: " << viscosityCoeff << "\n"
+ << "powerLawCoefficient: " << powerLawCoefficient << "\n"
<< "powerLawExponent: " << powerLawExponent << "\n";
throw std::runtime_error(msg.str());
} // if
@@ -227,10 +227,14 @@
} // if
assert(mu > 0);
+ const double factor = powerLawCoefficient *
+ pow(sqrt(3.0), (powerLawExponent + 1.0));
+ const double viscosityCoefficient = pow(factor, -1.0/powerLawExponent);
+
propValues[p_density] = density;
propValues[p_mu] = mu;
propValues[p_lambda] = lambda;
- propValues[p_viscosityCoeff] = viscosityCoeff;
+ propValues[p_viscosityCoefficient] = viscosityCoefficient;
propValues[p_powerLawExponent] = powerLawExponent;
PetscLogFlops(6);
@@ -259,8 +263,8 @@
_normalizer->nondimensionalize(values[p_mu], pressureScale);
values[p_lambda] =
_normalizer->nondimensionalize(values[p_lambda], pressureScale);
- values[p_viscosityCoeff] =
- _normalizer->nondimensionalize(values[p_viscosityCoeff],
+ values[p_viscosityCoefficient] =
+ _normalizer->nondimensionalize(values[p_viscosityCoefficient],
viscosityCoeffScale);
PetscLogFlops(7);
@@ -289,8 +293,8 @@
_normalizer->dimensionalize(values[p_mu], pressureScale);
values[p_lambda] =
_normalizer->dimensionalize(values[p_lambda], pressureScale);
- values[p_viscosityCoeff] =
- _normalizer->dimensionalize(values[p_viscosityCoeff],
+ values[p_viscosityCoefficient] =
+ _normalizer->dimensionalize(values[p_viscosityCoefficient],
viscosityCoeffScale);
PetscLogFlops(7);
@@ -380,7 +384,7 @@
assert(0 != stateVars);
assert(_numVarsQuadPt == numStateVars);
const double mu = properties[p_mu];
- const double viscosityCoeff = properties[p_viscosityCoeff];
+ const double viscosityCoefficient = properties[p_viscosityCoefficient];
const double powerLawExp = properties[p_powerLawExponent];
const double stress[] = {stateVars[s_stress],
@@ -401,8 +405,8 @@
double dtTest = 1.0;
if (effStress != 0.0)
dtTest = 0.1 *
- pow((viscosityCoeff/effStress), (powerLawExp - 1.0)) *
- (viscosityCoeff/mu);
+ pow((viscosityCoefficient/effStress), (powerLawExp - 1.0)) *
+ (viscosityCoefficient/mu);
const double dtStable = dtTest;
#if 0 // DEBUGGING
@@ -509,7 +513,7 @@
const double mu = properties[p_mu];
const double lambda = properties[p_lambda];
- const double viscosityCoeff = properties[p_viscosityCoeff];
+ const double viscosityCoefficient = properties[p_viscosityCoefficient];
const double powerLawExp = properties[p_powerLawExponent];
const double visStrainT[] = {stateVars[s_viscousStrain],
stateVars[s_viscousStrain + 1],
@@ -614,7 +618,7 @@
_effStressParams.dt = _dt;
_effStressParams.effStressT = effStressT;
_effStressParams.powerLawExp = powerLawExp;
- _effStressParams.viscosityCoeff = viscosityCoeff;
+ _effStressParams.viscosityCoefficient = viscosityCoefficient;
const double effStressInitialGuess = effStressT;
@@ -627,7 +631,8 @@
const double effStressTau = (1.0 - alpha) * effStressT +
alpha * effStressTpdt;
const double gammaTau = 0.5 *
- pow((effStressTau/viscosityCoeff), (powerLawExp - 1.0))/viscosityCoeff;
+ pow((effStressTau/viscosityCoefficient),
+ (powerLawExp - 1.0))/viscosityCoefficient;
const double factor1 = 1.0/(ae + alpha * _dt * gammaTau);
const double factor2 = timeFac * gammaTau;
double devStressTpdt = 0.0;
@@ -663,11 +668,11 @@
const double dt = _effStressParams.dt;
const double effStressT = _effStressParams.effStressT;
const double powerLawExp = _effStressParams.powerLawExp;
- const double viscosityCoeff = _effStressParams.viscosityCoeff;
+ const double viscosityCoefficient = _effStressParams.viscosityCoefficient;
const double factor1 = 1.0-alpha;
const double effStressTau = factor1 * effStressT + alpha * effStressTpdt;
- const double gammaTau = 0.5 * pow((effStressTau/viscosityCoeff),
- (powerLawExp - 1.0))/viscosityCoeff;
+ const double gammaTau = 0.5 * pow((effStressTau/viscosityCoefficient),
+ (powerLawExp - 1.0))/viscosityCoefficient;
const double a = ae + alpha * dt * gammaTau;
const double y = a * a * effStressTpdt * effStressTpdt - b +
c * gammaTau - d * d * gammaTau * gammaTau;
@@ -690,15 +695,15 @@
const double dt = _effStressParams.dt;
const double effStressT = _effStressParams.effStressT;
const double powerLawExp = _effStressParams.powerLawExp;
- const double viscosityCoeff = _effStressParams.viscosityCoeff;
+ const double viscosityCoefficient = _effStressParams.viscosityCoefficient;
const double factor1 = 1.0-alpha;
const double effStressTau = factor1 * effStressT + alpha * effStressTpdt;
- const double gammaTau = 0.5 * pow((effStressTau/viscosityCoeff),
- (powerLawExp - 1.0))/viscosityCoeff;
+ const double gammaTau = 0.5 * pow((effStressTau/viscosityCoefficient),
+ (powerLawExp - 1.0))/viscosityCoefficient;
const double a = ae + alpha * dt * gammaTau;
const double dGammaTau = 0.5 * alpha * (powerLawExp - 1.0) *
- pow((effStressTau/viscosityCoeff), (powerLawExp - 2.0))/
- (viscosityCoeff * viscosityCoeff);
+ pow((effStressTau/viscosityCoefficient), (powerLawExp - 2.0))/
+ (viscosityCoefficient * viscosityCoefficient);
const double dy = 2.0 * a * a * effStressTpdt + dGammaTau *
(2.0 * a * alpha * dt * effStressTpdt * effStressTpdt +
c - 2.0 * d * d * gammaTau);
@@ -726,14 +731,14 @@
const double dt = _effStressParams.dt;
const double effStressT = _effStressParams.effStressT;
const double powerLawExp = _effStressParams.powerLawExp;
- const double viscosityCoeff = _effStressParams.viscosityCoeff;
+ const double viscosityCoefficient = _effStressParams.viscosityCoefficient;
const double factor1 = 1.0-alpha;
const double effStressTau = factor1 * effStressT + alpha * effStressTpdt;
- const double gammaTau = 0.5 * pow((effStressTau/viscosityCoeff),
- (powerLawExp - 1.0))/viscosityCoeff;
+ const double gammaTau = 0.5 * pow((effStressTau/viscosityCoefficient),
+ (powerLawExp - 1.0))/viscosityCoefficient;
const double dGammaTau = 0.5 * alpha * (powerLawExp - 1.0) *
- pow((effStressTau/viscosityCoeff), (powerLawExp - 2.0))/
- (viscosityCoeff * viscosityCoeff);
+ pow((effStressTau/viscosityCoefficient), (powerLawExp - 2.0))/
+ (viscosityCoefficient * viscosityCoefficient);
const double a = ae + alpha * dt * gammaTau;
y = a * a * effStressTpdt * effStressTpdt -
b +
@@ -846,7 +851,7 @@
const int tensorSize = _tensorSize;
const double mu = properties[p_mu];
const double lambda = properties[p_lambda];
- const double viscosityCoeff = properties[p_viscosityCoeff];
+ const double viscosityCoefficient = properties[p_viscosityCoefficient];
const double powerLawExp = properties[p_powerLawExponent];
const double stress[] = {stateVars[s_stress],
stateVars[s_stress + 1],
@@ -868,12 +873,13 @@
stress[5]};
const double effStress = sqrt(0.5 * _scalarProduct(devStress, devStress));
const double gamma = 0.5 *
- pow((effStress/viscosityCoeff), (powerLawExp - 1.0))/viscosityCoeff;
+ pow((effStress/viscosityCoefficient),
+ (powerLawExp - 1.0))/viscosityCoefficient;
// const double visFac = 1.0/(3.0 * (ae + _dt * gamma));
const double visFac1 = 3.0 *(ae + _dt * gamma);
const double visFac2 = 0.75 * _dt * (powerLawExp - 1.0) *
- pow((effStress/viscosityCoeff), (powerLawExp - 2.0))/
- (viscosityCoeff * viscosityCoeff * effStress);
+ pow((effStress/viscosityCoefficient), (powerLawExp - 2.0))/
+ (viscosityCoefficient * viscosityCoefficient * effStress);
const double dStress11dStrain11 = 1.0/
(visFac1 + devStress[0] * devStress[0] * visFac2);
const double dStress22dStrain22 = 1.0/
@@ -948,7 +954,7 @@
const double mu = properties[p_mu];
const double lambda = properties[p_lambda];
- const double viscosityCoeff = properties[p_viscosityCoeff];
+ const double viscosityCoefficient = properties[p_viscosityCoefficient];
const double powerLawExp = properties[p_powerLawExponent];
// State variables.
@@ -1068,7 +1074,7 @@
_effStressParams.dt = _dt;
_effStressParams.effStressT = effStressT;
_effStressParams.powerLawExp = powerLawExp;
- _effStressParams.viscosityCoeff = viscosityCoeff;
+ _effStressParams.viscosityCoefficient = viscosityCoefficient;
const double effStressInitialGuess = effStressT;
@@ -1081,7 +1087,8 @@
const double effStressTau = (1.0 - alpha) * effStressT +
alpha * effStressTpdt;
const double gammaTau = 0.5 *
- pow((effStressTau/viscosityCoeff), (powerLawExp - 1.0))/viscosityCoeff;
+ pow((effStressTau/viscosityCoefficient),
+ (powerLawExp - 1.0))/viscosityCoefficient;
const double a = ae + alpha * _dt * gammaTau;
const double factor1 = 1.0/a;
const double factor2 = timeFac * gammaTau;
@@ -1106,8 +1113,8 @@
alpha * devStressT[4] + explicitFac * devStressTpdt[4],
alpha * devStressT[5] + explicitFac * devStressTpdt[5]};
const double factor3 = 0.25 * _dt * alpha * (powerLawExp - 1.0) *
- pow((effStressTau/viscosityCoeff), (powerLawExp - 2.0))/
- (viscosityCoeff * viscosityCoeff * effStressTpdt);
+ pow((effStressTau/viscosityCoefficient), (powerLawExp - 2.0))/
+ (viscosityCoefficient * viscosityCoefficient * effStressTpdt);
// Compute deviatoric derivatives
const double dStress11dStrain11 = 1.0/
@@ -1228,7 +1235,7 @@
// since otherwise we would have to redo a lot of calculations.
const double mu = properties[p_mu];
const double lambda = properties[p_lambda];
- const double viscosityCoeff = properties[p_viscosityCoeff];
+ const double viscosityCoefficient = properties[p_viscosityCoefficient];
const double powerLawExp = properties[p_powerLawExponent];
const double visStrainT[] = {stateVars[s_viscousStrain],
@@ -1330,7 +1337,7 @@
_effStressParams.dt = _dt;
_effStressParams.effStressT = effStressT;
_effStressParams.powerLawExp = powerLawExp;
- _effStressParams.viscosityCoeff = viscosityCoeff;
+ _effStressParams.viscosityCoefficient = viscosityCoefficient;
const double effStressInitialGuess = effStressT;
@@ -1344,7 +1351,8 @@
const double effStressTau = (1.0 - alpha) * effStressT +
alpha * effStressTpdt;
const double gammaTau = 0.5 *
- pow((effStressTau/viscosityCoeff), (powerLawExp - 1.0))/viscosityCoeff;
+ pow((effStressTau/viscosityCoefficient),
+ (powerLawExp - 1.0))/viscosityCoefficient;
const double factor1 = 1.0/(ae + alpha * _dt * gammaTau);
const double factor2 = timeFac * gammaTau;
double devStressTpdt = 0.0;
Modified: short/3D/PyLith/trunk/libsrc/materials/PowerLaw3D.hh
===================================================================
--- short/3D/PyLith/trunk/libsrc/materials/PowerLaw3D.hh 2009-08-18 23:00:18 UTC (rev 15556)
+++ short/3D/PyLith/trunk/libsrc/materials/PowerLaw3D.hh 2009-08-18 23:06:29 UTC (rev 15557)
@@ -16,11 +16,11 @@
*
* 3-D, isotropic, power-law Maxwell viscoelastic material. The
* physical properties are specified using density, shear-wave speed,
- * viscosity coefficient, power-law exponent, and compressional-wave speed.
+ * power-law coefficient, power-law exponent, and compressional-wave speed.
* The physical properties are stored internally using density, lambda, mu,
* which are directly related to the elasticity constants used in the
* finite-element integration. The viscosity information is retained as
- * specified.
+ * a viscosity coefficient and the power-law exponent.
*/
#if !defined(pylith_materials_powerlaw3d_hh)
@@ -510,7 +510,7 @@
double dt;
double effStressT;
double powerLawExp;
- double viscosityCoeff;
+ double viscosityCoefficient;
};
// PRIVATE MEMBERS ////////////////////////////////////////////////////
@@ -531,12 +531,12 @@
static const int p_density;
static const int p_mu;
static const int p_lambda;
- static const int p_viscosityCoeff;
+ static const int p_viscosityCoefficient;
static const int p_powerLawExponent;
static const int db_density;
static const int db_vs;
static const int db_vp;
- static const int db_viscosityCoeff;
+ static const int db_powerLawCoefficient;
static const int db_powerLawExponent;
static const int s_viscousStrain;
Modified: short/3D/PyLith/trunk/playpen/faultfaces/test/pylithapp.cfg
===================================================================
--- short/3D/PyLith/trunk/playpen/faultfaces/test/pylithapp.cfg 2009-08-18 23:00:18 UTC (rev 15556)
+++ short/3D/PyLith/trunk/playpen/faultfaces/test/pylithapp.cfg 2009-08-18 23:06:29 UTC (rev 15557)
@@ -1,5 +1,6 @@
# -*- Python -*-
[pylithapp]
+nodes = 4
# ----------------------------------------------------------------------
# journal
@@ -51,6 +52,8 @@
[pylithapp.timedependent.implicit]
output = [domain]
+# Uncomment to use MUMPS
+# matrix_type = aij
time_step.total_time = 0.0*s
time_step.dt = 1.0*s
@@ -130,8 +133,8 @@
filename = tets-f003.vtk
#Uncomment the following lines to include second fault.
-#[pylithapp.timedependent.interfaces.f005.output.writer]
-#filename = tets-f005.vtk
+[pylithapp.timedependent.interfaces.f005.output.writer]
+filename = tets-f005.vtk
[pylithapp.timedependent.materials.elastic_left.output]
cell_filter = pylith.meshio.CellFilterAvgMesh
Modified: short/3D/PyLith/trunk/pylith/materials/PowerLaw3D.py
===================================================================
--- short/3D/PyLith/trunk/pylith/materials/PowerLaw3D.py 2009-08-18 23:00:18 UTC (rev 15556)
+++ short/3D/PyLith/trunk/pylith/materials/PowerLaw3D.py 2009-08-18 23:06:29 UTC (rev 15557)
@@ -41,7 +41,7 @@
'data': []},
'cell': \
{'info': ["mu", "lambda", "density",
- "viscosity_coeff", "power_law_exponent"],
+ "viscosity_coefficient", "power_law_exponent"],
'data': ["total_strain", "stress", "viscous_strain"]}}
self._loggingPrefix = "MaPL3D "
return
Modified: short/3D/PyLith/trunk/unittests/libtests/materials/TestPowerLaw3D.cc
===================================================================
--- short/3D/PyLith/trunk/unittests/libtests/materials/TestPowerLaw3D.cc 2009-08-18 23:00:18 UTC (rev 15556)
+++ short/3D/PyLith/trunk/unittests/libtests/materials/TestPowerLaw3D.cc 2009-08-18 23:06:29 UTC (rev 15557)
@@ -195,7 +195,7 @@
CPPUNIT_ASSERT(material.hasProperty("mu"));
CPPUNIT_ASSERT(material.hasProperty("lambda"));
CPPUNIT_ASSERT(material.hasProperty("density"));
- CPPUNIT_ASSERT(material.hasProperty("viscosity_coeff"));
+ CPPUNIT_ASSERT(material.hasProperty("viscosity_coefficient"));
CPPUNIT_ASSERT(material.hasProperty("power_law_exponent"));
CPPUNIT_ASSERT(!material.hasProperty("aaa"));
} // testHasProperty
Modified: short/3D/PyLith/trunk/unittests/libtests/materials/data/PowerLaw3DElastic.py
===================================================================
--- short/3D/PyLith/trunk/unittests/libtests/materials/data/PowerLaw3DElastic.py 2009-08-18 23:00:18 UTC (rev 15556)
+++ short/3D/PyLith/trunk/unittests/libtests/materials/data/PowerLaw3DElastic.py 2009-08-18 23:06:29 UTC (rev 15557)
@@ -45,7 +45,7 @@
self.numLocs = numLocs
self.dbPropertyValues = ["density", "vs", "vp",
- "viscosity_coeff", "power_law_exponent"]
+ "power_law_coefficient", "power_law_exponent"]
self.numPropertyValues = numpy.array([1, 1, 1, 1, 1], dtype=numpy.int32)
self.dbStateVarValues = ["viscous-strain-xx",
@@ -66,24 +66,28 @@
densityA = 2500.0
vsA = 3000.0
vpA = vsA*3**0.5
- viscosityCoeffA = 1.0e18
+ powerLawCoeffA = 1.0/3.0e18
powerLawExponentA = 1.0
strainA = [1.1e-4, 1.2e-4, 1.3e-4, 1.4e-4, 1.5e-4, 1.6e-4]
initialStressA = [2.1e4, 2.2e4, 2.3e4, 2.4e4, 2.5e4, 2.6e4]
initialStrainA = [3.1e-4, 3.2e-4, 3.3e-4, 3.4e-4, 3.5e-4, 3.6e-4]
muA = vsA*vsA*densityA
lambdaA = vpA*vpA*densityA - 2.0*muA
+ viscosityCoeffA = (1.0/((3.0**0.5)**(powerLawExponentA + 1.0) \
+ * powerLawCoeffA))**(1.0/powerLawExponentA)
densityB = 2000.0
vsB = 1200.0
vpB = vsB*3**0.5
- viscosityCoeffB = 1.0e10
+ powerLawCoeffB = 1.0/9.0e30
powerLawExponentB = 3.0
strainB = [4.1e-4, 4.2e-4, 4.3e-4, 4.4e-4, 4.5e-4, 4.6e-4]
initialStressB = [5.1e4, 5.2e4, 5.3e4, 5.4e4, 5.5e4, 5.6e4]
initialStrainB = [6.1e-4, 6.2e-4, 6.3e-4, 6.4e-4, 6.5e-4, 6.6e-4]
muB = vsB*vsB*densityB
lambdaB = vpB*vpB*densityB - 2.0*muB
+ viscosityCoeffB = (1.0/((3.0**0.5)**(powerLawExponentB + 1.0) \
+ * powerLawCoeffB))**(1.0/powerLawExponentB)
self.lengthScale = 1.0e+3
self.pressureScale = muA
@@ -95,9 +99,9 @@
(self.timeScale**(1.0/powerLawExponentB)) * self.pressureScale
self.dbProperties = numpy.array([ [densityA, vsA, vpA, \
- viscosityCoeffA, powerLawExponentA],
+ powerLawCoeffA, powerLawExponentA],
[densityB, vsB, vpB, \
- viscosityCoeffB, powerLawExponentB] ],
+ powerLawCoeffB, powerLawExponentB] ],
dtype=numpy.float64)
self.properties = numpy.array([ [densityA, muA, lambdaA, \
viscosityCoeffA, powerLawExponentA],
Modified: short/3D/PyLith/trunk/unittests/libtests/materials/data/PowerLaw3DElasticData.cc
===================================================================
--- short/3D/PyLith/trunk/unittests/libtests/materials/data/PowerLaw3DElasticData.cc 2009-08-18 23:00:18 UTC (rev 15556)
+++ short/3D/PyLith/trunk/unittests/libtests/materials/data/PowerLaw3DElasticData.cc 2009-08-18 23:06:29 UTC (rev 15557)
@@ -58,7 +58,7 @@
"density",
"vs",
"vp",
-"viscosity_coeff",
+"power_law_coefficient",
"power_law_exponent",
};
@@ -81,12 +81,12 @@
2.50000000e+03,
3.00000000e+03,
5.19615242e+03,
- 1.00000000e+18,
+ 3.33333333e-19,
1.00000000e+00,
2.00000000e+03,
1.20000000e+03,
2.07846097e+03,
- 1.00000000e+10,
+ 1.11111111e-31,
3.00000000e+00,
};
Modified: short/3D/PyLith/trunk/unittests/libtests/materials/data/PowerLaw3DTimeDep.py
===================================================================
--- short/3D/PyLith/trunk/unittests/libtests/materials/data/PowerLaw3DTimeDep.py 2009-08-18 23:00:18 UTC (rev 15556)
+++ short/3D/PyLith/trunk/unittests/libtests/materials/data/PowerLaw3DTimeDep.py 2009-08-18 23:06:29 UTC (rev 15557)
@@ -48,9 +48,9 @@
self.numLocs = numLocs
self.dbPropertyValues = ["density", "vs", "vp",
- "viscosity_coeff", "power_law_exponent"]
+ "power_law_coefficient", "power_law_exponent"]
self.propertyValues = ["density", "mu", "lambda",
- "viscosity_coeff", "power_law_exponent"]
+ "viscosity_coefficient", "power_law_exponent"]
self.numPropertyValues = numpy.array([1, 1, 1, 1, 1], dtype=numpy.int32)
self.dbStateVarValues = ["viscous-strain-xx",
@@ -75,47 +75,48 @@
densityA = 2500.0
vsA = 3000.0
vpA = vsA*3**0.5
- viscosityCoeffA = 1.0e18
- powerLawExpA = 1.0
+ powerLawCoeffA = 1.0/3.0e18
+ powerLawExponentA = 1.0
strainA = [1.1e-4, 1.2e-4, 1.3e-4, 1.4e-4, 1.5e-4, 1.6e-4]
initialStressA = [2.1e4, 2.2e4, 2.3e4, 2.4e4, 2.5e4, 2.6e4]
initialStrainA = [3.6e-5, 3.5e-5, 3.4e-5, 3.3e-5, 3.2e-5, 3.1e-5]
muA = vsA*vsA*densityA
lambdaA = vpA*vpA*densityA - 2.0*muA
+ viscosityCoeffA = (1.0/((3.0**0.5)**(powerLawExponentA + 1.0) \
+ * powerLawCoeffA))**(1.0/powerLawExponentA)
+
densityB = 2000.0
vsB = 1200.0
vpB = vsB*3**0.5
- viscosityCoeffB = 1.0e12
- powerLawExpB = 3.0
+ powerLawCoeffB = 1.0/9.0e36
+ powerLawExponentB = 3.0
strainB = [4.1e-4, 4.2e-4, 4.3e-4, 4.4e-4, 4.5e-4, 4.6e-4]
initialStressB = [5.1e4, 5.2e4, 5.3e4, 5.4e4, 5.5e4, 5.6e4]
initialStrainB = [6.1e-5, 6.2e-5, 6.3e-5, 6.6e-5, 6.5e-5, 6.4e-5]
muB = vsB*vsB*densityB
lambdaB = vpB*vpB*densityB - 2.0*muB
+ viscosityCoeffB = (1.0/((3.0**0.5)**(powerLawExponentB + 1.0) \
+ * powerLawCoeffB))**(1.0/powerLawExponentB)
- # TEMPORARY, need to determine how to use initial strain
- # initialStrainA = [0.0, 0.0, 0.0, 0.0, 0.0, 0.0]
- # initialStrainB = [0.0, 0.0, 0.0, 0.0, 0.0, 0.0]
-
self.lengthScale = 1.0e+3
self.pressureScale = muA
self.timeScale = 1.0
self.densityScale = 1.0e+3
self.viscosityCoeffScaleA = \
- (self.timeScale**(1.0/powerLawExpA)) * self.pressureScale
+ (self.timeScale**(1.0/powerLawExponentA)) * self.pressureScale
self.viscosityCoeffScaleB = \
- (self.timeScale**(1.0/powerLawExpB)) * self.pressureScale
+ (self.timeScale**(1.0/powerLawExponentB)) * self.pressureScale
self.dbProperties = numpy.array([ [densityA, vsA, vpA, \
- viscosityCoeffA, powerLawExpA],
+ powerLawCoeffA, powerLawExponentA],
[densityB, vsB, vpB, \
- viscosityCoeffB, powerLawExpB] ],
+ powerLawCoeffB, powerLawExponentB] ],
dtype=numpy.float64)
self.properties = numpy.array([ [densityA, muA, lambdaA, \
- viscosityCoeffA, powerLawExpA],
+ viscosityCoeffA, powerLawExponentA],
[densityB, muB, lambdaB, \
- viscosityCoeffB, powerLawExpB] ],
+ viscosityCoeffB, powerLawExponentB] ],
dtype=numpy.float64)
# TEMPORARY, need to determine how to use initial state variables
@@ -131,11 +132,11 @@
numpy.array([ [densityA/density0, muA/mu0, \
lambdaA/mu0, \
viscosityCoeffA/viscosityCoeff0, \
- powerLawExpA], \
+ powerLawExponentA], \
[densityB/density0, muB/mu0, \
lambdaB/mu0, \
viscosityCoeffB/viscosityCoeff1, \
- powerLawExpB] ], \
+ powerLawExponentB] ], \
dtype=numpy.float64)
self.initialStress = numpy.array([initialStressA,
@@ -175,29 +176,29 @@
(self.elasticConsts[0,:], self.stress[0,:], self.stateVarsUpdated[0,:]) = \
self._calcStress(strainA,
muA, lambdaA, viscosityCoeffA,
- powerLawExpA,
+ powerLawExponentA,
visStrainA, stressA,
initialStressA, initialStrainA)
(self.elasticConsts[1,:], self.stress[1,:], self.stateVarsUpdated[1,:]) = \
self._calcStress(strainB,
muB, lambdaB, viscosityCoeffB,
- powerLawExpB,
+ powerLawExponentB,
visStrainB, stressB,
initialStressB, initialStrainB)
# Use state variables to compute Maxwell times (and stable time step size).
maxwellTimeA = self._getMaxwellTime(muA, viscosityCoeffA, \
- powerLawExpA, stressA)
+ powerLawExponentA, stressA)
maxwellTimeB = self._getMaxwellTime(muB, viscosityCoeffB, \
- powerLawExpB, stressB)
+ powerLawExponentB, stressB)
self.dtStableImplicit = 0.1*min(maxwellTimeA, maxwellTimeB)
return
def _bracket(self, effStressInitialGuess, ae, b, c, d, alpha, dt, effStressT,
- powerLawExpV, viscosityCoeffV):
+ powerLawExponentV, viscosityCoeffV):
"""
Function to bracket the effective stress.
"""
@@ -215,9 +216,9 @@
x2 = 1500.0
funcValue1 = self._effStressFunc(x1, ae, b, c, d, alpha, dt,
- effStressT, powerLawExpV, viscosityCoeffV)
+ effStressT, powerLawExponentV, viscosityCoeffV)
funcValue2 = self._effStressFunc(x2, ae, b, c, d, alpha, dt,
- effStressT, powerLawExpV, viscosityCoeffV)
+ effStressT, powerLawExponentV, viscosityCoeffV)
iteration = 0
bracketed = False
@@ -231,14 +232,14 @@
x1 = max(x1, 0.0)
funcValue1 = self._effStressFunc(x1, ae, b, c, d, alpha, dt,
- effStressT, powerLawExpV,
+ effStressT, powerLawExponentV,
viscosityCoeffV)
else:
x2 += bracketFactor * (x1 - x2)
x2 = max(x2, 0.0)
funcValue2 = self._effStressFunc(x2, ae, b, c, d, alpha, dt,
- effStressT, powerLawExpV,
+ effStressT, powerLawExponentV,
viscosityCoeffV)
iteration += 1
@@ -248,7 +249,7 @@
return x1, x2
- def _getMaxwellTime(self, mu, viscosityCoeff, powerLawExp, stress):
+ def _getMaxwellTime(self, mu, viscosityCoeff, powerLawExponent, stress):
"""
Compute Maxwell time from stress, viscosity coefficient, shear modulus, and
power-law exponent.
@@ -267,7 +268,7 @@
effStress = (0.5 * devStressProd)**0.5
maxwellTime = 1.0
if (effStress != 0.0):
- maxwellTime = (viscosityCoeff/effStress)**(powerLawExp - 1.0) * \
+ maxwellTime = (viscosityCoeff/effStress)**(powerLawExponent - 1.0) * \
(viscosityCoeff/mu)
return maxwellTime
@@ -288,7 +289,7 @@
def _calcStressComponent(self, strainVal, strainComp, stressComp, strainTpdt,
muV, lambdaV, viscosityCoeffV,
- powerLawExpV, visStrainT, stressT,
+ powerLawExponentV, visStrainT, stressT,
initialStress, initialStrain):
"""
Function to compute a particular stress component as a function of a
@@ -298,7 +299,8 @@
strainTest[strainComp] = strainVal
stressTpdt, visStrainTpdt = self._computeStress(strainTest, muV, lambdaV,
viscosityCoeffV,
- powerLawExpV, visStrainT,
+ powerLawExponentV,
+ visStrainT,
stressT,
initialStress,
initialStrain)
@@ -306,15 +308,15 @@
def _effStressFunc(self, effStressTpdt, ae, b, c, d, alpha, dt, effStressT,
- powerLawExpV, viscosityCoeffV):
+ powerLawExponentV, viscosityCoeffV):
"""
Function to compute effective stress function for a given effective stress.
"""
factor1 = 1.0 - alpha
effStressTau = factor1 * effStressT + alpha * effStressTpdt
- gammaTau = 0.5 * (effStressTau/viscosityCoeffV)**(powerLawExpV - 1.0)/ \
- viscosityCoeffV
+ gammaTau = 0.5 * (effStressTau/viscosityCoeffV)** \
+ (powerLawExponentV - 1.0) / viscosityCoeffV
a = ae + alpha * dt * gammaTau
effStressFunc = a * a * effStressTpdt * effStressTpdt - b + \
c * gammaTau - d * d * gammaTau * gammaTau
@@ -323,7 +325,7 @@
def _computeStress(self, strainTpdt, muV, lambdaV, viscosityCoeffV,
- powerLawExpV, visStrainT, stressT,
+ powerLawExponentV, visStrainT, stressT,
initialStress, initialStrain):
"""
Function to compute stresses and viscous strains using the
@@ -376,21 +378,22 @@
effStressInitialGuess = effStressT
x1, x2 = self._bracket(effStressInitialGuess, ae, b, c, d, self.alpha,
- self.dt, effStressT, powerLawExpV, viscosityCoeffV)
+ self.dt, effStressT, powerLawExponentV,
+ viscosityCoeffV)
# Find the root using Brent's method (from scipy)
rootTolerance = 1.0e-14
effStressTpdt = scipy.optimize.brentq(self._effStressFunc, x1, x2,
args=(ae, b, c, d, self.alpha,
self.dt, effStressT,
- powerLawExpV,
+ powerLawExponentV,
viscosityCoeffV),
xtol=rootTolerance)
# Compute stresses from the effective stress.
effStressTau = (1.0 - self.alpha) * effStressT + self.alpha * effStressTpdt
- gammaTau = 0.5 * ((effStressTau/viscosityCoeffV)**(powerLawExpV - 1.0)) / \
- viscosityCoeffV
+ gammaTau = 0.5 * ((effStressTau/viscosityCoeffV)** \
+ (powerLawExponentV - 1.0)) / viscosityCoeffV
factor1 = 1.0/(ae + self.alpha * self.dt * gammaTau)
factor2 = timeFac * gammaTau
devStressTpdt = 0.0
@@ -412,7 +415,7 @@
def _calcStress(self, strainV, muV, lambdaV, viscosityCoeffV,
- powerLawExpV,visStrainV, stressV,
+ powerLawExponentV,visStrainV, stressV,
initialStressV, initialStrainV):
"""
Compute stress, updated state variables and derivative of elasticity matrix.
@@ -429,7 +432,7 @@
stressTpdt, visStrainTpdt = self._computeStress(strainTpdt, muV, lambdaV,
viscosityCoeffV,
- powerLawExpV,
+ powerLawExponentV,
visStrainT, stressT,
initialStress,
initialStrain)
@@ -452,7 +455,8 @@
stressComp,
strainTpdt, muV, lambdaV,
viscosityCoeffV,
- powerLawExpV, visStrainT,
+ powerLawExponentV,
+ visStrainT,
stressT, initialStress,
initialStrain),
order=derivOrder)
Modified: short/3D/PyLith/trunk/unittests/libtests/materials/data/PowerLaw3DTimeDepData.cc
===================================================================
--- short/3D/PyLith/trunk/unittests/libtests/materials/data/PowerLaw3DTimeDepData.cc 2009-08-18 23:00:18 UTC (rev 15556)
+++ short/3D/PyLith/trunk/unittests/libtests/materials/data/PowerLaw3DTimeDepData.cc 2009-08-18 23:06:29 UTC (rev 15557)
@@ -58,7 +58,7 @@
"density",
"vs",
"vp",
-"viscosity_coeff",
+"power_law_coefficient",
"power_law_exponent",
};
@@ -81,12 +81,12 @@
2.50000000e+03,
3.00000000e+03,
5.19615242e+03,
- 1.00000000e+18,
+ 3.33333333e-19,
1.00000000e+00,
2.00000000e+03,
1.20000000e+03,
2.07846097e+03,
- 1.00000000e+12,
+ 1.11111111e-37,
3.00000000e+00,
};
More information about the CIG-COMMITS
mailing list