# [cig-commits] r13139 - in doc/snac: . figures

echoi at geodynamics.org echoi at geodynamics.org
Tue Oct 28 10:38:45 PDT 2008

Author: echoi
Date: 2008-10-28 10:38:45 -0700 (Tue, 28 Oct 2008)
New Revision: 13139

doc/snac/figures/snac-mesh.pdf
Modified:
doc/snac/snac.lyx
Log:
Updated the section 1.1 SNAC Implementation.

===================================================================
(Binary files differ)

Property changes on: doc/snac/figures/snac-mesh.pdf
___________________________________________________________________
Name: svn:mime-type
+ application/octet-stream

Modified: doc/snac/snac.lyx
===================================================================
--- doc/snac/snac.lyx	2008-10-28 16:14:29 UTC (rev 13138)
+++ doc/snac/snac.lyx	2008-10-28 17:38:45 UTC (rev 13139)
@@ -1,4 +1,4 @@
\lyxformat 276
\begin_document
@@ -259,13 +259,17 @@

\begin_layout Itemize
E.
- Choi, Thoutireddy, P., Lavier, L., Quenette, S., Tan, E., Gurnis, M., Aivazis
- M., and Appelbe, B., GeoFramework Part II: Coupling models of crustal deformation
- and mantle convection with a computational framework, to be submitted to
-
+ Choi, Lavier, L., Gurnis, M., Thermomechanics of mid-ocean ridge segmentation,
+ accepted by Physics of the Earth and Planetary Interiors, 2008,
\emph on
-Geochem., Geophys., Geosys.

+\begin_inset LatexCommand htmlurl
+name "http://dx.doi.org/10.1016/j.pepi.2008.08.010"
+target "http://dx.doi.org/10.1016/j.pepi.2008.08.010"
+
+\end_inset
+
+.
\end_layout

\begin_layout Standard
@@ -319,7 +323,7 @@

\begin_layout Standard
SNAC (StGermaiN Analysis of Continua) is an updated Lagrangian explicit
- finite element code for modeling a finitely deforming viscoplastic solid
+ finite difference code for modeling a finitely deforming viscoplastic solid
in 3D, released under the GNU General Public License (see Appendix
\begin_inset LatexCommand vref
@@ -329,6 +333,9 @@
).
In this code, nodal velocities satisfying a weak-form of the momentum balance
are obtained as the nodal solution.
+ SNAC shares a mathematical foundation, and thus major advantages, with
+ a standard finite element method (FEM).
+ However, it departs from the FEM by not making explicit use of shape functions.
A Cartesian mesh consisting of 4-node linear or constant-strain tetrahedral
elements is used to represent a discretized domain, although a spherical
domain can also be used.
@@ -348,7 +355,7 @@
\end_layout

\begin_layout Subsection
-Weak Formulation
+Governing Equations
\end_layout

\begin_layout Standard
@@ -359,16 +366,12 @@
\begin_layout Standard
\begin_inset Formula $$\begin{array}{c} -\frac{\partial\sigma_{ij}}{\partial x_{j}}+\rho g_{i}=\frac{Dv_{i}}{Dt}\,\, in\,\,\Omega\\ -v_{i}=g_{i}\,\, on\,\,\Gamma_{g}\\ -t_{i}=h_{i}\,\, on\,\,\Gamma_{h}\end{array}\label{eq:Momentum balance for continuum}$$
+\frac{\partial\sigma_{ij}}{\partial x_{j}}+\rho g_{i}=\frac{Dv_{i}}{Dt}\,\, in\,\,\Omega,\\
+v_{i}=g_{i}\,\, on\,\,\Gamma_{g},\\
+t_{i}=h_{i}\,\, on\,\,\Gamma_{h},\end{array}\label{eq:Momentum balance for continuum}

\end_inset

-
-\end_layout
-
-\begin_layout Standard
where
\begin_inset Formula $\partial\Omega$
\end_inset
@@ -401,329 +404,734 @@

is the material or total derivative, and equal to partial derivative with
respect to time in SNAC since the updated Lagrangian viewpoint is taken.
- The corresponding weak form of the momentum balance is
+
\end_layout

+\begin_layout Subsection
+Spatial Descritization
+\end_layout
+
\begin_layout Standard
-\begin_inset Formula $$-\int\delta\nu_{i,j}\sigma_{ij}d\Omega-\int\delta\nu_{i}\rho g_{i}d\Omega+\int\delta\nu_{i}\rho\frac{Dv_{i}}{Dt}d\Omega-\int_{\Gamma_{h}}\delta\nu_{i}\left(t_{i}-h_{i}\right)d\Gamma=0\label{eq:weak form of momentum balance}$$
+A 3-D domain is discretized into hexahedral elements, each of which is filled
+ with two sets of 5 tetrahedra (Fig.
+
+\begin_inset LatexCommand ref
+reference "fig:snac_mesh"

\end_inset

+a).
+ In this mesh hierarchy, called the mixed discretization
+\begin_inset ERT
+status collapsed

+\begin_layout Standard
+
+
+\backslash
+citep{MartCund1982}
\end_layout

+\end_inset
+
+, hexahedral elements are used only as an averaging unit for volumetric
+ strain.
+ The averaging is enforced at all times, for incompressible viscoelastic
+ or plastic constitutive laws.
+ The use of two equivalent sets of tetrahedra is required to ensure a symmetric
+ response.
+ For a given loading, responses of one set of tetrahedra can be different
+ from those of the other set because of the differently orientated faces
+ of tetrahedra in each set
+\begin_inset ERT
+status collapsed
+
\begin_layout Standard
-where
-\begin_inset Formula $\delta\nu_{i}$
+
+
+\backslash
+citep[e.g.,]{Zien_etal1995}
+\end_layout
+
\end_inset

- are components of virtual velocity and a comma stands for partial differentiati
-on.
- This weak form is also known as principle of virtual power.
- Introducing the isoparametric discretization for the velocity in each element
+.
\end_layout

\begin_layout Standard
-\begin_inset Formula $$-v_{i}=v_{i}^{a}N_{a}\left(X\right)\label{eq:isoparametric discretization for the velocity}$$
+\begin_inset Float figure
+wide false
+sideways false
+status open

+\begin_layout Standard
+\begin_inset ERT
+status collapsed
+
+\begin_layout Standard
+
+
+\backslash
+centering
+\end_layout
+
\end_inset

+
+\begin_inset Graphics
+	filename figures/snac-mesh.pdf
+	scale 60

+\end_inset
+
+
\end_layout

\begin_layout Standard
-where
-\begin_inset Formula $v_{i}^{a}$
+\begin_inset Caption
+
+\begin_layout Standard
+\begin_inset OptArg
+status open
+
+\begin_layout Standard
+Configurations of tetrahedra and conventions for the notation
+\end_layout
+
\end_inset

- is the
-\emph on
-i
-\emph default
--th component velocity of node
-\begin_inset Formula $a$
+(a) Two configurations of five tetrahedra in a hexahedral element used in
+ the mixed discretization.
+ Numbers next to apexes indicate the local node numbering.
+ (b) Conventions for the notation.
+
+\begin_inset Formula $A_{l}$
\end_inset

-,
-\begin_inset Formula $N_{a}$
+ and
+\begin_inset Formula $n_{l}$
\end_inset

- is the shape function associated with node
-\begin_inset Formula $a$
+ denote the face and the unit normal vector, respectively, associated with
+ a local node
+\begin_inset Formula $l$
\end_inset

-, and
-\begin_inset Formula $X$
+.
+\end_layout
+
\end_inset

- material coordinate.
- Since the discretization is isoparametric, the same discretization holds
- for virtual velocities.
- This discretization then provides a discrete version of the principle of
- virtual power:
+
\end_layout

\begin_layout Standard
-\begin_inset Formula $$-\delta\nu_{i}^{a}\left(\int\sigma_{ij}N_{a,j}d\Omega-\int N_{a}\rho g_{i}d\Omega+\int N_{a}\rho\frac{Dv_{i}}{Dt}d\Omega-\int_{\Gamma_{h}}N_{a}\left(t_{i}-h_{i}\right)d\Gamma\right)=0\label{eq:discrete version of virtual power}$$
+\begin_inset LatexCommand label
+name "fig:snac_mesh"

\end_inset

+
+\end_layout

+\end_inset
+
+
\end_layout

+\begin_layout Subsection
+Approximation of partial derivatives
+\end_layout
+
\begin_layout Standard
-Since the virtual velocities,
-\begin_inset Formula $\delta\nu_{i}$
-\end_inset
+The approximation of partial derivatives with respect to spatial variables
+ follows the integral definitions
+\begin_inset ERT
+status collapsed

-, are arbitrary everywhere except on Dirichlet boundary
+\begin_layout Standard
+
+
+\backslash
+citep[e.g.,]{Wilkins1964}
\end_layout

-\begin_layout Standard
+\end_inset
+
+:
\begin_inset Formula $$-r_{i}^{a}=\int N_{a}\rho g_{i}d\Omega-\int\sigma_{ij}N_{a,j}d\Omega-\int N_{a}\rho\frac{Dv_{i}}{Dt}d\Omega-\int_{\Gamma_{h}}N_{a}\left(t_{i}-h_{i}\right)d\Gamma=0\label{eq:5}$$
+\int_{\Omega}f_{,i}dV=\int_{\partial\Omega}fn_{i}d\Gamma,\label{eq:partialapprox}

\end_inset

+where
+\shape italic

-\end_layout
+\begin_inset Formula $\Omega$
+\end_inset

-\begin_layout Standard
+
+\shape default
+ represent a tetrahedron as an integration domain,
+\shape italic
+
+\begin_inset Formula $\partial\Omega$
+\end_inset
+
+
+\shape default
+ is the boundary surfaces of the tetrahedron,
+\begin_inset Formula $f_{,i}$
+\end_inset
+
+ is the partial derivative of a variable
+\begin_inset Formula $f$
+\end_inset
+
+ with respect to
+\begin_inset Formula $i$
+\end_inset
+
+-th spatial coordinate,
+\begin_inset Formula $n_{i}$
+\end_inset
+
+ is the
+\begin_inset Formula $i$
+\end_inset
+
+-th component of the unit normal vector of the surface.
+ If the partial derivative is constant within a tetrahedron, it is evaluated
+ as
+\begin_inset Formula $$+f_{,i}=\frac{1}{V}\int_{\partial\Omega}fn_{i}d\Gamma,\label{eq:def_fi}$$
+
+\end_inset
+
where
-\begin_inset Formula $r_{i}^{a}$
+\begin_inset Formula $V$
\end_inset

- is the residual force vector.
- The above set of equations corresponds to discrete momentum balance and
- we solve them for the nodal velocities,
-\begin_inset Formula $v_{i}^{a}$
+ is the volume of the tetrahedron.
+ By further substituting an algebraic expression for the surface integral,
+ reordering terms, and using
+\begin_inset Formula $\int_{\partial\Omega}n_{i}d\Gamma=0$
\end_inset

-.
- Because SNAC seeks a quasistatic solution, the inertia term in Eq.
-
+ (when
+\begin_inset Formula $f=1$
+\end_inset
+
+ in (
\begin_inset LatexCommand ref
-reference "eq:5"
+reference "eq:def_fi"

\end_inset

- can be dropped.
- This weak form of momentum balance equation, Eq.
+)),
+\begin_inset Formula $$+\begin{split}f_{,i} & =\frac{1}{V}\sum_{l=1}^{4}\bar{f}^{l}n_{i}^{l}A^{l}=\frac{1}{V}\sum_{l=1}^{4}\frac{1}{3}\sum_{m=1,\neq l}^{4}f^{m}n_{i}^{l}A^{l}\\ + & =\frac{1}{3V}\sum_{m=1}^{4}f^{m}\sum_{l=1,\neq m}^{4}n_{i}^{l}A^{l}\\ + & =-\frac{1}{3V}\sum_{m=1}^{4}f^{m}n_{i}^{m}A^{m},\end{split} +\label{eq:formula_fi}$$
+
+\end_inset
+
+ where
+\begin_inset Formula $l$
+\end_inset
+
+ is the local node index varying from 1 to 4,
+\begin_inset Formula $A^{l}$
+\end_inset
+
+ and
+\begin_inset Formula $n^{l}$
+\end_inset
+
+ are the area and the unit normal vector of the triangular surface not having
+ the node
+\begin_inset Formula $l$
+\end_inset
+
+ as one of its apexes (Fig.

\begin_inset LatexCommand ref
-reference "eq:5"
+reference "fig:snac_mesh"

\end_inset

-, is equivalent to that acquired from the original Finite Difference formulation
- of FLAC [Cundall, 1989].
- In the finite difference formulation, all the partial derivatives are replaced
- by the integral definitions [Wilkins, 1964], so that the partial derivatives
- of shape functions are not involved explicitly.
+b).
+ Hereafter, we call such a face a
+\emph on
+corresponding
+\emph default
+ face to node
+\begin_inset Formula $l$
+\end_inset
+
+.

+\begin_inset Formula $\bar{f}^{l}$
+\end_inset
+
+ is the averaged
+\begin_inset Formula $f$
+\end_inset
+
+ on the surface
+\begin_inset Formula $l$
+\end_inset
+
+.
\end_layout

\begin_layout Subsection
-Solution Scheme
+Nodal Assemblage
\end_layout

\begin_layout Standard
-In SNAC, the solution to velocity is obtained by an explicit version of
- dynamic relaxation.
- Consider the autonomous set of equations
-\end_layout
+We can convert the differential equation for momentum balance (
+\begin_inset LatexCommand ref
+reference "eq:Momentum balance for continuum"

-\begin_layout Standard
+\end_inset
+
+) (the following description is applied to the heat equation in the same
+ fashion) to a principle of minimum work rate as in the standard finite
+ element formulation:
\begin_inset Formula $$-M\frac{D^{2}\nu}{Dt^{2}}+C\frac{Dv}{Dt}+K\nu=f\label{eq:autonomous set of equations}$$
+\int_{\Omega}\delta v_{i}\rho\frac{Dv_{i}}{Dt}dV=\int_{\Omega}\delta v_{i}\rho g_{i}dV+\int_{\Omega}\delta\xi_{ij}\sigma_{ij}dV,\label{eq:momentumweak}

\end_inset

+where
+\shape italic

-\end_layout
+\begin_inset Formula $\xi_{ij}$
+\end_inset

-\begin_layout Standard
-where
-\series bold
-\emph on
-M
-\series default
-\emph default
-,
-\series bold
-\emph on
-C
-\series default
-\emph default
-, and
-\series bold
-\emph on
-K
-\series default
-\emph default
- are mass, damping and stiffness matrix respectively and
-\series bold
-\emph on

-\begin_inset Formula $f$
+\shape default
+ are components of the strain rate tensor,
+\shape italic
+
+\begin_inset Formula $\delta v_{i}$
\end_inset

-\series default
-\emph default
- is independent of time.
- For
-\series bold
-\emph on
-C
-\series default
-\emph default
-, positive definite, the transient velocity solution converges to that of
- equilibrium velocity solution
-\begin_inset Formula $K\nu=f$
+\shape default
+ and
+\shape italic
+
+\begin_inset Formula $\delta\xi_{ij}$
\end_inset

- as
-\begin_inset Formula $t\rightarrow\infty$
+
+\shape default
+ represent variations of velocity and strain rate, and
+\shape italic
+
+\begin_inset Formula $\Omega$
\end_inset

-.
- Critical damping is desired for the fastest convergence to equilibrium,
- but the accurate evaluation of critical damping involves the expensive
- computation of eigenvalues.
- So, we adopt a local non-viscous damping [Cundall, 1987].
- In this approach, the magnitude of the damping force on a node is proportional
- to the magnitude of the unbalanced force, and the direction of damping
- force is such that energy is always dissipated.
- The velocity update can then be written as follows:
-\end_layout

-\begin_layout Standard
+\shape default
+ here corresponds to the whole domain.
+ The local contribution to nodes corresponding to each term can be computed
+ by following the standard finite element procedure for linear tetrahedral
+ elements.
+ However, our method does not need to construct coefficient matrices such
+ as mass and stiffness matrices since it adopts an explicit time discretization.
+ The resultant momentum equation is
\begin_inset Formula $$-\nu_{ia}^{t+\frac{\Delta t}{2}}=\nu_{ia}^{t-\frac{\Delta t}{2}}+\left(r_{ia}-f_{ia}^{d}\right)\frac{\Delta t}{M_{a}}\label{eq:velocity update}$$
+M^{n}\frac{Dv_{i}^{n}}{Dt}=\frac{1}{3}T_{i}^{[n]}+\frac{1}{4}\rho^{[n]}g_{i}V^{[n]},\label{eq:momentumdiscrete}

\end_inset

+where the superscript
+\begin_inset Formula $n$
+\end_inset

-\end_layout
+ represents values evaluated at the global node
+\begin_inset Formula $n$
+\end_inset

-\begin_layout Standard
-where
-\begin_inset Formula $M_{a}$
+, the superscript
+\begin_inset Formula $[n]$
\end_inset

- is the fictitious nodal mass for node
-\begin_inset Formula $a$
+ means the sum of contributions from all the tetrahedra having the global
+ node n as an apex,
+\begin_inset Formula $T_{i}$
\end_inset

-, and
-\begin_inset Formula $f_{ia}^{d}$
+ is the traction that is defined as
+\shape italic
+
+\begin_inset Formula $\sigma_{ij}n_{j}$
\end_inset

- is the
-\emph on
-i
-\emph default
--th component damping force on node
-\begin_inset Formula $a$
+
+\shape default
+ and evaluated on a face of one of the contributing tetrahedra.
+ Any applied traction boundary conditions should be explicitly computed
+ and added on the right hand side of (
+\begin_inset LatexCommand ref
+reference "eq:momentumdiscrete"
+
\end_inset

+).
+ The nodal mass
+\begin_inset Formula $M^{n}$
+\end_inset
+
+ is not the actual inertial mass but an adjusted one to satisfy a local
+ stability criterion discussed in the section
+\begin_inset LatexCommand ref
+reference "sec:appendix_massscaling"
+
+\end_inset
+
.
- The damping force,
-\begin_inset Formula $f_{ia}^{d}$
+ The correspondence between an apex and a face for the traction calculation
+ is determined as in the derivation of the expression, (
+\begin_inset LatexCommand ref
+reference "eq:formula_fi"
+
\end_inset

-, is given by
+).
+ Note that the factor of 1/3 in the traction term is inherited from (
+\begin_inset LatexCommand ref
+reference "eq:formula_fi"
+
+\end_inset
+
+) and the factor of 1/4 in the body force term implies that the nodal contributi
+on takes one quarter of a tetrahedron's volume-dependent quantity.
\end_layout

\begin_layout Standard
+While looping over the entire set of nodes, mass and nodal forces are assembled
+ by adding up the contributions from boundary conditions and all the tetrahedra
+ sharing that node as one of their apexes.
+ The structured mesh of SNAC renders the assemblage step conveniently static.
+ The acquired net force (or heat flux) at each node is used to update velocities
+ and node coordinates (or temperature).
+\end_layout
+
+\begin_layout Subsection
+Solution Scheme
+\end_layout
+
+\begin_layout Standard
+We seek static or quasi-static solutions through a dynamic relaxation method.
+ local non-viscous damping scheme
+\begin_inset ERT
+status collapsed
+
+\begin_layout Standard
+
+
+\backslash
+citep{Cundall1987}
+\end_layout
+
+\end_inset
+
+:
\begin_inset Formula $$-f_{ia}^{d}=\alpha|r_{ia}|\mathrm{sgn}\left(\nu_{ia}^{t-\frac{\Delta t}{2}}\right)\label{eq:damping forcec}$$
+F_{i}^{damped}=F_{i}-\alpha\text{ sgn}(v_{i})|F_{i}|,\label{eq:forcedamped}

\end_inset

+ where
+\begin_inset Formula $F_{i}$
+\end_inset

-\end_layout
+ is the
+\begin_inset Formula $i$
+\end_inset

-\begin_layout Standard
-where
+-th component of the residual force vector (the right hand side of (
+\begin_inset LatexCommand ref
+reference "eq:momentumdiscrete"
+
+\end_inset
+
+)),
+\shape italic
+
\begin_inset Formula $\alpha$
\end_inset

- is the proportionality constant for damping and ??? (ends in G3 paper)
- TODO: what need to happen here?
+
+\shape default
+ is a positive coefficient less than 1, sgn
+\begin_inset Formula $(v_{i})$
+\end_inset
+
+ returns the sign of the
+\begin_inset Formula $i$
+\end_inset
+
+-th component of velocity,
+\begin_inset Formula $v_{i}$
+\end_inset
+
+.
+ Once net forces are assembled and damped, velocity at that node is updated
+ using a forward Euler method:
+\begin_inset Formula $$+v(t+\frac{\Delta t}{2})=v(t-\frac{\Delta t}{2})+\Delta t\frac{F_{i}^{damped}}{M}\label{eq:velupdate}$$
+
+\end_inset
+
+
+\begin_inset Formula $$+x(t+\Delta t)=x(t)+\Delta tv(t+\frac{\Delta t}{2}).\label{eq:posupdate}$$
+
+\end_inset
+
+ Damping is irrelevant to the update of temperature field, but the same
+ forward Euler method is used.
\end_layout

-\begin_layout Standard
-Due to the explicit nature of the solution scheme, the solution is conditionally
- stable when the numerical information speed of the mesh is more than the
- physical information speed.
- This condition can be expressed as
+\begin_layout Subsection
+Mass scaling for numerical stability
\end_layout

\begin_layout Standard
-\begin_inset Formula $-\Delta t\leq\frac{\Delta x}{C_{p}}$
+The conventional Courant-Friedrichs-Lewy (CFL) condition imposes a stringent
+ upper limit for the time step size such that dynamic relaxation takes long
+ time to get quasi-static solution over a geological time scale.
+ To overcome this limit, a mass scaling technique is applied.
+ This technique adjusts each nodal mass such that the stability condition
+ for a user-specified time step can be locally satisfied.
+ The stability condition to be satisfied, however, is not the same as in
+ the CFL condition, i.e.,
+\shape italic

+\begin_inset Formula $\Delta t$
\end_inset

+\shape default
+
+\begin_inset Formula $\leq$
+\end_inset
+
+ (
+\begin_inset Formula $l_{min}/v_{p}$
+\end_inset
+
+), where
+\shape italic
+
+\begin_inset Formula $\Delta t$
+\end_inset
+
+
+\shape default
+ is the time step,
+\begin_inset Formula $l_{min}$
+\end_inset
+
+ is the minimum element size, and
+\begin_inset Formula $v_{p}$
+\end_inset
+
+ is the P wave velocity.
+ Instead, through an analogy of continuum to an infinite mass-spring system,
+ we use a criterion that does not explicitly include length scale and P
+ wave velocity
+\begin_inset ERT
+status collapsed
+
+\begin_layout Standard
+
+
+\backslash
+citep[see Ch.
+ 9
\end_layout

\begin_layout Standard
+
+in]{Bathe1996}
+\end_layout
+
+\end_inset
+
+:
+\begin_inset Formula $$+\Delta t\leq\frac{T}{\pi},\label{eq:dtcrit}$$
+
+\end_inset
+
+ where
+\begin_inset Formula $T$
+\end_inset
+
+ is the period of system,
+\begin_inset Formula $2\pi(m/k)^{1/2}$
+\end_inset
+
+,
+\begin_inset Formula $m$
+\end_inset
+
+ is a point mass, and
+\begin_inset Formula $k$
+\end_inset
+
+ is the stiffness of the spring attached to the point mass.
+ Now, reducing the infinite series of mass and springs in one dimension
+ to a single mass-spring system, the stiffness of that single system becomes
+
+\begin_inset Formula $4k$
+\end_inset
+
+, leading to an expression for the mass scaling:
+\begin_inset Formula $$+m\geq k(\Delta t)^{2}.\label{eq:critmass}$$
+
+\end_inset
+
+ For a given size of
+\begin_inset Formula $\Delta t$
+\end_inset
+
+, the nodal mass is adjusted according to (
+\begin_inset LatexCommand ref
+reference "eq:critmass"
+
+\end_inset
+
+) to automatically satisfy the stability critetion, (
+\begin_inset LatexCommand ref
+reference "eq:dtcrit"
+
+\end_inset
+
+).
+ The value of
+\begin_inset Formula $k$
+\end_inset
+
+ is computed by equating internal force contribution at a node with
+\begin_inset Formula $-ku_{i}$
+\end_inset
+
+:
+\begin_inset Formula $$+\begin{split}\frac{1}{3}T_{i} & =-ku_{i}\Rightarrow\\ +\frac{1}{3}(\lambda+2\mu)(\dot{\epsilon}_{ii}dt)n_{i}S & =-k(v_{i}dt)\text{ (no sum)},\end{split} +\label{eq:eqfork}$$
+
+\end_inset
+
+where only the volumetric contribution from internal forces is taken into
+ account.
+ By substituting the approximation for the partial derivative (
+\begin_inset LatexCommand ref
+reference "eq:formula_fi"
+
+\end_inset
+
+) into the above equation and dividing both sides by
+\begin_inset Formula $v_{i}dt$
+\end_inset
+
+, we obtain
+\begin_inset Formula $$+k_{i}^{l}=\frac{1}{9V}(\lambda+2\mu)(n_{i}^{l}S^{l})^{2},\label{eq:kdef}$$
+
+\end_inset
+
where
-\begin_inset Formula $C_{p}=\sqrt{\frac{\kappa+4G/3}{\rho}}$
+\begin_inset Formula $l$
\end_inset

- is the P-wave speed and
-\begin_inset Formula $\Delta x$
+ is the local index for apexes of a tetrahedron and the surface-related
+ quantities are computed on the corresponding face of the tetrahedron.
+ Finally, a tetrahedron's contribution to the scaled mass is given as
+\begin_inset Formula $$+m^{l}=\frac{\lambda+2\mu}{9V}\max[(n_{i}^{l}S^{l})^{2},i=1,\dots,3].\label{eq:scaledmass}$$
+
\end_inset

- is the element size.
- Since the mass matrix is fictitious, we can choose a mass matrix in such
- a way that convergence is optimal.
- The optimal convergence is obtained when local values of critical time
- steps for each element are equal.
- The condition for this can be obtained as
+
\end_layout

\begin_layout Standard
-\begin_inset Formula $$-M_{ia}=\sum\frac{\left(\kappa+4G/3\right)\Delta x_{\mathrm{max}}^{2}}{6V}\label{eq:condition for optimal convergence}$$
+As in the standard FEM, appropriate mappings between local and global indices
+ are required.
+\end_layout

-\end_inset
+\begin_layout Subsection
+Constitutive update
+\end_layout

+\begin_layout Standard
+SNAC uses a general elasto-visco-plastic rheological model to update the
+ Cauchy stress tensor
+\begin_inset ERT
+status collapsed

+\begin_layout Standard
+
+
+\backslash
+citep[e.g.,]{Albe_etal2000}
\end_layout

-\begin_layout Standard
-where summation is over all the tetrahedral elements sharing the node a,
- and the
-\begin_inset Formula $\Delta x_{\mathrm{max}}$
\end_inset

- is the maximum propagation distance.
+.
+ First, the initial guess of stress is acquired by the Maxwell viscoelastic
+ constitutive law
+\begin_inset ERT
+status collapsed
+
+\begin_layout Standard
+
+
+\backslash
+citep{Poli_etal1993a}
\end_layout

-\begin_layout Subsection
-Constitutive Update
+\end_inset
+
+.
+ If this initial guess exceeds a given yield stress, it is projected onto
+ the yield surface using a return mapping method
+\begin_inset ERT
+status collapsed
+
+\begin_layout Standard
+
+
+\backslash
+citep{SimoHugh2004}
\end_layout

-\begin_layout Standard
-SNAC uses a general elasto-visco-plastic (or viscoplastic) rheological model
- to update Cauchy stress based on strain-rate.
- This viscoplastic model can deal with various constitutive laws that are
- typically used for the Earth's crustal and mantle as its limiting cases.
+\end_inset
+
+; otherwise, the viscoelastic stress update is retained.
+ This elasto-visco-plastic model can deal with various constitutive laws
+ that are typically used for the Earth's crustal and mantle material as
+ its limiting cases.
For example, elastic, viscoelastic and elastoplastic material are realized
in the following cases:
\end_layout

\begin_layout Enumerate
-Elastic material corresponds to the limit of infinite viscosity and infinite
- yield stress.
+Elastic material corresponds to the limit of infinite viscosity and yield
+ stress.

\end_layout

@@ -733,76 +1141,67 @@
\end_layout

\begin_layout Enumerate
-Elasoplastic material corresponds to the infinite viscosity.
+Elasoplastic material corresponds to the inifinte viscosity.

\end_layout

\begin_layout Standard
Using the viscoplastic rheology is physically more realistic than using
one of the limiting cases listed above since all materials have dissipative
- mechanisms, however small they might be, and hence have viscosity.
+ mechanisms and hence viscosity.
This viscosity also provides a length scale for the problem of localization,
which in turn enables physically meaningful mesh independent solution when
- the mesh size is smaller than this length scale [e.g., Needleman, 1988].
-
+ the mesh size is smaller than this length scale.
\end_layout

\begin_layout Standard
Since the nodal variables are velocities and whose spatial gradients are
- deformation rates, we formulate constitutive update in terms of strain
+ deformation rates, we formulate the constitutive update in terms of strain
rate.
- The objective stress rate of choice is Jaumann or the corotational stress
- rate (
-\begin_inset Formula $\sigma^{\Delta J}$
+ The objective stress rate of choice is the Jaumann or the corotational
+ stress rate (
+\begin_inset Formula $\Delta\sigma^{\Delta J}$
\end_inset

- ),
-\begin_inset Formula $\sigma^{\Delta J}=C^{\Delta J}D$
-\end_inset
+)
+\begin_inset ERT
+status collapsed

-
-\end_layout
-
\begin_layout Standard
-where
-\end_layout

-\begin_layout Standard
-\begin_inset Formula $$-D\,(D_{i,j}=\left(\frac{\partial\nu_{i}}{\partial x_{j}}+\frac{\partial\nu_{j}}{\partial x_{i}}\right)),\,\, W\,(W_{i,j}=\frac{1}{2}\left(\frac{\partial\nu_{i}}{\partial x_{j}}-\frac{\partial\nu_{j}}{\partial x_{i}}\right))\label{eq:deformation and spin tensors}$$

-\end_inset
+\backslash
+citep{RudnRice1975}
+\end_layout

- are the deformation and spin tensors and
-\begin_inset Formula $C^{\Delta J}$
\end_inset

- is Jaumann tangent modulus.
- For the isotropic material Jaumann tangent modulus is the same as the elastic
- tangent modulus, in which case correction to the stresses due to rotation
- can be given as
-\end_layout
-
-\begin_layout Standard
+
\begin_inset Formula $$-\Delta\sigma=\Delta\sigma^{C}+\Delta\sigma^{R}=\Delta\sigma^{C}+\left(W\sigma-\sigma W\right)\Delta t,\label{eq:lots of deltas}$$
+\Delta\sigma^{\Delta J}=\frac{\partial(\Delta\sigma)}{\partial t}-W\cdot\Delta\sigma-\Delta\sigma\cdot W^{T},\label{eq:corotstress}

\end_inset

+ where
+\begin_inset Formula $W_{ij}=(1/2)(\partial v_{i}/\partial x_{j}$
+\end_inset

-\end_layout
+-
+\begin_inset Formula $\partial v_{j}/\partial x_{i})$
+\end_inset

-\begin_layout Standard
-where
-\begin_inset Formula $\Delta\sigma^{C}$
+ are the components of spin tensor and
+\begin_inset Formula $\Delta\sigma$
\end_inset

- denotes the stress update by constitutive law and
-\begin_inset Formula $\Delta\sigma^{R}$
+ is the increment of stress tensor.
+ Correction to the stresses due to rotation can be given as
+\begin_inset Formula $$+\sigma^{t+\Delta t}=\sigma^{t}+\Delta\sigma^{\Delta J}.\Delta t\label{eq:stressupdate}$$
+
\end_inset

- stands for the correction due to finite rotation, of the stress update.
-
+
\end_layout

\begin_layout Standard
@@ -897,8 +1296,8 @@
\begin_inset Formula $\sigma^{n+1}$
\end_inset

- onto the yield surface using a return-mapping algorithm [TODO need citation:
- Ortiz and Simo, 1986].
+ onto the yield surface using a return-mapping algorithm [Simo and Hughes,
+ 1999].

\end_layout

@@ -979,8 +1378,7 @@
\begin_layout Standard
In general, flow rule for frictional materials is non-associative, i.e., flow
direction differs from the normal of the yield surface normal.
- The plastic flow potential for plastic flow [TODO -- is there some duplication
- in the foregoing phrase?] can be given as
+ The plastic flow potential can be given as
\end_layout

\begin_layout Standard
@@ -1314,8 +1712,7 @@

\begin_layout Standard
SNAC is built up from the StGermain framework, an environment for the developmen
-t of physically based scientific codes [Quenette, et al.
- TODO--need year].
+t of physically based scientific codes [Quenette, et al., 2005].
As frameworks, StGermain and Pyre differ in that StGermain provides the
infrastructural components needed to build complete codes while Pyre provides
the superstructure to couple codes together; for this reason the former