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

[echoi at geodynamics.org]

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

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



Added: doc/snac/figures/snac-mesh.pdf
===================================================================
(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 @@
-#LyX 1.5.1 created this file. For more info see http://www.lyx.org/
+#LyX 1.5.3 created this file. For more info see http://www.lyx.org/
 \lyxformat 276
 \begin_document
 \begin_header
@@ -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
 reference "cha:License"
@@ -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{equation}
 \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}\end{equation}
+\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{equation}
 
 \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 \begin{equation}
-\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}\end{equation}
+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 \begin{equation}
-v_{i}=v_{i}^{a}N_{a}\left(X\right)\label{eq:isoparametric discretization for the velocity}\end{equation}
+\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 \begin{equation}
-\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}\end{equation}
+\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 \begin{equation}
-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}\end{equation}
+\int_{\Omega}f_{,i}dV=\int_{\partial\Omega}fn_{i}d\Gamma,\label{eq:partialapprox}\end{equation}
 
 \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 \begin{equation}
+f_{,i}=\frac{1}{V}\int_{\partial\Omega}fn_{i}d\Gamma,\label{eq:def_fi}\end{equation}
+
+\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{equation}
+\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{equation}
+
+\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 \begin{equation}
-M\frac{D^{2}\nu}{Dt^{2}}+C\frac{Dv}{Dt}+K\nu=f\label{eq:autonomous set of equations}\end{equation}
+\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{equation}
 
 \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 \begin{equation}
-\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}\end{equation}
+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{equation}
 
 \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.
+ Instead of adding a usual velocity-dependent friction term, we adopt a
+ local non-viscous damping scheme 
+\begin_inset ERT
+status collapsed
+
+\begin_layout Standard
+
+
+\backslash
+citep{Cundall1987}
+\end_layout
+
+\end_inset
+
+: 
 \begin_inset Formula \begin{equation}
-f_{ia}^{d}=\alpha|r_{ia}|\mathrm{sgn}\left(\nu_{ia}^{t-\frac{\Delta t}{2}}\right)\label{eq:damping forcec}\end{equation}
+F_{i}^{damped}=F_{i}-\alpha\text{ sgn}(v_{i})|F_{i}|,\label{eq:forcedamped}\end{equation}
 
 \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 \begin{equation}
+v(t+\frac{\Delta t}{2})=v(t-\frac{\Delta t}{2})+\Delta t\frac{F_{i}^{damped}}{M}\label{eq:velupdate}\end{equation}
+
+\end_inset
+
+ 
+\begin_inset Formula \begin{equation}
+x(t+\Delta t)=x(t)+\Delta tv(t+\frac{\Delta t}{2}).\label{eq:posupdate}\end{equation}
+
+\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 \begin{equation}
+\Delta t\leq\frac{T}{\pi},\label{eq:dtcrit}\end{equation}
+
+\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 \begin{equation}
+m\geq k(\Delta t)^{2}.\label{eq:critmass}\end{equation}
+
+\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{equation}
+\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{equation}
+
+\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 \begin{equation}
+k_{i}^{l}=\frac{1}{9V}(\lambda+2\mu)(n_{i}^{l}S^{l})^{2},\label{eq:kdef}\end{equation}
+
+\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 \begin{equation}
+m^{l}=\frac{\lambda+2\mu}{9V}\max[(n_{i}^{l}S^{l})^{2},i=1,\dots,3].\label{eq:scaledmass}\end{equation}
+
 \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 \begin{equation}
-M_{ia}=\sum\frac{\left(\kappa+4G/3\right)\Delta x_{\mathrm{max}}^{2}}{6V}\label{eq:condition for optimal convergence}\end{equation}
+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 \begin{equation}
-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{equation}
 
-\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 \begin{equation}
-\Delta\sigma=\Delta\sigma^{C}+\Delta\sigma^{R}=\Delta\sigma^{C}+\left(W\sigma-\sigma W\right)\Delta t,\label{eq:lots of deltas}\end{equation}
+\Delta\sigma^{\Delta J}=\frac{\partial(\Delta\sigma)}{\partial t}-W\cdot\Delta\sigma-\Delta\sigma\cdot W^{T},\label{eq:corotstress}\end{equation}
 
 \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 \begin{equation}
+\sigma^{t+\Delta t}=\sigma^{t}+\Delta\sigma^{\Delta J}.\Delta t\label{eq:stressupdate}\end{equation}
+
 \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