[cig-commits] commit: Merge
Mercurial
hg at geodynamics.org
Thu Sep 1 18:33:54 PDT 2011
changeset: 72:26022948ab53
tag: tip
parent: 71:4dde491c41e8
parent: 69:a1b4751101fa
user: Matthew G. Knepley <knepley at gmail.com>
date: Thu Sep 01 20:33:15 2011 -0500
files: faultRup.tex references.bib
description:
Merge
diff -r 4dde491c41e8 -r 26022948ab53 faultRup.tex
--- a/faultRup.tex Thu Sep 01 13:17:57 2011 -0500
+++ b/faultRup.tex Thu Sep 01 20:33:15 2011 -0500
@@ -1,5 +1,5 @@
-%\documentclass[draftstyle]{bssa}
-\documentclass{bssa}
+%\documentclass[draftstyle]{bssa} % 2 column, BSSA print style
+\documentclass{bssa} % double-spaced, submission style
\usepackage{array}
\usepackage{epsfig}
@@ -22,7 +22,8 @@
\clubpenalty=999999
\brokenpenalty=999999
-\title{ADD TITLE HERE}
+\title{A Domain Decomposition Approach to Implementing Fault Slip in
+ Finite-Element Models of Quasi-static and Dynamic Crustal Deformation}
\author{Brad T. Aagaard, Matthew G. Knepley, Charles A. Williams}
\affiliation{U.S. Geological Survey, MS977\\
@@ -62,7 +63,7 @@ this vast range of scales generally lead
this vast range of scales generally leads most researchers to focus on
a narrow space-time window in order to isolate just one or a few
processes; the limited spatial and temporal coverage of observations
-also justifies this narrow focus.
+also often justifies this narrow focus.
Researchers have recognized for some time, though, that interseismic
deformation and fault interactions influence earthquake rupture
@@ -111,20 +112,20 @@ deformation, they examined the effects o
deformation, they examined the effects of low-rigidity layers and a
fault damaged zone on rupture dynamics. In addition to purely dynamic
effects, such as amplified slip rates during dynamic rupture, they
-found several effects that would be almost impossible to include
-without resolving both the interseismic deformation and the rapid slip
-during dynamic rupture; the low-rigidity layers reduced the nucleation
-size, amplified slip rates during dynamic rupture, increased the
-recurrent interval, and reduced the amount of aseismic slip
+found several effects that required resolving both the interseismic
+deformation and the rapid slip during dynamic rupture; the
+low-rigidity layers reduced the nucleation size, amplified slip rates
+during dynamic rupture, increased the recurrent interval, and reduced
+the amount of aseismic slip
Collectively, these studies suggest a set of desirable features for
models of the earthquake cycle in order to capture both the slow
deformation associated with interseismic behavior and the rapid
deformation associated with earthquake rupture propagation. These
features include the general capabilities of modeling elasticity with
-elastic, viscoelastic, and viscoelastoplastic deformation, as well as
+elastic, viscoelastic, and viscoelastoplastic rheologies, as well as
slip on faults via either prescribed ruptures or spontaneous ruptures
-controlled by a fault constitutive model. Additionally, a model might
+controlled by a fault constitutive model. Additionally, a model could
also include the coupling of elasticity to fluid and/or heat flow.
Our long-term objective is to develop modeling software with these
@@ -137,40 +138,51 @@ propagation. We plan to seamlessly coupl
propagation. We plan to seamlessly couple these two types of
simulations together to resolve the earthquake cycle.
-Other compelling reasons led us to develop such a code with the
-capability to model both interseismic deformation and earthquake
-rupture propagation. Both types of simulations require the same basic
+Other compelling reasons led us to develop such a code capable of
+modeling both interseismic deformation and earthquake rupture
+propagation. Both types of simulations require the same basic
infrastructure: importing of models from mesh generators, parallel
data structures for finite-elements, bulk constitutive models for
elasticity, fault implementations for prescribed slip and fault
constitutive models, and output of the solution and state variables
-over the domain. The two types of simulations tend to use different
+over the domain. The two types of simulations do tend to use different
boundary conditions, with interseismic deformation usually driven by
Dirichlet (displacement/velocity) or Neumann (traction) boundary
conditions and rupture propagation simulations using absorbing
-boundaries to truncate the domains. However, these features constitute
-a small fraction of the code. The primary different between the two
-types of simulations is the time integration scheme, with an implicit
-scheme used in the quasi-static simulations and an explicit scheme
-used in the dynamic simulations; this also results in using different
-solvers as we will discuss later.
+boundaries to truncate the laterl and bottom boundaries of the
+domains. However, these features constitute a small fraction of the
+code. The primary difference between the two types of simulations is
+the time integration scheme, with an implicit scheme used in the
+quasi-static simulations and an explicit scheme used in the dynamic
+simulations; this also results in using different solvers as we will
+discuss later.
-\brad{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 simulating
+quasi-static and dynamic deformation. Maintaining two
+seperate code bases would require a considerably greater development
+effort.
-Implementing slip on the potentially nonplanar fault surface
-differentiates these types of problems from many other elasticity
-problems. Complexities arise because earthquakes potentially involve
-offset on multiple, intersecting irregularly shaped fault surfaces in
-the interior of the domain. Furthermore, we want the flexibility to
-either prescribe the slip on the fault or have the fault slip evolve
-according to a fault constitutive model that specifies the friction on
-the fault surface. Here, we describe a robust yet flexible method for
-implementing fault slip, its affect on the overall design of PyLith,
+\brad{Insert transition sentence : modeling fault slip is a key
+ feature} Implementing slip on the potentially nonplanar fault
+surface differentiates these types of problems from many other
+elasticity problems. Complexities arise because earthquakes may
+involve offset on multiple, intersecting irregularly shaped fault
+surfaces in the interior of a modeling domain. Furthermore, we want
+the flexibility to either prescribe the slip on the fault or have the
+fault slip evolve according to a fault constitutive model that
+specifies the friction on the fault surface. Here, we describe a
+robust, yet flexible method for implementing fault slip with a domain
+decomposition approach, its affect on the overall design of PyLith,
and verification of its implementation using a few benchmarks.
% ------------------------------------------------------------------
-\section{Numerical Model of Fault Slip}\brad{Rough draft}
+\section{Numerical Model of Fault Slip}
In this section we summarize the formulation of the governing
equations using the finite-element method. We augment the conventional
@@ -185,7 +197,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, \label{eq:faultDisp}
\end{gather}
where $\vec{u}$ is the displacement vector, $\rho$ is the mass
@@ -193,7 +205,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
@@ -250,54 +262,55 @@ 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 trial solution $\vec{u}$, weighting function
-$\vec{\phi}$, and Lagrange multipliers $\vec{l}$ as linear combinations
+We express the weighting function $\vec{\phi}$, trial solution
+$\vec{u}$, and Lagrange multipliers $\vec{l}$ as linear combinations
of basis functions,
\begin{gather}
-\vec{u} = \sum_{m} \vec{u}_m N_m, \\
-\vec{\phi} = \sum_{n} \vec{a}_n N_n, \\
+\vec{\phi} = \sum_{n} \vec{a}_m N_m, \\
+\vec{u} = \sum_{m} \vec{u}_n N_n, \\
\vec{l} = \sum_{p} \vec{l}_p N_p.
\end{gather}
Because the weighting function is zero on $S_u$, the number of basis
functions for the trial solution $\vec{u}$ is generally greater than
the number of basis functions for the weighting function $\vec{\phi}$,
-i.e., $m > n$. The basis functions for the Lagrange multipliers are
+i.e., $n > m$. The basis functions for the Lagrange multipliers are
associated with the fault surface, which is a lower dimension than the
domain, so $p \ll n$ in most cases. If we express the linear
combination of basis functions in terms of a matrix-vector product, we
have
\begin{gather}
-\vec{u} = \tensor{N}_m \vec{u}_m, \\
-\vec{\phi} = \tensor{N}_n \vec{a}_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}_n$, which
+arbitrary, so the integrands must be zero for all $\vec{a}_m$, which
leads to
\begin{gather}
-- \int_{V} \nabla \tensor{N}_n^T \cdot \tensor{\sigma} \, dV
-+ \int_{S_T} \tensor{N}_n^T \vec{T} \, dS
-+ \int_{S_f^{+}} \tensor{N}_n^T \tensor{N}_p \vec{l}_p \, dS
-- \int_{S_f^{-}} \tensor{N}_n^T \tensor{N}_p \vec{l}_p \, dS
-+ \int_{V} \tensor{N}_n^T \vec{f} \, dV
-- \int_{V} \rho \tensor{N}_n^T \tensor{N}_m \frac{\partial^2 \vec{u}_m}{\partial
+ \begin{split}
+- \int_{V} \nabla \tensor{N}_m^T \cdot \tensor{\sigma} \, dV
++ \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}
- - \tensor{N}_{m^+} \vec{u}_{m^+} + \tensor{N}_{m^-} \vec{u}_{m^-}
+ \left( \tensor{R} \cdot \vec{d}
+ - \tensor{N}_{n^+} \vec{u}_{n^+} + \tensor{N}_{n^-} \vec{u}_{n^-}
\right)
\, dS = \vec{0}.
\end{gather}
-We want to solve these equations for the coefficients $\vec{u}_m$ and
+We want to solve these equations for the coefficients $\vec{u}_n$ and
$\vec{l}_p$ subject to $\vec{u} = \vec{u}_0 \text{ on }S_u$.
For nonlinear bulk rheologies it is convenient to work with the
@@ -317,18 +330,22 @@ and the loading conditions. Considering
and the loading conditions. Considering the deformation at time
$\tpdt$,
\begin{gather}
-\label{eqn:quasi-static:residual:domain}
-- \int_{V} \nabla \tensor{N}_n^T \cdot \tensor{\sigma}(\tpdt) \, dV
-+ \int_{S_T} \tensor{N}_n^T \vec{T}(\tpdt) \, dS
-+ \int_{S_f^{+}} \tensor{N}_n^T \tensor{N}_p \vec{l}_p(\tpdt) \, dS
-- \int_{S_f^{-}} \tensor{N}_n^T \tensor{N}_p \vec{l}_p(\tpdt) \, dS
-+ \int_{V} \tensor{N}_n^T \vec{f}(\tpdt) \, dV
-=\vec{0}, \\
+ \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 \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 \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}_{m^+} \vec{u}_{m^+}(\tpdt) + \tensor{N}_{m^-} \vec{u}_{m^-}(\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}
@@ -337,27 +354,28 @@ solution. We solve these equations usin
solution. We solve these equations using the Portable, Extensible
Toolkit for Scientific Computation (PETSc), which provides a suite of
tools for solving linear systems of algebraic equations with parallel
-processing \cite{petsc-user-ref,petsc-web-page}. In solving the system, we compute the
-residual (i.e., $\vec{r} = \vec{b} - \tensor{A} \vec{u}$ and the
-Jacobian of the system ($A$). In our case the solution is $\vec{u} =
-\left( \begin{smallmatrix} \vec{u}_m \\ \vec{l}_m \end{smallmatrix} \right)$,
-and the residual is simply the left hand sides of
+processing \cite{PETSc:web:page,PETSc:manual,PETSC:efficient}. In
+solving the system, we compute the residual (i.e., $\vec{r} = \vec{b}
+- \tensor{A} \vec{u}$ and the Jacobian of the system ($A$). In our
+case the solution is $\vec{u} = \left( \begin{smallmatrix} \vec{u}_n
+ \\ \vec{l}_n \end{smallmatrix} \right)$, and the residual is
+simply the left hand sides of
equations~(\ref{eqn:quasi-static:residual:domain})
and~(\ref{eqn:quasi-static:residual:fault}). That is,
\begin{equation}
\vec{r} = \left( \begin{array}{c}
\begin{aligned}
- - \int_{V} \nabla \tensor{N}_n^T \cdot \tensor{\sigma}(\tpdt) \, dV
- + \int_{S_T} \tensor{N}_n^T \vec{T}(\tpdt) \, dS
- + \int_{S_f^{+}} \tensor{N}_n^T \tensor{N}_p \vec{l}_p(\tpdt) \, dS
- &- \int_{S_f^{-}} \tensor{N}_n^T \tensor{N}_p \vec{l}_p(\tpdt) \, dS
+ - \int_{V} \nabla \tensor{N}_m^T \cdot \tensor{\sigma}(\tpdt) \, dV
+ + \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}_n^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}_{m^+} \vec{u}_{m^+}(\tpdt) + \tensor{N}_{m^-} \vec{u}_{m^-}(\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).
@@ -381,14 +399,16 @@ displacement vector as a linear combinat
displacement vector as a linear combination of basis functions, after
some algebra we find this portion of the Jacobian is
\begin{equation}
+ \label{eqn:jacobian:implicit:stiffness}
\tensor{K} = \frac{1}{4} \int_v
- (\nabla^T + \nabla) \tensor{N}_n^T \cdot
- \tensorfour{C} \cdot (\nabla + \nabla^T) \tensor{N}_m \, dV.
+ (\nabla^T + \nabla) \tensor{N}_m^T \cdot
+ \tensorfour{C} \cdot (\nabla + \nabla^T) \tensor{N}_n \, dV.
\end{equation}
Following a similar procedure, we find the portion of the Jacobian
associated with equation~(\ref{eqn:quasi-static:residual:fault}) is
\begin{equation}
- \tensor{L} = \int_{S_f} \tensor{N}_p (-\tensor{N}_{m^+}+\tensor{N}_{m^-}) \, dS.
+ \label{eqn:jacobian:constraint}
+ \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}\label{eqn:saddle-point}
@@ -411,19 +431,19 @@ inertial term:
\begin{equation}
\vec{r} = \left( \begin{array}{c}
\begin{aligned}
- - \int_{V} \nabla \tensor{N}_n^T \cdot \tensor{\sigma}(\tpdt) \, dV
- + \int_{S_T} \tensor{N}_n^T \vec{T}(\tpdt) \, dS
- &+ \int_{S_f^{+}} \tensor{N}_n^T \tensor{N}_p \vec{l}_p(\tpdt) \, dS
- - \int_{S_f^{-}} \tensor{N}_n^T \tensor{N}_p \vec{l}_p(\tpdt) \, dS
+ - \int_{V} \nabla \tensor{N}_m^T \cdot \tensor{\sigma}(\tpdt) \, dV
+ + \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}_n^T \vec{f}(\tpdt) \, dV
- - \int_{V} \rho \tensor{N}_n^T \tensor{N}_m
- \frac{\partial^2 \vec{u}_m}{\partial t^2} \, dV
+ &+ \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}_{m^+} \vec{u}_{m^+}(\tpdt) + \tensor{N}_{m^-} \vec{u}_{m^-}(\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).
@@ -447,16 +467,17 @@ difference scheme, so that the accelerat
\end{gather}
Expanding the inertial term yields
\begin{equation}
- - \int_{V} \rho \tensor{N}_n^T \tensor{N}_m \frac{\partial^2 \vec{u}_m}{\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}_n^T \tensor{N}_m \left(
- d\vec{u}_m(t) - \vec{u}_m(t) + \vec{u}_m(\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}_n^T \tensor{N}_m \, dV.
+ \frac{1}{\Delta t^2} \int_{V} \rho \tensor{N}_m^T\ \cdot \tensor{N}_n \, dV.
\end{equation}
@@ -490,8 +511,8 @@ the change in the tractions on the fault
the change in the tractions on the fault surfaces corresponding to the
slip. PyLith includes a variety of slip-time histories, including a
step function, constant slip rate, the integral of Brune's far-field
-time function, a sine-cosine function developed by
-\citeN{Liu:etal:200?}, and a user-defined time history. These are
+time function \cite{Brune:1970}, a sine-cosine function developed by
+\citeN{Liu:etal:2006}, and a user-defined time history. These are
discussed in detail in the PyLith manual
\cite{PyLith:manual:1.6.2}. PyLith allows multiple earthquake ruptures
to be superimposed on each other, thereby permitting rather complex
@@ -527,55 +548,61 @@ of the general form of a linear system o
\vec{u} = \vec{b}$), our subset of equations has the form
\begin{equation}
\begin{pmatrix}
- \tensor{K}_{m^+m^+} & 0 & \tensor{L}_p \\
- 0 & \tensor{K}_{m^-m^-} & -\tensor{L}_p \\
+ \tensor{K}_{n^+n^+} & 0 & \tensor{L}_p^T \\
+ 0 & \tensor{K}_{n^-n^-} & -\tensor{L}_p^T \\
\tensor{L}_p & -\tensor{L}_p & 0
\end{pmatrix}
\begin{pmatrix}
- \vec{u}_{m^+} \\
- \vec{u}_{m^-} \\
+ \vec{u}_{n^+} \\
+ \vec{u}_{n^-} \\
\vec{l}_p \\
\end{pmatrix}
=
\begin{pmatrix}
- \vec{b}_{m^+} \\
- \vec{b}_{m^-} \\
+ \vec{b}_{n^+} \\
+ \vec{b}_{n^-} \\
\vec{b}_p \\
\end{pmatrix},
\end{equation}
-where $m^+$ and $m^-$ refer to the degrees of freedom associated with
+where $n^+$ and $n^-$ refer to the degrees of freedom associated with
the positive and negative sides of the fault,
-respectively. Furthermore, we can ignore the terms $\vec{b}_{m^+}$ and
-$\vec{b}_{m^-}$ because they do not change as we change the Lagrange
+respectively. Furthermore, we can ignore the terms $\vec{b}_{n^+}$ and
+$\vec{b}_{n^-}$ because they do not change as we change the Lagrange
multipliers or fault slip. This means we can solve the following
equations to estimate the change in fault slip $\partial \vec{d}$
corresponding to a perturbation in the Lagrange multipliers $\partial
\vec{l}_p$:
\begin{gather}
- \tensor{K}_{m^+m^+} \vec{u}_{m^+} = - \tensor{L}_p \partial \vec{l}_p, \\
- \tensor{K}_{m^-m^-} \vec{u}_{m^-} = - \tensor{L}_p \partial \vec{l}_p, \\
- \partial \vec{d} = \tensor{R} \left( \vec{u}_{m^+} - \vec{u}_{m^-} \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
if the deformation associated with the fault slip is confined close to
-the fault (which occurs in most simulations because the fault slip is
-limited to the middle of the domain), and poor if the deformation
-associated with fault slip extends to the boundaries (which occurs in
-races cases in which the domain has a through-going fault).
+the fault (which occurs in most simulations because the deformation
+from fault slip is limited to the middle of the domain), and poor if
+the deformation associated with fault slip extends to the boundaries
+(which occurs in the rare cases in which the domain has a through-going
+fault).
PyLith includes several commonly used fault constitutive models, all
of which specify the shear traction on the fault $T_f$ as a function
of the cohesive stress $T_c$, coefficient of friction, $\mu_f$, and
-normal traction $T_n$,
+normal traction $T_m$,
\begin{equation}
- T_f = T_c - \mu_f T_n.
+ 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 Dieterich-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}
@@ -699,12 +726,184 @@ PETSc. Figure~(\ref{fig:solver-options})
\matt{Table of convergence for a fault example with ASM, block tridiagonal, block tridiagonal+custom preconditioner, LSC}
-\subsection{Dynamic-static Simulations}\brad{Rough draft}
+\subsection{Dynamic Simulations}
- \begin{itemize}
- \item Lumped Jacobian
- \item Schur complement
- \end{itemize}
+In the dynamic simulations the the Courant-Friderichs-Lewy condition
+controls the stability of the time integration. In most dynamic
+problems this dictates a relatively small time step so that a typical
+simulation involves tens of thousands of time steps. Hence, we want a
+very efficient solver in order to run dynamic simulations in a
+reasonable amount of time.
+
+% Lumped Jacobian
+
+The Jacobian for our system of equations involves two terms: the
+inertial term given by equation~(\ref{eqn:jacobian:explicit:inertia})
+and the fault slip constraint term given by
+equation~(\ref{eqn:jacobian:constraint}). Using conventional
+finite-element basis functions in these integrations results in a
+sparse matrix with off-diagonal terms. Although we can use the same
+solvers as we do for quasi-static simulations to find the solution,
+eliminating the off-diagonal terms so that the Jacobian is diagonal,
+permits use of a much faster solver. With a diagonal Jacobian the
+number of operations required for the solve is proportional to the
+number of degrees of freedom, and the memory requirements are greatly
+reduced by storing the diagonal of the matrix as a vector rather
+than as a sparse matrix. However, the block structure of our Jacobian
+matrix, with the fault slip constraints occupying off-diagonal blocks,
+requires a two step approach in order to solve the linear system
+of equations without forming a sparse matrix.
+
+First, we eliminate the off-diagonal entries in each block of the
+matrix during the finite-element integrations. The current best
+available option for eliminating the off-diagonal terms formed during
+the integration of the inertial term focuses on choosing a set of
+orthogonal basis functions, such as the Legendre polynomials with
+Gauss-Lobatto-Legendre quadrature points
+\cite{Komatitsch:Vilotte:1998}. This discretization (often called the
+spectral element method) naturally produces a diagonal block for each
+finite-element cell without introducing any additional
+approximations. Because the fault slip constraint term also involves
+integration of the products of the basis functions over
+lower-dimension cells, orthogonal basis functions also produce a
+diagonal block for this integration.
+
+In contrast, traditional finite-element approaches do introduce
+additional approximations when constructing a diagonal
+approximation. In PyLith we employ one of these traditional
+approaches, because it produces good approximations for many different
+choices of basis functions and quadrature points. For each
+finite-element cell, we construct a diagonal approximation of the
+integral such that the action on rigid body motion is the same for the
+diagonal approximation of the integral as it is for the original
+integral,
+\begin{equation}
+ \tensor{A} \cdot \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 \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
+$\Omega$ may be the domain volume (as in the case of the inertial
+term) or a boundary (as in the case of the fault slip constraint
+term).
+
+The errors associated with this approximation are small as
+long as the deformation occurs at length scales significantly larger
+than the discretization size, which is consistent with resolving
+seismic wave propagation accurately. Furthermore, in contrast to other
+approaches that choose basis functions or quadrature points that
+affect the accuracy of all of the finite-element integrations, such as
+choosing quadrature points coincident with the vertices of a cell,
+this approach only affects the accuracy of the terms involved in the
+Jacobian. For consistency in the formulation of the system of
+equations, these approximations are also applied to the inertial term
+and fault slip constraint term when computing the residual.
+
+
+% Schur complement
+
+Second, we leverage the structure of the off-diagonal blocks
+associated with the fault slip constraint in solving the system of
+equations via a Schur's complement algorithm. We compute an initial
+residual assuming the increment in the solution is zero (i.e.,
+$d\vec{u}_n = \vec{0}$ and $d\vec{l}_p = \vec{0}$,
+\begin{equation}
+ \vec{r}^* = \begin{pmatrix} \vec{r}_n^* \\ \vec{r}_p^* \end{pmatrix} =
+ \begin{pmatrix} \vec{b}_n \\ \vec{b}_p \end{pmatrix}
+ - \begin{pmatrix}
+ \tensor{K} & \tensor{L}^T \\ \tensor{L} & 0
+ \end{pmatrix}
+ \begin{pmatrix} \vec{u}_n \\ \vec{l}_n \end{pmatrix}.
+\end{equation}
+ We compute a corresponding initial solution to the system of equations
+$\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} \cdot \vec{r}_n,
+\end{equation}
+taking advantage of the fact that we construct $\tensor{K}$ so that it
+is diagonal.
+
+We next compute the increment in the Lagrange multipliers in order to
+correct this initial solution so that the true residual is zero.
+Making use of the initial residual, the expression for the true
+residual is
+\begin{equation}
+ \label{eqn:lumped:jacobian:residual}
+ \vec{r} = \begin{pmatrix} \vec{r}_n \\ \vec{r}_p \end{pmatrix} =
+ \begin{pmatrix} \vec{r}_n^* \\ \vec{r}_p^* \end{pmatrix}
+ - \begin{pmatrix}
+ \tensor{K} & \tensor{L}^T \\ \tensor{L} & 0
+ \end{pmatrix}
+ \begin{pmatrix} d\vec{u}_n \\ d\vec{l}_n \end{pmatrix}.
+\end{equation}
+Solving the first row of equation~(\ref{eqn:lumped:jacobian:residual})
+for the increment in the solution and accounting for the structure of
+$\tensor{L}$ as we write the expressions for degrees of freedom on
+each side of the fault, we have
+\begin{gather}
+ d\vec{u}_{n^+} =
+ 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} \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 \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 \cdot
+ \left( \tensor{K}_{n^+n^+}^{-1} + \tensor{K}_{n^-n^-}^{-1} \right) \cdot
+ \tensor{L}^T_p,
+\end{equation}
+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} \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} \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
+the degrees of freedom attached to the fault and the Lagrange
+multipliers.
+
+% Spontaneous rupture model and lumped Jacobian
+
+We also leverage the elimination of off-diagonal entries from the
+blocks of the Jacobian in dynamic simulations when updating the slip
+in spontaneous rupture models. Because $\tensor{K}$ is diagonal in
+this case, the expression for the change in slip for a perturbation in
+the Lagrange multipliers
+(equation~(\ref{eqn:spontaneous:rupture:slip:update})) is exact and
+simplifies to
+\begin{equation}
+ \partial \vec{d} = - \tensor{R}^{-1} \cdot
+ \left( \tensor{K}_{n^+n^+}^{-1} + \tensor{K}_{n^-n^-}^{-1} \right)
+ \cdot \tensor{L}_p^T \cdot \partial \vec{l}_p.
+\end{equation}
+
% ------------------------------------------------------------------
@@ -721,6 +920,12 @@ PETSc. Figure~(\ref{fig:solver-options})
% ------------------------------------------------------------------
\begin{acknowledgments}
+
+\begin{itemize}
+ \item Funding from SCEC, CIG, GNS Science, USGS.
+ \item SCEC contribution number.
+ \end{itemize}
+
\end{acknowledgments}
MGK acknowledges partial support from NSF grant EAR-0949446.
diff -r 4dde491c41e8 -r 26022948ab53 pylith.bib
--- a/pylith.bib Thu Sep 01 13:17:57 2011 -0500
+++ b/pylith.bib Thu Sep 01 20:33:15 2011 -0500
@@ -31,28 +31,3 @@
month = may,
}
- at Article{Liu:etal:2006,
- author = {Liu, P. and Archuleta, R.~J. and Hartzell, S.~H.},
- title = {Prediction of Broadband Ground-Motion Time
- Histories: {Hybrid} Low/High- Frequency Method with
- Correlated Random Source Parameters},
- journal = BSSA,
- year = {2006},
- volume = {96},
- number = {6},
- pages = {2118--2130},
- month = dec,
- doi = {10.1785/0120060036},
-}
-
- at article{Brune:1970,
- author = {Brune, James N.},
- title = {Tectonic stress and spectra of seismic shear waves
- from earthquakes},
- journal = JGR,
- year = 1970,
- volume = 75,
- month = sep # {~10},
- pages = {4997--5009},
-}
-
diff -r 4dde491c41e8 -r 26022948ab53 references.bib
--- a/references.bib Thu Sep 01 13:17:57 2011 -0500
+++ b/references.bib Thu Sep 01 20:33:15 2011 -0500
@@ -600,7 +600,7 @@
± 10\%. Events M ≥ 6 have as much as a 60\%
probability of recurrence within 5 years due to the
clustering of small earthquakes in foreshocks and
- aftershocks. This probability drops to less than 15%
+ aftershocks. This probability drops to less than 15\%
for M ≥ 7 events. Increasing gap time generally
increases conditional probability of earthquake
occurrence, but the effect is weak. For the MAT, the
@@ -621,6 +621,20 @@
year = {1996}
}
+ at Article{Liu:etal:2006,
+ author = {Liu, P. and Archuleta, R.~J. and Hartzell, S.~H.},
+ title = {Prediction of Broadband Ground-Motion Time
+ Histories: {Hybrid} Low/High- Frequency Method with
+ Correlated Random Source Parameters},
+ journal = BSSA,
+ year = {2006},
+ volume = {96},
+ number = {6},
+ pages = {2118--2130},
+ month = dec,
+ doi = {10.1785/0120060036},
+}
+
@techreport{petsc-user-ref,
author = {Satish Balay and Jed Brown and Kris Buschelman and Victor Eijkhout and William~D. Gropp and Dinesh Kaushik and Matthew~G. Knepley and
Lois Curfman McInnes and Barry~F. Smith and Hong Zhang},
@@ -631,12 +645,84 @@
url = {http://www.mcs.anl.gov/petsc}
}
+ at article{Brune:1970,
+ author = {Brune, James N.},
+ title = {Tectonic stress and spectra of seismic shear waves
+ from earthquakes},
+ journal = JGR,
+ year = 1970,
+ volume = 75,
+ month = sep # {~10},
+ pages = {4997--5009},
+}
+
@misc{petsc-web-page,
author = {Satish Balay and Jed Brown and Kris Buschelman and Victor Eijkhout and William~D. Gropp and Dinesh Kaushik and Matthew~G. Knepley and Lois Curfman McInnes and Barry~F. Smith and Hong Zhang},
title = {{PETS}c {W}eb page},
url = {http://www.mcs.anl.gov/petsc},
howpublished = {\url{http://www.mcs.anl.gov/petsc}},
year = {2011}
+}
+
+ at article{Komatitsch:Vilotte:1998,
+ author = {Komatitsch, D. and Vilotte, J.~P.},
+ title = {The spectral element method: {An} efficient tool to
+ simulate the seismic response of 2D and 3D
+ geological structures},
+ journal = BSSA,
+ year = 1998,
+ volume = 88,
+ number = 2,
+ month = apr,
+ pages = {368--392},
+ abstract = {We present the spectral element method to simulate
+ elastic-wave propagation in realistic geological
+ structures involving complieated free-surface
+ topography and material interfaces for two- and
+ three-dimensional geometries. The spectral element
+ method introduced here is a high-order variational
+ method for the spatial approximation of elastic-wave
+ equations. The mass matrix is diagonal by
+ construction in this method, which drastically
+ reduces the computational cost and allows an
+ efficient parallel implementation. Absorbing
+ boundary conditions are introduced in variational
+ form to simulate unbounded physical domains. The
+ time discretization is based on an energy-momentum
+ conserving scheme that can be put into a classical
+ explicit-implicit predictor/multi-corrector
+ format. Long-term energy conservation and stability
+ properties are illustrated as well as the efficiency
+ of the absorbing conditions. The associated Courant
+ condition behaves as {Delta}tC < O (nel–1/ndN–2),
+ with nel the number of elements, nd the spatial
+ dimension, and N the polynomial order. In practice,
+ a spatial sampling of approximately 5 points per
+ wavelength is found to be very accurate when working
+ with a polynomial degree of N = 8. The accuracy of
+ the method is shown by comparing the spectral
+ element solution to analytical solutions of the
+ classical two-dimensional (2D) problems of Lamb and
+ Garvin. The flexibility of the method is then
+ illustrated by studying more realistic 2D models
+ involving realistic geometries and complex
+ free-boundary conditions. Very accurate modeling of
+ Rayleigh-wave propagation, surface diffraction, and
+ Rayleigh-to-body-wave mode conversion associated
+ with the free-surface curvature are obtained at low
+ computational cost. The method is shown to provide
+ an efficient tool to study the diffraction of
+ elastic waves by three-dimensional (3D) surface
+ topographies and the associated local effects on
+ strong ground motion. Complex amplification
+ patterns, both in space and time, are shown to occur
+ even for a gentle hill topography. Extension to a
+ heterogeneous hill structure is considered. The
+ efficient implementation on parallel distributed
+ memory architectures will allow to perform real-time
+ visualization and interactive physical
+ investigations of 3D amplification phenomena for
+ seismic risk assessment.},
}
@article{KnepleyKarpeev09,
@@ -649,3 +735,60 @@
year = {2009},
doi = {10.3233/SPR-2009-0249}
}
+
+ at Book{Zienkiewicz:Taylor:2005,
+ author = {Ziekiewicz, O.~C. and Taylor, R.~L. and Zhu, J.~Z.},
+ title = {The Finite-Element Method: Its Basis and Fundamentals},
+ publisher = {Butterworth-Heinemann},
+ year = {2005},
+ edition = {6th},
+}
+
+ at Misc{PETSc:web:page,
+ Author = "Balay, S. and Brown, J. and Buschelman, K. and Gropp,
+ W.~D. and Kaushik, D. and Knepley, M.~G. and
+ McInnes, L.~C. and Smith, B.~F. and Zhang, H.",
+ Title = "{PETSc} {Web} page",
+ Note = "http://www.mcs.anl.gov/petsc",
+ Year = 2011,
+}
+
+ at TechReport{PETSc:manual,
+ Author = "Balay, S. and Brown, J. and Buschelman, K. and Gropp,
+ W.~D. and Kaushik, D. and Knepley, M.~G. and
+ McInnes, L.~C. and Smith, B.~F. and Zhang, H.",
+ Title = "{PETSc} Users Manual",
+ Number = "ANL-95/11 - Revision 3.1",
+ Institution = "Argonne National Laboratory",
+ Year = 2010,
+}
+
+ at InProceedings{PETSc:efficient,
+ Author = "Balay, S. and Gropp, W.~D. and McInnes, L.~C. and
+ Smith, B.~F. ",
+ Title = "Efficient Management of Parallelism in Object
+ Oriented Numerical Software Libraries",
+ Booktitle = "Modern Software Tools in Scientific Computing",
+ Editor = "E. Arge and A. M. Bruaset and H. P. Langtangen",
+ Pages = "163--202",
+ Publisher = "Birkh{\"{a}}user Press",
+ Year = 1997,
+}
+
+ at Manual{PyLith:manual:1.6.2,
+ title = {PyLith User Manual, Version 1.6.2},
+ author = {Aagaard, B. and Kientz, S. and Knepley, M. and
+ Somala, S. and Strand, L. and Williams, C.},
+ organization = {Computational Infrastructure for Geodynamics (CIG)},
+ address = {University of California, Davis},
+ 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