[cig-commits] commit: Added short paragraph on role of CIG. Cleanup of notation.
Mercurial
hg at geodynamics.org
Thu Sep 1 13:40:16 PDT 2011
changeset: 66:0e97c1c8a9ef
tag: tip
user: Brad Aagaard <baagaard at usgs.gov>
date: Thu Sep 01 13:40:13 2011 -0700
files: faultRup.tex references.bib
description:
Added short paragraph on role of CIG. Cleanup of notation.
diff -r 736e9ef822bf -r 0e97c1c8a9ef faultRup.tex
--- a/faultRup.tex Thu Sep 01 12:58:29 2011 -0700
+++ b/faultRup.tex Thu Sep 01 13:40:13 2011 -0700
@@ -153,8 +153,16 @@ used in the dynamic simulations; this al
used in the dynamic simulations; this also results in using different
solvers as we will discuss later.
-\brad{TODO: Add paragraph on CIG (open-source, community code so needs
- to be robust, flexible, easy to maintain)}
+Additional motivation for developing PyLith arose from the geophysics
+community as part of the Computational Infrastructure for Geodynamics
+(CIG) project \cite{CIG:web:page}. CIG provides guidelines for
+developing robust, open-source code as well as a forum for gathering
+feature requests from the community. Serving the broad needs of the
+community with limited resources generated further incentives for
+designing PyLith to leverage common infrastructure for simulation
+quasi-static deformation and dynamic deformation. Maintaining two
+seperate code bases would a require considerably greater development
+effort.
Implementing slip on the potentially nonplanar fault surface
differentiates these types of problems from many other elasticity
@@ -183,7 +191,7 @@ We solve the elasticity equation includi
- \tensor{\nabla} \cdot \tensor{\sigma} = \vec{0} \text{ in }V, \\
\tensor{\sigma} \cdot \vec{n} = \vec{T} \text{ on }S_T, \\
\vec{u} = \vec{u}_0 \text{ on }S_u, \\
- \vec{d} - \tensor{R} \left( \vec{u}_{+} - \vec{u}_{-}\right) = \vec{0}
+ \tensor{R} \cdot \vec{d} - \left( \vec{u}_{+} - \vec{u}_{-}\right) = \vec{0}
\text{ on }S_f,
\end{gather}
where $\vec{u}$ is the displacement vector, $\rho$ is the mass
@@ -191,7 +199,7 @@ stress tensor, and $t$ is time. We speci
stress tensor, and $t$ is time. We specify tractions $\vec{T}$ on
surface $S_T$, displacements $\vec{u_0}$ on surface $S_u$, and slip
$\vec{d}$ on fault surface $S_f$. The rotation matrix $\tensor{R}$
-transforms vectors from the global coordinate system to the fault
+transforms vectors from the fault coordinate system to the global
coordinate system. Because both $\vec{T}$ and $\vec{u}$ are vector
quantities, there can be some spatial overlap of the surfaces $S_T$
and $S_u$; however, a degree of freedom cannot be associated with both
@@ -248,11 +256,9 @@ surface to zero,
surface to zero,
\begin{equation}
\int_{S_f} \vec{\phi} \cdot
- \left( \tensor{R}_T \vec{d} - \vec{u}_{+} + \vec{u}_{-} \right) \, dS = 0.
+ \left( \tensor{R} \cdot \vec{d} - \vec{u}_{+} + \vec{u}_{-} \right) \, dS = 0.
\end{equation}
-The rotation matrix $\tensor{R}_T$ is the inverse of
-$\tensor{R}$ and is composed of blocks that are the transpose of the
-corresponding blocks in $\tensor{R}$.
+The rotation matrix $\tensor{R}$ is comprised of blocks of direction cosines.
We express the weighting function $\vec{\phi}$, trial solution
$\vec{u}$, and Lagrange multipliers $\vec{l}$ as linear combinations
@@ -271,26 +277,29 @@ combination of basis functions in terms
combination of basis functions in terms of a matrix-vector product, we
have
\begin{gather}
-\vec{\phi} = \tensor{N}_m \vec{a}_m, \\
-\vec{u} = \tensor{N}_n \vec{u}_n, \\
-\vec{l} = \tensor{N}_p \vec{l}_p.
+\vec{\phi} = \tensor{N}_m \cdot \vec{a}_m, \\
+\vec{u} = \tensor{N}_n \cdot \vec{u}_n, \\
+\vec{l} = \tensor{N}_p \cdot \vec{l}_p.
\end{gather}
The weighting function is
arbitrary, so the integrands must be zero for all $\vec{a}_m$, which
leads to
\begin{gather}
+ \begin{split}
- \int_{V} \nabla \tensor{N}_m^T \cdot \tensor{\sigma} \, dV
-+ \int_{S_T} \tensor{N}_m^T \vec{T} \, dS
-+ \int_{S_f^{+}} \tensor{N}_m^T \tensor{N}_p \vec{l}_p \, dS
-- \int_{S_f^{-}} \tensor{N}_m^T \tensor{N}_p \vec{l}_p \, dS
-+ \int_{V} \tensor{N}_m^T \vec{f} \, dV
-- \int_{V} \rho \tensor{N}_m^T \tensor{N}_n \frac{\partial^2 \vec{u}_n}{\partial
++ \int_{S_T} \tensor{N}_m^T \cdot \vec{T} \, dS
++ \int_{S_f^{+}} \tensor{N}_m^T \cdot \tensor{N}_p \cdot \vec{l}_p \, dS
+- \int_{S_f^{-}} \tensor{N}_m^T \cdot \tensor{N}_p \cdot \vec{l}_p \, dS
++ \int_{V} \tensor{N}_m^T \cdot \vec{f} \, dV \\
+- \int_{V} \rho \tensor{N}_m^T \cdot \tensor{N}_n \cdot \frac{\partial^2 \vec{u}_n}{\partial
t^2} \, dV
-=\vec{0}, \\
-%
+=\vec{0},
+\end{split}
+\\
+%
\int_{S_f} \tensor{N}_p^T
- \left( \tensor{R}_T \vec{d}
+ \left( \tensor{R} \cdot \vec{d}
- \tensor{N}_{n^+} \vec{u}_{n^+} + \tensor{N}_{n^-} \vec{u}_{n^-}
\right)
\, dS = \vec{0}.
@@ -318,19 +327,19 @@ and the loading conditions. Considering
\label{eqn:quasi-static:residual:domain}
\begin{split}
- \int_{V} \nabla \tensor{N}_m^T \cdot \tensor{\sigma}(\tpdt) \, dV
- + \int_{S_T} \tensor{N}_m^T \vec{T}(\tpdt) \, dS
- + \int_{S_f^{+}} \tensor{N}_m^T \tensor{N}_p \vec{l}_p(\tpdt) \, dS
- - \int_{S_f^{-}} \tensor{N}_m^T \tensor{N}_p \vec{l}_p(\tpdt) \, dS
+ + \int_{S_T} \tensor{N}_m^T \cdot \vec{T}(\tpdt) \, dS
+ + \int_{S_f^{+}} \tensor{N}_m^T \cdot \tensor{N}_p \cdot \vec{l}_p(\tpdt) \, dS
+ - \int_{S_f^{-}} \tensor{N}_m^T \cdot \tensor{N}_p \cdot \vec{l}_p(\tpdt) \, dS
\\
- + \int_{V} \tensor{N}_m^T \vec{f}(\tpdt) \, dV
+ + \int_{V} \tensor{N}_m^T \cdot \vec{f}(\tpdt) \, dV
=\vec{0},
\end{split}
\\
%
\label{eqn:quasi-static:residual:fault}
\int_{S_f} \tensor{N}_p^T
- \left( \tensor{R}_T \vec{d}(\tpdt)
- - \tensor{N}_{n^+} \vec{u}_{n^+}(\tpdt) + \tensor{N}_{n^-} \vec{u}_{n^-}(\tpdt)
+ \left( \tensor{R} \cdot \vec{d}(\tpdt)
+ - \tensor{N}_{n^+} \cdot \vec{u}_{n^+}(\tpdt) + \tensor{N}_{n^-} \cdot \vec{u}_{n^-}(\tpdt)
\right)
\, dS = \vec{0}.
\end{gather}
@@ -351,16 +360,16 @@ and~(\ref{eqn:quasi-static:residual:faul
\vec{r} = \left( \begin{array}{c}
\begin{aligned}
- \int_{V} \nabla \tensor{N}_m^T \cdot \tensor{\sigma}(\tpdt) \, dV
- + \int_{S_T} \tensor{N}_m^T \vec{T}(\tpdt) \, dS
- + \int_{S_f^{+}} \tensor{N}_m^T \tensor{N}_p \vec{l}_p(\tpdt) \, dS
- &- \int_{S_f^{-}} \tensor{N}_m^T \tensor{N}_p \vec{l}_p(\tpdt) \, dS
+ + \int_{S_T} \tensor{N}_m^T \cdot \vec{T}(\tpdt) \, dS
+ + \int_{S_f^{+}} \tensor{N}_m^T \cdot \tensor{N}_p \cdot \vec{l}_p(\tpdt) \, dS
+ &- \int_{S_f^{-}} \tensor{N}_m^T \cdot \tensor{N}_p \cdot \vec{l}_p(\tpdt) \, dS
\\
- &+ \int_{V} \tensor{N}_m^T \vec{f}(\tpdt) \, dV
+ &+ \int_{V} \tensor{N}_m^T \cdot \vec{f}(\tpdt) \, dV
\end{aligned}
\\
\int_{S_f} \tensor{N}_p^T
- \left( \tensor{R}_T \vec{d}(\tpdt)
- - \tensor{N}_{n^+} \vec{u}_{n^+}(\tpdt) + \tensor{N}_{n^-} \vec{u}_{n^-}(\tpdt)
+ \left( \tensor{R} \cdot \vec{d}(\tpdt)
+ - \tensor{N}_{n^+} \cdot \vec{u}_{n^+}(\tpdt) + \tensor{N}_{n^-} \cdot \vec{u}_{n^-}(\tpdt)
\right)
\, dS
\end{array} \right).
@@ -393,7 +402,7 @@ associated with equation~(\ref{eqn:quasi
associated with equation~(\ref{eqn:quasi-static:residual:fault}) is
\begin{equation}
\label{eqn:jacobian:constraint}
- \tensor{L} = \int_{S_f} \tensor{N}_p^T (-\tensor{N}_{n^+}+\tensor{N}_{n^-}) \, dS.
+ \tensor{L} = \int_{S_f} \tensor{N}_p^T \cdot (-\tensor{N}_{n^+}+\tensor{N}_{n^-}) \, dS.
\end{equation}
Thus, the Jacobian of the entire system has the form,
\begin{equation}
@@ -417,18 +426,18 @@ inertial term:
\vec{r} = \left( \begin{array}{c}
\begin{aligned}
- \int_{V} \nabla \tensor{N}_m^T \cdot \tensor{\sigma}(\tpdt) \, dV
- + \int_{S_T} \tensor{N}_m^T \vec{T}(\tpdt) \, dS
- &+ \int_{S_f^{+}} \tensor{N}_m^T \tensor{N}_p \vec{l}_p(\tpdt) \, dS
- - \int_{S_f^{-}} \tensor{N}_m^T \tensor{N}_p \vec{l}_p(\tpdt) \, dS
+ + \int_{S_T} \tensor{N}_m^T \cdot \vec{T}(\tpdt) \, dS
+ &+ \int_{S_f^{+}} \tensor{N}_m^T \cdot \tensor{N}_p \cdot \vec{l}_p(\tpdt) \, dS
+ - \int_{S_f^{-}} \tensor{N}_m^T \cdot \tensor{N}_p \cdot \vec{l}_p(\tpdt) \, dS
\\
- &+ \int_{V} \tensor{N}_m^T \vec{f}(\tpdt) \, dV
- - \int_{V} \rho \tensor{N}_m^T \tensor{N}_n
+ &+ \int_{V} \tensor{N}_m^T \cdot \vec{f}(\tpdt) \, dV
+ - \int_{V} \rho \tensor{N}_m^T \cdot \tensor{N}_n \cdot
\frac{\partial^2 \vec{u}_n}{\partial t^2} \, dV
\end{aligned}
\\
\int_{S_f} \tensor{N}_p^T
- \left( \tensor{R}_T \vec{d}(\tpdt)
- - \tensor{N}_{n^+} \vec{u}_{n^+}(\tpdt) + \tensor{N}_{n^-} \vec{u}_{n^-}(\tpdt)
+ \left( \tensor{R} \cdot \vec{d}(\tpdt)
+ - \tensor{N}_{n^+} \cdot \vec{u}_{n^+}(\tpdt) + \tensor{N}_{n^-} \cdot \vec{u}_{n^-}(\tpdt)
\right)
\, dS
\end{array} \right).
@@ -452,17 +461,17 @@ difference scheme, so that the accelerat
\end{gather}
Expanding the inertial term yields
\begin{equation}
- - \int_{V} \rho \tensor{N}_m^T \tensor{N}_n \frac{\partial^2 \vec{u}_n}{\partial
+ - \int_{V} \rho \tensor{N}_m^T \cdot \tensor{N}_n \cdot \frac{\partial^2 \vec{u}_n}{\partial
t^2} \, dV =
- - \frac{1}{\Delta t^2} \int_{V} \rho \tensor{N}_m^T \tensor{N}_n \left(
- d\vec{u}_n(t) - \vec{u}_n(t) + \vec{u}_n(\tmdt)
- \right) \, dV,
+ - \frac{1}{\Delta t^2} \int_{V} \rho \tensor{N}_m^T \cdot
+ \tensor{N}_n \cdot
+ \left( d\vec{u}_n(t) - \vec{u}_n(t) + \vec{u}_n(\tmdt) \right) \, dV,
\end{equation}
so that the upper portion of the Jacobian is
\begin{equation}
\label{eqn:jacobian:explicit:inertia}
\tensor{K} =
- \frac{1}{\Delta t^2} \int_{V} \rho \tensor{N}_m^T \tensor{N}_n \, dV.
+ \frac{1}{\Delta t^2} \int_{V} \rho \tensor{N}_m^T\ \cdot \tensor{N}_n \, dV.
\end{equation}
@@ -558,9 +567,12 @@ corresponding to a perturbation in the L
corresponding to a perturbation in the Lagrange multipliers $\partial
\vec{l}_p$:
\begin{gather}
- \tensor{K}_{n^+n^+} \partial \vec{u}_{n^+} = - \tensor{L}_p^T \partial \vec{l}_p, \\
- \tensor{K}_{n^-n^-} \partial \vec{u}_{n^-} = \tensor{L}_p^T \partial \vec{l}_p, \\
- \partial \vec{d} = \tensor{R} \left( \partial \vec{u}_{n^+} - \partial \vec{u}_{n^-} \right).
+ \tensor{K}_{n^+n^+} \cdot \partial \vec{u}_{n^+} =
+ - \tensor{L}_p^T \cdot \partial \vec{l}_p, \\
+ \tensor{K}_{n^-n^-} \cdot \partial \vec{u}_{n^-} =
+ \tensor{L}_p^T \cdot \partial \vec{l}_p, \\
+ \partial \vec{d} = \tensor{R}^{-1} \cdot
+ \left( \partial \vec{u}_{n^+} - \partial \vec{u}_{n^-} \right).
\end{gather}
This estimate is exact if all other degrees of freedom are constrained
(which is generally only the case for simple toy problems), quite good
@@ -578,11 +590,13 @@ normal traction $T_m$,
\begin{equation}
T_f = T_c - \mu_f T_m.
\end{equation}
-We use the sign convention that a compressive normal tractions are
-negative. The fault constitutive models include static friction,
-linear slip-weakening, linear time-weakening, and Dietrich-Ruina
-rate-state friction with an aging law (see \citeN{PyLith:manual:1.6.2}
-for details).
+$T_f$ in the this equation corresponds to the magnitude of the shear
+traction vector; the shear traction vector is resolved into the
+direction of the slip rate. We use the sign convention that a
+compressive normal tractions are negative. The fault constitutive
+models include static friction, linear slip-weakening, linear
+time-weakening, and Dietrich-Ruina rate-state friction with an aging
+law (see \citeN{PyLith:manual:1.6.2} for details).
% ------------------------------------------------------------------
\section{Finite-Element Mesh Processing}\matt{Rough draft}
@@ -655,12 +669,12 @@ integral,
integral,
\begin{equation}
\tensor{A} \cdot \vec{u}_\mathrm{rigid} =
- \tensor{A}_\mathit{diagonal} \vec{u}_\mathrm{rigid}.
+ \tensor{A}_\mathit{diagonal} \cdot \vec{u}_\mathrm{rigid}.
\end{equation}
Expressing the diagonal block of the Jacobian matrix as a vector and
the matrix of basis functions as a vector we have,
\begin{equation}
- \tensor{A} = \int_\Omega \tensor{N}^T \tensor{N} \, d\Omega \rightarrow
+ \tensor{A} = \int_\Omega \tensor{N}^T \cdot \tensor{N} \, d\Omega \rightarrow
\tensor{A}_\mathit{diagonal} = \int_\Omega \vec{N} \sum_i N_i \, d\Omega,
\end{equation}
where $N_i$ is the scalar basis function for degree of freedom $i$ and
@@ -700,7 +714,7 @@ residual assuming the increment in the s
$\vec{u}_n^*$ ignoring the off-diagonal blocks in the Jacobian and
the increment in the Lagrange multipliers.
\begin{equation}
-d\vec{u}_n^* = \tensor{K}^{-1} \vec{r}_n,
+d\vec{u}_n^* = \tensor{K}^{-1} \cdot \vec{r}_n,
\end{equation}
taking advantage of the fact that we construct $\tensor{K}$ so that it
is diagonal.
@@ -724,37 +738,41 @@ each side of the fault, we have
each side of the fault, we have
\begin{gather}
d\vec{u}_{n^+} =
- d\vec{u}_{n^+}^* - \tensor{K}_{n^+n^+}^{-1} \tensor{L}_p^T d\vec{l}_p, \\
+ d\vec{u}_{n^+}^* - \tensor{K}_{n^+n^+}^{-1} \cdot \tensor{L}_p^T \cdot d\vec{l}_p, \\
d\vec{u}_{n^-} =
- d\vec{u}_{n^-}^* + \tensor{K}_{n^-n^-}^{-1} \tensor{L}_p^T d\vec{l}_p.
+ d\vec{u}_{n^-}^* + \tensor{K}_{n^-n^-}^{-1} \cdot \tensor{L}_p^T \cdot d\vec{l}_p.
\end{gather}
Substituting into the second row of
equation~(\ref{eqn:lumped:jacobian:residual}) and isolating the term
with the increment in the Lagrange multipliers yields
\begin{equation}
- \tensor{L}_p \left( \tensor{K}_{n^+n^+}^{-1} + \tensor{K}_{n^-n^-}^{-1} \right)
- \tensor{L}^T_p d\vec{l}_p =
- -\vec{r}_p^* + \tensor{L}_p \left( d\vec{u}_{n^+}^* - d\vec{u}_{n^-}^* \right).
+ \tensor{L}_p \cdot
+ \left( \tensor{K}_{n^+n^+}^{-1} + \tensor{K}_{n^-n^-}^{-1} \right) \cdot
+ \tensor{L}^T_p \cdot d\vec{l}_p =
+ -\vec{r}_p^* + \tensor{L}_p \cdot
+ \left( d\vec{u}_{n^+}^* - d\vec{u}_{n^-}^* \right).
\end{equation}
Letting
\begin{equation}
- \tensor{S}_p = \tensor{L}_p
- \left( \tensor{K}_{n^+n^+}^{-1} + \tensor{K}_{n^-n^-}^{-1} \right)
+ \tensor{S}_p = \tensor{L}_p \cdot
+ \left( \tensor{K}_{n^+n^+}^{-1} + \tensor{K}_{n^-n^-}^{-1} \right) \cdot
\tensor{L}^T_p,
\end{equation}
-and recognizing that $S_p$ is diagonal because ${K}$ and ${L}_p$ are
+and recognizing that $\tensor{S}_p$ is diagonal because ${K}$ and ${L}_p$ are
diagonal allows us to solve for the increment in the Lagrange
multipliers,
\begin{equation}
- d\vec{l}_p = \tensor{S}_p^{-1} \left(
- -\vec{r}_p^* + \tensor{L}_p \left( d\vec{u}_{n^+}^* - d\vec{u}_{n^-}^* \right)
+ d\vec{l}_p = \tensor{S}_p^{-1} \cdot \left(
+ -\vec{r}_p^* + \tensor{L}_p \cdot
+ \left( d\vec{u}_{n^+}^* - d\vec{u}_{n^-}^* \right)
\right).
\end{equation}
Now that we have the increment in the Lagrange multipliers, we can
correct our initial solution $d\vec{u}_n^*$ so that the true residual
is zero,
\begin{equation}
- d\vec{u}_n = d\vec{u}_n^* - \tensor{K}^{-1} \tensor{L}^T d\vec{l}_p.
+ d\vec{u}_n =
+ d\vec{u}_n^* - \tensor{K}^{-1} \cdot \tensor{L}^T \cdot d\vec{l}_p.
\end{equation}
Because $\tensor{K}$ and $\tensor{L}$ are comprised of diagonal
blocks, this expression for the updates to the solution are local to
@@ -771,9 +789,9 @@ the Lagrange multipliers
(equation~(\ref{eqn:spontaneous:rupture:slip:update})) is exact and
simplifies to
\begin{equation}
- \partial \vec{d} = - \tensor{R}
+ \partial \vec{d} = - \tensor{R}^{-1} \cdot
\left( \tensor{K}_{n^+n^+}^{-1} + \tensor{K}_{n^-n^-}^{-1} \right)
- \tensor{L}_p^T \partial \vec{l}_p.
+ \cdot \tensor{L}_p^T \cdot \partial \vec{l}_p.
\end{equation}
diff -r 736e9ef822bf -r 0e97c1c8a9ef references.bib
--- a/references.bib Thu Sep 01 12:58:29 2011 -0700
+++ b/references.bib Thu Sep 01 13:40:13 2011 -0700
@@ -748,3 +748,12 @@
year = 2011,
note = {http://www.geodynamics.org/cig/software/pylith/pylith\_manual-1.6.2.pdf}
}
+
+
+ at Misc{CIG:web:page,
+ Author = "Kellogg, L.",
+ Title = "{CIG} {Web} page",
+ Note = "http://geodynamics.org",
+ Year = 2011,
+}
+
More information about the CIG-COMMITS
mailing list