[cig-commits] commit: Edits based on internal reviews.
Mercurial
hg at geodynamics.org
Tue Aug 28 12:03:17 PDT 2012
changeset: 140:6b10724d5da2
parent: 138:91049a93b5e2
user: Brad Aagaard <baagaard at usgs.gov>
date: Tue Aug 28 12:02:36 2012 -0700
files: faultRup.tex references.bib
description:
Edits based on internal reviews.
diff -r 91049a93b5e2 -r 6b10724d5da2 faultRup.tex
--- a/faultRup.tex Tue Aug 14 13:29:49 2012 -0700
+++ b/faultRup.tex Tue Aug 28 12:02:36 2012 -0700
@@ -75,14 +75,10 @@
schemes, and bulk and fault rheologies. We have developed a custom
preconditioner for the Lagrange multiplier portion of the system of
equations that provides excellent scalability with problem size
- compared to conventional Additive Schwarz methods. We demonstrate
+ compared to conventional additive Schwarz methods. We demonstrate
application of this approach using benchmarks for both quasi-static
viscoelastic deformation and dynamic spontaneous rupture propagation
- that verify the numerical implementation. Future work will focus on
- linking the quasi-static and dynamic simulations together to capture
- both the slow strain accumulation and post-seismic relaxation at
- long time scales and the dynamic rupture propagation and radiation
- of seismic waves at short time scales.
+ that verify the numerical implementation in PyLith.
\end{abstract}
% ------------------------------------------------------------------
@@ -91,19 +87,19 @@
% ------------------------------------------------------------------
\section{Introduction}
-Understanding the earthquake cycle, from slow deformation associated
-with interseismic behavior to rapid deformation associated with
-earthquake rupture, spans spatial scales ranging from tiny fractions
-of a meter associated with the size of contact asperities on faults
-and individual grains to thousands of kilometers associated with plate
+The earthquake cycle, from slow deformation associated with
+interseismic behavior to rapid deformation associated with earthquake
+rupture, spans spatial scales ranging from fractions of a meter
+associated with the size of contact asperities on faults and
+individual grains to thousands of kilometers associated with plate
boundaries. Similarly, temporal scales range from fractions of a
second associated with slip at a point during earthquake rupture to
thousands of years of strain accumulation between earthquakes. The
-complexity of dealing with the many physical processes operating over
-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 often justifies this narrow focus.
+complexity of the many physical processes operating over this vast
+range of scales 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 often justifies
+this narrow focus.
Researchers have recognized for some time, though, that interseismic
deformation and fault interactions influence earthquake rupture
@@ -118,14 +114,14 @@ rupture propagation often approximate th
rupture propagation often approximate the loading of the crust via
simplistic assumptions about the stress field at the beginning of a
rupture
-\citep{Aagaard:etal:BSSA:2001,Peyrat:etal:2001,Oglesby:Day:2001,Dunham:Archuleta:2004}. Earthquake
-simulators, which attempt to model multiple earthquake cycles,
+\citep{Mikumo:etal:1998,Harris:Day:1999,Aagaard:etal:BSSA:2001,Peyrat:etal:2001,Oglesby:Day:2001,Dunham:Archuleta:2004}. Numerical
+seismicity models that attempt to model multiple earthquake cycles,
generally simplify not only the fault loading and rupture propagation
but also the physical properties in order to make the calculations
tractable
-\citep{Ward:1992,Rundle:etal:2006,Pollitz:Schwartz:2008,Dieterich:Richards-Dinger:2010}.
+\citep{Ward:1992,Robinson:Benites:1995,Rundle:etal:2006,Pollitz:Schwartz:2008,Dieterich:Richards-Dinger:2010}.
-Several recent dynamic spontaneous rupture modeling studies have
+Some dynamic spontaneous rupture modeling studies have
attempted to examine a broader space-time window in order to remove
simplifying assumptions and more accurately capture the complex
interactions over the earthquake cycle. For example,
@@ -169,50 +165,14 @@ controlled by a fault constitutive model
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
-capabilities in order to simulate the earthquake cycle, resolving the
-deformation across as wide of a range of temporal and spatial scales
-as possible to avoid simplifications that affect earthquake cycle
-behavior. Our current focus is developing a modeling code, PyLith,
-that supports both quasi-static simulations of interseismic and
-coseismic deformation and dynamic simulations of earthquake rupture
-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 capable of
-modeling both interseismic deformation and earthquake rupture
-propagation. Both types of simulations require the same basic
-infrastructure: importing 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 generally 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 lateral 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.
-
-Additional motivation for developing PyLith arose from the geophysics
-community as part of the Computational Infrastructure for Geodynamics
-(CIG) project (\url{http://www.geodynamics.org}). 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
-separate code bases would require a considerably greater development
-effort.
-
With the goal of modeling the entire earthquake cycle with as few
-simplifications as possible, much of the code development has focused
-on modeling fault slip. Implementing slip on the potentially nonplanar
-fault surface differentiates these types of problems from many other
+simplifications as possible, much of our work in developing PyLith has
+focused on modeling fault slip with application to quasi-static
+simulations of interseismic and coseismic deformation and dynamic
+simulations of earthquake rupture propagation. We plan to seamlessly
+couple these two types of simulations together to resolve the
+earthquake cycle. 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
@@ -417,29 +377,31 @@ system, we compute the residual (i.e., $
($\bm{A}$). In our case the solution is $\bm{u} =
\left( \begin{smallmatrix} \bm{u}_n \\
\bm{l}_n \end{smallmatrix} \right)$, and the residual is
-simply the left hand sides of
+simply the left sides of
equations~(\ref{eqn:quasi-static:residual:elasticity})
and~(\ref{eqn:quasi-static:residual:fault}).
-The Jacobian of the system is the action (operation) on the increment
-in the solution. To find the portion of the Jacobian associated with
+The Jacobian of the system, $\bm{A}$, is the action (operation) that
+we apply to the increment of the solution, $\bm{du}$. To find the
+portion of the Jacobian associated with
equation~(\ref{eqn:quasi-static:residual:elasticity}), we let
-$\bm{\sigma}(t+\Delta t) = \bm{\sigma}(t) +
-\bm{d\sigma}(t)$. The action on the increment of the solution is
-associated with the increment in the stress tensor
-$\bm{d\sigma}(t)$. We approximate the increment in the stress
-tensor using linear elasticity and infinitesimal strains,
+$\bm{\sigma}(t+\Delta t) = \bm{\sigma}(t) + \bm{d\sigma}(t)$. The
+action on the increment of the solution is associated with the
+increment in the stress tensor $\bm{d\sigma}(t)$. We approximate the
+increment in the stress tensor using linear elasticity and
+infinitesimal strains,
\begin{equation}
\bm{d\sigma}(t) = \frac{1}{2} \bm{C}(t) \cdot (\nabla + \nabla^T)
\bm{u}(t),
\end{equation}
-where $\bm{C}$ is the fourth order tensor of elastic
-constants. Substituting into the first term in
-equation~(\ref{eqn:quasi-static:residual:elasticity}) and expressing the
-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}
+where $\bm{C}$ is the fourth order tensor of elastic constants. For
+bulk consitutive models with a linear response, $\bm{C}$ is constant
+in time. For other constitutive models we form $\bm{C}(t)$ from the
+current solution and state variables. Substituting into the first term
+in equation~(\ref{eqn:quasi-static:residual:elasticity}) and
+expressing the displacement vector as a linear combination of basis
+functions, we find this portion of the Jacobian is
+\begin{equation}\label{eqn:jacobian:implicit:stiffness}
\bm{K} = \frac{1}{4} \int_V
(\nabla^T + \nabla) \bm{N}_m^T \cdot
\bm{C} \cdot (\nabla + \nabla^T) \bm{N}_n \, dV.
@@ -448,8 +410,7 @@ finite-element formulations. Following a
finite-element formulations. Following a similar procedure, we find
the portion of the Jacobian associated with the constraints,
equation~(\ref{eqn:quasi-static:residual:fault}), is
-\begin{equation}
- \label{eqn:jacobian:constraint}
+\begin{equation}\label{eqn:jacobian:constraint}
\bm{L} = \int_{S_f} \bm{N}_p^T \cdot (\bm{N}_{n^+} - \bm{N}_{n^-}) \, dS.
\end{equation}
Thus, the Jacobian of the entire system has the form,
@@ -465,8 +426,7 @@ we compute the terms for the positive si
we compute the terms for the positive side of the fault and assemble
the terms into the appropriate DOF for both sides of
the fault. Hence, we compute
-\begin{equation}
- \label{eqn:jacobian:constraint:code}
+\begin{equation}\label{eqn:jacobian:constraint:code}
\bm{L_p} = \int_{S_f} \bm{N}_p^T \cdot \bm{N}_{n^+} \, dS,
\end{equation}
with the Jacobian of the entire system taking the form,
@@ -489,13 +449,15 @@ associated with the Lagrange multipliers
\subsection{Dynamic Simulations}
In dynamic simulations we include the inertial term in order to
-resolve the propagation of seismic waves. The integral equation for
-the fault slip constraint remains unchanged, so the corresponding
-portions of the Jacobian ($\bm{L}$) and residual ($\bm{r}_p$)
-are exactly the same as in the quasi-static simulations. Including the
-inertial term in equation~(\ref{eqn:quasi-static:residual:elasticity}) yields
-\begin{equation}
- \label{eqn:dynamic:residual:elasticity}
+resolve the propagation of seismic waves. The general form of the
+system Jacobian remains the same as in quasi-static
+simulations given in equation~(\ref{eqn:saddle:point}). The integral
+equation for the fault slip constraint remains unchanged, so the
+corresponding portions of the Jacobian ($\bm{L}$) and residual
+($\bm{r}_p$) are also exactly the same as in the quasi-static
+simulations. Including the inertial term in
+equation~(\ref{eqn:quasi-static:residual:elasticity}) yields
+\begin{equation}\label{eqn:dynamic:residual:elasticity}
\begin{split}
- \int_{V} \nabla \bm{N}_m^T \cdot \bm{\sigma}(t+\Delta t) \, dV
+ \int_{S_T} \bm{N}_m^T \cdot \bm{T}(t+\Delta t) \, dS \\
@@ -508,6 +470,7 @@ inertial term in equation~(\ref{eqn:quas
=\bm{0}.
\end{split}
\end{equation}
+
We find the upper portion of the Jacobian of the system by considering
the action on the increment in the solution, just as we did for the
quasi-static simulations. In this case we associate the increment in
@@ -563,25 +526,6 @@ 0.1--1.0.
0.1--1.0.
% ------------------------------------------------------------------
-\subsection{Leveraging Common Components}
-
-Many of the advantages of designing PyLith to handle both quasi-static
-and dynamic simulations should now be apparent. There are only minor
-differences between the formulations for the quasi-static and dynamic
-simulations. As a result both types of simulations can use the same
-data structures and many of the same routines for assembling the
-system of equations. In particular, the integrals for the bulk
-rheologies and fault behavior are identical. This means that all of
-the various elastic, viscoelastic, and viscoelastoplastic bulk
-constitutive models and fault rupture behavior via slip-time functions
-and fault constitutive models need only a single implementation in
-order to be used for both quasi-static and dynamic simulations. This
-significantly reduces development time and simplifies maintenance
-compared to maintaining different codes for quasi-static and dynamic
-simulations.
-
-
-% ------------------------------------------------------------------
\subsection{Nondimensionalization}
The domain decomposition approach for implementing fault slip with
@@ -594,9 +538,9 @@ magnitude. This results in a rather ill-
magnitude. This results in a rather ill-conditioned system. We avoid
this ill-conditioning by nondimensionalizing all of the quantities
involved in the problem based upon user-specified scales
-\citep{PyLith:manual:1.7.1}, facilitating forming well-conditioned
-systems of equations for problems across a wide range of spatial and
-temporal scales.
+\citep{PyLith:manual:1.7.1}, facilitating the formation of
+well-conditioned systems of equations for problems across a wide range
+of spatial and temporal scales.
% ------------------------------------------------------------------
@@ -624,17 +568,17 @@ times, thereby permitting complex slip b
% ------------------------------------------------------------------
\subsection{Spontaneous Fault Rupture}
-In a spontaneous fault rupture a constitutive model controls the
-tractions on the fault surface. The fault slip evolves based on the
-fault tractions as driven by the initial conditions, boundary
-conditions and deformation. In our formulation of fault slip, we
-specify the fault slip, and the fault tractions (Lagrange multipliers)
-are part of the solution (unknowns). The fault constitutive model
-places bounds on the Lagrange multipliers and the system of equations
-is nonlinear-- when a location on the fault is slipping, we must solve
-for both the fault slip (which is known in the prescribed ruptures)
-and the Lagrange multipliers to find values consistent with the fault
-constitutive model.
+In contrast to prescribed fault rupture, in spontaneous fault rupture
+a constitutive model controls the tractions on the fault surface. The
+fault slip evolves based on the fault tractions as driven by the
+initial conditions, boundary conditions and deformation. In our
+formulation of fault slip, fault slip is assumed to be known and the
+fault tractions (Lagrange multipliers) are part of the solution
+(unknowns). The fault constitutive model places bounds on the Lagrange
+multipliers and the system of equations is nonlinear-- when a location
+on the fault is slipping, we must solve for both the fault slip (which
+is known in the prescribed ruptures) and the Lagrange multipliers to
+find values consistent with the fault constitutive model.
At each time step, we first assume the increment in the fault slip is
zero, so that the Lagrange multipliers correspond to the fault
@@ -718,8 +662,8 @@ where $f(\bm{d})$ corresponds to the fau
where $f(\bm{d})$ corresponds to the fault constitutive model and
$\|x\|_2$ denotes the L$^2$-norm of $x$. Performing this search in
logarithmic space rather than linear space greatly accelerates the
-convergence in rate-state fault constitutive models in which the
-coefficient of friction depends on the logarithm of the slip rate.
+convergence in constitutive models in which the coefficient of
+friction depends on the logarithm of the slip rate.
PyLith includes several commonly used fault constitutive models, all
of which specify the shear traction on the fault $T_f$ as a function
@@ -735,33 +679,36 @@ compression, we prevent interpenetration
compression, we prevent interpenetration, and when the fault is under
tension, the fault opens ($d_n > 0$) and the fault traction vector is
zero. The fault constitutive models include static friction, linear
-slip-weakening, linear time-weakening, and Dieterich-Ruina rate-state
-friction with an aging law (see the PyLith manual
-\citep{PyLith:manual:1.7.1} for details).
+slip-weakening \citep{Ida:1972}, linear time-weakening
+\citep{Andrews:2004}, and Dieterich-Ruina rate-state friction with an
+aging law \citep{Dieterich:1979}. See the PyLith manual
+\citep{PyLith:manual:1.7.1} for details.
% ------------------------------------------------------------------
\section{Finite-Element Mesh Processing}
-Like all finite-element engines, PyLith must be able to calculate cell
-and face integrals necessary to evaluate weak forms, assemble local
-cell vectors and matrices into global vector and matrix objects,
-impose Dirichlet boundary conditions on the algebraic system, and
-solve the resulting system of nonlinear algebraic equations. In
-PyLith, these operations are accomplished using PETSc, and in
-particular the Sieve package for finite-element support
-\citep{Knepley:Karpeev:2009}.
+Like all finite-element engines, PyLith performs operations on cells
+and vertices comprising the discretized domain (finite-element
+mesh). These operations include calculating cell and face integrals to
+evaluate weak forms, assemble local cell vectors and matrices into
+global vector and matrix objects, impose Dirichlet boundary conditions
+on the algebraic system, and solve the resulting system of nonlinear
+algebraic equations. In PyLith, these operations are accomplished
+using PETSc, and in particular the Sieve package for finite-element
+support \citep{Knepley:Karpeev:2009}.
The Sieve application programming interface (API) for mesh
representation and manipulation is based upon a direct acyclic graph
-representation of the \textit{covering} relation in a mesh. For
-example, faces cover cells, edges cover faces, and points cover
-edges. By focusing on the key topological relations, the interface can
-be both concise and quite general. Using this generic API, PyLith is
-able to support one, two, and three dimensional meshes, with
-simplicial, hex, and even prismatic cell shapes, while using very
-little dimension or shape specific code. However, in order to include
-faults, we include additional operations in Sieve beyond those
-necessary for conventional finite-element operations.
+representation of the \textit{covering} relation in a mesh. \brad{ADD
+ REFERENCE TO NEW FIGURE HERE} For example, faces cover cells, edges
+cover faces, and points cover edges. By focusing on the key
+topological relations, the interface can be both concise and quite
+general. Using this generic API, PyLith is able to support one, two,
+and three dimensional meshes, with simplicial, hex, and even prismatic
+cell shapes, while using very little dimension or shape specific
+code. However, in order to include faults, we include additional
+operations in Sieve beyond those necessary for conventional
+finite-element operations.
In our domain decomposition approach, the finite-element mesh includes
the fault as an interior surface. This forces alignment of the element
@@ -870,7 +817,7 @@ The introduction of Lagrange multipliers
The introduction of Lagrange multipliers to implement the fault slip
constraints produces the saddle-point problem shown in
equation~(\ref{eqn:saddle:point}). Traditional black-box parallel
-preconditioners, such as the Additive Schwarz Method (ASM)
+preconditioners, such as the additive Schwarz Method (ASM)
\citep{Smith:etal:1996}, are not very effective for this type of
problem and produce slow convergence. However, PETSc provides tools to
construct many variations of effective parallel preconditioners for
@@ -928,7 +875,7 @@ recompilation or special code.
We now turn our attention to evaluating the fault portion of the
preconditioning matrix associated with the Lagrange multipliers, since
-stock PETSc preconditioners can handle the elastic portion as
+PETSc preconditioners can handle the elastic portion as
discussed in the previous paragraph. In computing
$\bm{P_\mathit{fault}}$ we approximate $\bm{K}^{-1}$ with
the inverse of the diagonal portion of $\bm{K}$. $\bm{L}$ which
@@ -960,18 +907,18 @@ combine an AMG preconditioner for the el
combine an AMG preconditioner for the elasticity submatrix with our
custom fault preconditioner for the Lagrange multiplier submatrix. See
Section~\ref{sec:performance:benchmark} for a comparison of
-preconditioner performance for an application involved a static
+preconditioner performance for an application involving a static
simulation with multiple faults. It shows the clear superiority of
this setup over several other possible preconditioning strategies.
\subsection{Dynamic Simulations}
-In dynamic simulations the Courant-Friderichs-Lewy condition controls
-the stability of the explicit 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.
+In dynamic simulations the Courant-Friderichs-Lewy condition
+\citep{Courant:etal:1967} controls the stability of the explicit 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
@@ -1169,11 +1116,12 @@ We generate both hexahedral meshes and t
(available from http://cubit.sandia.gov) and construct meshes so that
the problem size (number of DOF) for the two different cell types
(hexahedra and tetrahedra) are nearly the same. The suite of
-simulations examine problem sizes increasing by about a factor of two
-from $1.8\times 10^5$ DOF (1 process) to $1.1\times 10^7$ DOF (64
-processes). The corresponding discretization sizes are 2033 m to 500 m
-for the hexahedral meshes and 2326 m to 581 m for the tetrahedral
-meshes. Figure~\ref{fig:solvertest:mesh} shows the 1846 m resolution
+simulations examine increasing larger problem sizes as we increase the
+number of processes, with $1.8\times 10^5$ DOF for 1 process up to
+$1.1\times 10^7$ DOF for 64 processes. The corresponding
+discretization sizes are 2033 m to 500 m for the hexahedral meshes and
+2326 m to 581 m for the tetrahedral meshes.
+Figure~\ref{fig:solvertest:mesh} shows the 1846 m resolution
tetrahedral mesh. As we will see in
Section~\ref{sec:verification:quasi-static}, the hexahedral mesh for a
given resolution in a quasi-static problem is slightly more accurate,
@@ -1195,7 +1143,7 @@ faulting and post-seismic and interseism
faulting and post-seismic and interseismic deformation.
For this benchmark of preconditioner performance, we examine the
-number of iterations required for convergence using the PETSc Additive
+number of iterations required for convergence using the PETSc additive
Schwarz (ASM), field split (with and without our custom
preconditioner), and Schur complement preconditioners discussed in
section~\ref{sec:solver:quasi-static}. We characterize the dependence
@@ -1275,11 +1223,12 @@ operations reported in the log summary \
operations reported in the log summary \cite{PETSc:manual}. These
operations are limited by available memory bandwidth rather than the
rate at which a processor can perform floating points operations. From
-Table~\ref{tab:solvertest:memory:events}, we see that we saturate the memory
-system using two processes per processor, since scaling plateaus from
-2 to 4 processes, but shows good scaling from 8 to 16 processes. This
-lack of memory bandwidth will depress overall performance, but should
-not affect the inter-node scaling of the application.
+Table~\ref{tab:solvertest:memory:events}, we see that we saturate the
+memory bandwidth using two processes per processor, since scaling
+plateaus from 2 to 4 processes, but shows good scaling from 8 to 16
+processes. This lack of memory bandwidth will depress overall
+performance, but should not affect the inter-node scaling of the
+application.
Machine network performance can be elucidated by the \texttt{VecMDot}
operation for vector reductions, and \texttt{MatMult} for
@@ -1302,14 +1251,14 @@ scale very well. However, that would in
scale very well. However, that would increase overall solver time and
thus would not be the choice to maximize scientific output.
-Figure~\ref{fig:solvertest:scaling} illustrates excellent the parallel
+Figure~\ref{fig:solvertest:scaling} illustrates the excellent parallel
performance for the finite-element assembly routines (reforming the
Jacobian sparse matrix and computing the residual). From the ASM
performance, we see that the basic solver building-blocks, like
parallel sparse matrix-vector multiplication, scale well. However, the
-ASM preconditioner itself is not scalable, and the number of iterates
-increases strongly with the number of processes. The introduction of Schur
-complement methods and an AMG preconditione slows the growth
+ASM preconditioner itself is not scalable, and the number of iterations
+increases significantly with the number of processes. The introduction
+of Schur complement methods and an AMG preconditioner slows the growth
considerably, but future work will pursue the ultimate goal of
iteration counts independent of the number of processes.
@@ -1317,40 +1266,40 @@ iteration counts independent of the numb
\section{Code Verification Benchmarks}
\label{sec:verification:benchmarks}
-In developing PyLith we verify the numerical implementation using a
-number of techniques. We employ unit testing to verify correct
-implementation of nearly all of the individual routines. Having a test
-for most object methods or functions isolates bugs at their origin
-during code development and prevents new bugs from occurring as code
-is modified or optimized. We also rely on full-scale benchmarks to
-verify that the code properly solves the numerical problem. These
-benchmark problems include quasi-static strike-slip and reverse
-viscoelastic simulations and various problems in the suite of
-dynamic spontaneous rupture benchmarks developed by the Southern
-California Earthquake Center (SCEC) and the United States Geological
-Survey \citep{Harris:etal:SRL:2009}. In this section we focus on two
-benchmarks that test different scientific applications: quasi-static
-relaxation of a Maxwell viscoelastic material subjected to multiple
-earthquake cycles involving slip and steady creep on a vertical
-strike-slip fault \citep{Savage:Prescott:1978} and supershear
+In developing PyLith we verify the numerical implementation of various
+features using a number of techniques. We employ unit testing to
+verify correct implementation of nearly all of the individual
+routines. Having a test for most object methods or functions isolates
+bugs at their origin during code development and prevents new bugs
+from occurring as code is modified or optimized. We also rely on
+full-scale benchmarks to verify that the code properly solves the
+numerical problem. These benchmarks include quasi-static strike-slip
+and reverse viscoelastic simulations and various exercises in the
+suite of dynamic spontaneous rupture benchmarks developed by the
+Southern California Earthquake Center (SCEC) and the United States
+Geological Survey \citep{Harris:etal:SRL:2009}. In this section we
+focus on two benchmarks that test different scientific applications:
+quasi-static relaxation of a Maxwell viscoelastic material subjected
+to multiple earthquake cycles involving slip and steady creep on a
+vertical strike-slip fault \citep{Savage:Prescott:1978} and supershear
dynamic spontaneous rupture of a 60 degree dipping normal fault in a
Drucker-Prager elastoplastic medium. This second benchmark corresponds
to benchmark TPV13 in the SCEC suite of dynamic spontaneous rupture
-benchmarks.
+benchmarks \citep{Harris:etal:SRL:2011}.
\subsection{Quasi-static Benchmark}
\label{sec:verification:quasi-static}
As a test of our quasi-static solution, we compare our numerical
results against the analytical solution of
-\citep{Savage:Prescott:1978}. This problem consists of an infinitely
+\citet{Savage:Prescott:1978}. This problem consists of an infinitely
long strike-slip fault in an elastic layer overlying a Maxwell
viscoelastic half-space. Figure~\ref{fig:savage:prescott::solution}
illustrates the geometry of the problem with an exaggerated view of
the deformation during the tenth earthquake cycle. Between earthquakes
the upper portion of the fault is locked, while the lower portion
slips at the plate velocity. At regular intervals (the earthquake
-recurrence time) the upper portion of the fault ruptures such that the
+recurrence time) the upper portion of the fault slides such that the
slip on the locked portion exactly complements the slip on the
creeping portion so the cumulative slip over an earthquake cycle is
uniform.
@@ -1368,17 +1317,19 @@ viscosity.
For this benchmark we use a locking depth of 20 km, an elastic layer
thickness of 40 km, an earthquake recurrence time of 200 years, a
-shear modulus of 30 GPa, a viscosity of $2.37\times 10^{19}$ Pa-s,
-and a relative plate velocity of 2 cm/year, implying a coseismic
-offset of 4 m every 200 years. The viscosity and shear modulus values
-yield a viscoelastic relaxation time of 50 years, and $\tau_0 =4$.
-We employ a 3-D model (2000 km by 1000 km by 400 km) with Dirichlet
-boundary conditions enforcing symmetry to approximate an infinitely
-long strike-slip fault. We apply velocity boundary conditions in the
-y-direction to the -x and +x faces with zero x-displacement. We
-constrain the vertical displacements on the bottom of the
-domain. Finally, we fix the x-displacements on the -y and +y faces to
-enforce symmetry consistent with an infinitely long strike-slip fault.
+shear modulus of 30 GPa, a viscosity of $2.37\times 10^{19}$ Pa-s, and
+a relative plate velocity of 2 cm/year, implying a coseismic offset of
+4 m every 200 years (see
+Table~\ref{tab:Savage:Prescott:parameters}). The viscosity and shear
+modulus values yield a viscoelastic relaxation time of 50 years, and
+$\tau_0 =4$. We employ a 3-D model (2000 km by 1000 km by 400 km)
+with Dirichlet boundary conditions enforcing symmetry to approximate
+an infinitely long strike-slip fault. We apply velocity boundary
+conditions in the y-direction to the -x and +x faces with zero
+x-displacement. We constrain the vertical displacements on the bottom
+of the domain to be zero. Finally, we fix the x-displacements on the
+-y and +y faces to enforce symmetry consistent with an infinitely long
+strike-slip fault.
We examine four different numerical solutions considering the effects
of cell type (hexahedral versus tetrahedral) and discretization
@@ -1404,23 +1355,24 @@ plate motion is approached after several
plate motion is approached after several earthquake cycles, once the
applied fault slip and velocity boundary conditions have produced
nearly steady flow in the viscoelastic half-space. It is therefore
-necessary to spin-up both solutions for several earthquake cycles to
-allow a comparison between the two. In this way, both models will have
-achieved steady-state behavior, and both models will have
-approximately the same component of steady plate motion. We simulate
-ten earthquake cycles for both the analytical and numerical models for
-a total duration of 2000 years. For the numerical solution we use a
-constant time step size of five years. This time step corresponds to
-one tenth of the viscoelastic relaxation time; hence it tests the
-accuracy of the viscoelastic solution for moderately large time steps
-relative to the relaxation time.
+necessary to spin-up both solutions to their steady-state solution
+over several earthquake cycles to allow a comparison between the
+two. In this way, the transient behavior present in both models will
+have nearly disappeared, and both models will have approximately the
+same component of steady plate motion. We simulate ten earthquake
+cycles for both the analytical and numerical models for a total
+duration of 2000 years. For the numerical solution we use a constant
+time step size of five years. This time step corresponds to one tenth
+of the viscoelastic relaxation time; hence it tests the accuracy of
+the viscoelastic solution for moderately large time steps relative to
+the relaxation time.
Figure~\ref{fig:savage:prescott:profiles} compares the numerical
results extracted on the ground surface along the center of the model
perpendicular to the fault with the analytic solution. Using a
logarithmic scale with distance from the fault facilitates examining
-the solution both close to the fault and in the far-field. For the
-second earthquake cycle, the far-field numerical solution does not yet
+the solution both close to and far from the fault. For the second
+earthquake cycle, the far-field numerical solution does not yet
accurately represent steady plate motion and the numerical simulations
underpredict the displacement. By the tenth earthquake cycle, steady
plate motion is accurately simulated and the numerical results match
@@ -1444,14 +1396,16 @@ allows the model to capture a more compl
\subsection{Dynamic Benchmark}
\label{sec:verification:dynamic}
-SCEC Spontaneous Rupture Benchmark TPV13 focuses on modeling an
-earthquake that produces extreme (very large) ground motions
-associated with supershear rupture towards the ground surface on a
-dipping fault with a large stress drop to generate large slip and fast
-slip rates \citep{Harris:etal:SRL:2011}. It uses a Drucker-Prager
-elastoplastic bulk rheology and a slip-weakening friction model in a
-depth-dependent initial stress field. Figure~\ref{fig:tpv13:geometry}
-show the geometry of the benchmark problem and the size of the domain
+As a test of PyLith's dynamic spontaneous rupture solutions, we use
+SCEC Spontaneous Rupture Benchmark TPV13 that models a high
+stress-drop, supershear, dip-slip earthquake that produces extreme
+(very large) ground motions, large slip, and fast slip rates
+\citep{Harris:etal:SRL:2011}. It uses a Drucker-Prager elastoplastic
+bulk rheology and a slip-weakening friction model in a depth-dependent
+initial stress field.
+
+Figure~\ref{fig:tpv13:geometry}
+show the geometry of the benchmark and the size of the domain
we used in our verification test. The benchmark includes both 2-D
(TPV13-2D is a vertical slice through the fault center-line with plane
strain conditions) and 3-D versions (TPV13). This benchmark specifies
@@ -1462,18 +1416,18 @@ for TPV13. We gradually coarsen the mesh
for TPV13. We gradually coarsen the mesh with distance from the fault
by increasing the discretization size at a geometric rate of 2\%. This
provides high resolution at the fault surface to resolve the small
-scale features of the rupture process with less resolution at the
+scale features of the rupture process and less resolution at the
edges of the boundary where the solution is much
smoother. Figure~\ref{fig:tpv13-2d:mesh} shows the triangular mesh for
a discretization size of 100 m on the fault.
Rupture initiates due to a low static coefficient of friction in the
-nucleation region. Figure~\ref{fig:tpv13-2d:stress:slip} illustrates
-the depth dependence of the stress field in terms of the fault
-tractions and Table~\ref{tab:tpv13:parameters} summarizes the
+nucleation region. Figure~\ref{fig:tpv13-2d:stress:slip}(a)
+illustrates the depth dependence of the stress field in terms of the
+fault tractions and Table~\ref{tab:tpv13:parameters} summarizes the
benchmark parameters. \citet{Harris:etal:SRL:2011} provides a more
complete description with all of the details available from
-\url{http://scecdata.usc.edu/cvws/cgi-bin/cvws.cgi}. An unfortunate
+\url{http://scecdata.usc.edu/cvws/cgi-bin/cvws.cgi}. A challenging
feature of this, and many other benchmarks in the SCEC Spontaneous
Rupture Code Verification Exercise, is the use of parameters with
spatial variations that are not continuous. This includes the
@@ -1483,16 +1437,16 @@ construction of the finite-element mesh
construction of the finite-element mesh and use the spatial average of
the parameters where they are discontinuous. This decreases the
sensitivity of the numerical solution to the discretization size. This
-benchmark also includes fluid pressures. Because PyLith does not
-include fluid pressure, we formulate the simulation parameters in
-terms of effective stresses.
+SCEC benchmark also includes fluid pressures. Because PyLith does not
+include fluid pressure, we instead formulate the simulation parameters
+in terms of effective stresses.
-Figure~\ref{fig:tpv13-2d:stress:slip} displays the final slip
+Figure~\ref{fig:tpv13-2d:stress:slip}(b) displays the final slip
distribution in the TPV13-2D simulation with triangular cells at a
resolution of 100 m. The large dynamic stress drop and supershear
rupture generate 20 m of slip at a depth of about 7
-km. Figure~\ref{fig:tpv13-2d:slip:rate}(a) demonstrates the convergence
-of the solution as the discretization size decreases as evidence in
+km. Figure~\ref{fig:tpv13-2d:slip:rate}(a)--(d) demonstrates the convergence
+of the solution as the discretization size decreases as evident in
the normal faulting component of fault slip rate time histories. For a
resolution of 200 m on the fault, the solution contains some
high-frequency oscillation due to insufficient resolution of the
@@ -1501,20 +1455,21 @@ oscillation in the slip rate time histor
oscillation in the slip rate time histories. The triangular cells
generate less oscillation compared with quadrilateral cells.
-In this problem without an analytical solution, as in all of the
-benchmarks in the SCEC spontaneous rupture benchmark suite, we rely on
+In this benchmark without an analytical solution, as in all of the
+exercises in the SCEC spontaneous rupture benchmark suite, we rely on
comparison with other dynamic spontaneous rupture modeling codes to
verify the numerical implementation in
-PyLith. Figure~\ref{fig:tpv13-2d:slip:rate}(b) compares the slip rate
-time histories from PyLith with four other codes (see
-\citep{Harris:etal:SRL:2011}, \citep{Andrews:etal:2007},
-\citep{Barall:2009}, and \citep{Dunham:etal:2011} for a discussion of
-these other finite-element and finite-difference codes). The slip
-rate time histories agree very well, although some codes yield more
-oscillation than others. We attribute this to variations in the amount
-of numerical damping used across the various codes.
+PyLith. Figure~\ref{fig:tpv13-2d:slip:rate}(e)--(h) compares the slip
+rate time histories from PyLith with four other codes (see
+\citet{Harris:etal:SRL:2011}, \citet{Andrews:etal:2007},
+\citet{Barall:2009}, \citet{Ma:2009}, and \citet{Dunham:etal:2011} for
+a discussion of these other finite-element and finite-difference
+codes). The slip rate time histories agree very well, although some
+codes yield more oscillation than others. We attribute this to
+variations in the amount of numerical damping used in the various
+codes.
-The results in for the 3-D version of the TPV13 benchmark yield
+The results for the 3-D version of the TPV13 benchmark yield
similar results. Figure~\ref{fig:tpv13:rupture:time}(a) shows the same
trends in rupture speed with discretization size that we observed in
the 2-D version. In both cases models with insufficient resolution to
@@ -1540,16 +1495,16 @@ consistent with the trends observed in t
consistent with the trends observed in the comparison of rupture
times. Furthermore, the codes all produce consistent results
throughout the entire time histories. The small differences in rupture
-time in the mode-II (along-strike) direction between the two groups of
+time in the mode-III (along-strike) direction between the two groups of
codes is evident in the slip rate time histories at a depth of 7.5 km
and 12 km along strike
(Figure~\ref{fig:tpv13:slip:rate}(f)). Nevertheless, this simply
produces a small time shift in the time history.
From the 2-D and 3-D versions of the SCEC spontaneous rupture
-benchmark TPV13, we conclude that the PyLith performs very similarly
-to other finite-element and finite-difference dynamic
-spontaneous rupture modeling codes. In particular it is well-suited to problems
+benchmark TPV13, we conclude that PyLith performs similarly
+to other finite-element and finite-difference dynamic spontaneous
+rupture modeling codes. In particular it is well-suited to problems
with complex geometry as we are able to vary the discretization size
while simulating a dipping normal fault. The code accurately captures
supershear rupture and properly implements a Drucker-Prager
@@ -1559,25 +1514,24 @@ elastoplastic bulk rheology and slip-wea
\section{Conclusions}
\label{sec:conclusions}
-The domain decomposition approach to implementing fault slip in PyLith
-provides a flexible numerical implementation as evident by its
-application in quasi-static and dynamic simulations and kinematic
-(prescribed) slip and spontaneous rupture (frictional sliding) models
-of earthquake rupture. We find that algebraic multigrid
-preconditioners for elasticity combined with a custom preconditioner
-for the fault block associated with the Lagrange multipliers
-accelerates the convergence of the Krylov solver with the fewest
-number of iterations and the least sensitivity to problem
-size. Benchmark tests demonstrate the accuracy of our implementation
-in PyLith with excellent agreement with the analytical solution for
-viscoelastic relaxation and strike-slip faulting over multiply
-earthquake cycle and excellent agreement with other codes for
-supershear dynamic spontaneous rupture on a dipping normal fault
-embedded in an elastoplastic domain. Consequently, we believe this
-methodology provides a promising avenue for modeling the
-earthquake cycle through coupling of quasi-static simulations of the
-interseismic and postseismic deformation and dynamic rupture
-simulations of earthquake rupture propagation.
+PyLith provides a flexible numerical implementation of fault slip
+using a domain decomposition approach. We have evaluated the
+efficiency of several preconditioners for use of this fault
+implementation in quasi-static simulations. We find that algebraic
+multigrid preconditioners for elasticity combined with a custom
+preconditioner for the fault block associated with the Lagrange
+multipliers accelerates the convergence of the Krylov solver with the
+fewest number of iterations and the least sensitivity to problem
+size. Benchmark tests demonstrate the accuracy of our fault slip
+implementation in PyLith with excellent agreement to (1) an analytical
+solution for viscoelastic relaxation and strike-slip faulting over
+multiple earthquake cycles and (2) other codes for supershear dynamic
+spontaneous rupture on a dipping normal fault embedded in an
+elastoplastic domain. Consequently, we believe this methodology
+provides a promising avenue for modeling the earthquake cycle through
+coupling of quasi-static simulations of the interseismic and
+postseismic deformation and dynamic rupture simulations of earthquake
+rupture propagation.
% ----------------------------------------------------------------------
% Notation -- End each entry with a period.
@@ -1673,10 +1627,10 @@ simulations of earthquake rupture propag
the algebraic multigrid preconditioner and fault block custom
preconditioner. The finite-element integrations for the Jacobian
and residual exhibit good weak scaling with minimal sensitivity to
- the problem size. The linear solve does not scale as well, which
- we attribute to the poor scaling of the algebraic multigrid setup
- and application as well as limited memory and interconnect
- bandwidth.}
+ the problem size. The linear solve (solid lines in the top panel)
+ does not scale as well, which we attribute to the poor scaling of
+ the algebraic multigrid setup and application as well as limited
+ memory and interconnect bandwidth.}
\label{fig:solvertest:scaling}
\end{figure}
@@ -1726,7 +1680,7 @@ simulations of earthquake rupture propag
a Drucker-Prager elastoplastic bulk rheology, slip-weakening
friction, a depth-dependent stress field, and normal fault with a
60 degree dip angle. The 2-D version corresponds to the vertical
- slice shown by the dashed line. The red dotes denote locations on
+ slice shown by the dashed line. The red dots denote locations on
the fault used in the comparison of the vertical slip rates
(Figures~\ref{fig:tpv13-2d:slip:rate}
and~\ref{fig:tpv13:slip:rate}). The blue dots indicate locations
@@ -1911,7 +1865,7 @@ simulations of earthquake rupture propag
& Hex8 & 35 & 39 & 43 \\
\hline
\end{tabular}
-\tablenotetext{a}{Number of iterations for Additive Schwarz (ASM),
+\tablenotetext{a}{Number of iterations for additive Schwarz (ASM),
Schur complement (Schur), and field split (additive, multiplicative,
and multiplicative with custom fault block preconditioner),
preconditioners for tetrahedral and hexahedral discretizations and
@@ -2052,6 +2006,36 @@ VecMAXPY & 1 & 1.0 & 1733 \\
\clearpage
\begin{table}
+\caption{Savage and Prescott Benchmark Parameters\tablenotemark{a}}
+\label{tab:Savage:Prescott:parameters}
+\centering
+\begin{tabular}{llc}
+ \hline
+ \multicolumn{2}{l}{Parameter} & Value \\
+ \hline
+ \multicolumn{2}{l}{Domain} & \\
+ & Length & 2000 km \\
+ & Width & 1000 km \\
+ & Height & 400 km \\
+ & Fault dip angle & 90 $\deg$ \\
+ \multicolumn{2}{l}{Physical properties} & \\
+ & Shear modulus & 30 GPa \\
+ & Viscosity & $2.37\times10^{19}$ Pa-s \\
+ & Elastic thickness & 40 km \\
+ \multicolumn{2}{l}{Fault slip} & \\
+ & Locking depth & 20 km \\
+ & Earthquake recurrence & 200 yr \\
+ & Earthquake slip & 4 m \\
+ & Creep rate & 2 cm/yr \\
+ \hline
+\end{tabular}
+\tablenotetext{a}{Simulation parameters for the Savage and Prescott benchmark
+ with multiple earthquake cycles on a vertical strike-slip fault
+ embedded in an elastic layer over a viscoelastic half-space.}
+\end{table}
+
+
+\begin{table}
\caption{SCEC Benchmark TPV13 Parameters\tablenotemark{a}}
\label{tab:tpv13:parameters}
\centering
@@ -2070,10 +2054,9 @@ VecMAXPY & 1 & 1.0 & 1733 \\
& Density ($\rho$) & 2700. kg/m$^3$ \\
\hline
\end{tabular}
-\tablenotetext{a}{Simulation parameters for the performance benchmark
- with three faults embedded in a volume domain as shown in
- Figure~\ref{fig:solvertest:geometry}. We prescribe right-lateral
- (RL) slip on the middle fault and left-lateral (LL) slip on the end faults.}
+\tablenotetext{a}{Basic simulation parameters for the SCEC
+ dynamic spontaneous rupture benchmark TPV13. Full specification of
+ the parameters can be found in \citet{Harris:etal:SRL:2011}.}
\end{table}
diff -r 91049a93b5e2 -r 6b10724d5da2 references.bib
--- a/references.bib Tue Aug 14 13:29:49 2012 -0700
+++ b/references.bib Tue Aug 28 12:02:36 2012 -0700
@@ -293,9 +293,92 @@
data},
}
+ at article{Harris:Day:1999,
+ author = {Harris, R.~A. and Day, S.~M.},
+ title = {Dynamic 3-{D} simulations of earthquakes on en
+ echelon faults},
+ journal = GRL,
+ volume = 26,
+ number = 14,
+ year = 1999,
+ month = jul # {~15},
+ pages = {2089--2092},
+ abstract = {One of the mysteries of earthquake mechanics is why
+ earthquakes stop. This process determines the
+ difference between small and devastating
+ ruptures. One possibility is that fault geometry
+ controls earthquake size. We test this hypothesis
+ using a numerical algorithm that simulates
+ spontaneous rupture propagation in a
+ three-dimensional medium and apply our knowledge to
+ two California fault zones. We find that the size
+ difference between the 1934 and 1966 Parkfield,
+ California, earthquakes may be the product of a
+ stepover at the southern end of the 1934 earthquake
+ and show how the 1992 Landers, California,
+ earthquake followed physically reasonable
+ expectations when it jumped across en echelon faults
+ to become a large event. If there are no linking
+ structures, such as transfer faults, then
+ strike-slip earthquakes are unlikely to propagate
+ through stepovers $>$5 km wide.},
+}
+
+ at article{Mikumo:etal:1998,
+ author = {Mikumo, T. and Miyatake, T. and Santoyo, M.~A.},
+ title = {Dynamic rupture of asperities and stress change
+ during a sequence of large interplate earthquakes in
+ the Mexican subduction zone},
+ journal = BSSA,
+ volume = 88,
+ number = 3,
+ month = jun,
+ year = 1998,
+ pages = {686--702},
+ abstract = {We investigate the spatial and temporal variations
+ of shear stress due to the successive failures over
+ an extensive segment of the Mexican subduction zone
+ during a sequence of large interplate earthquakes
+ that occurred over a period of 13 yr. For this
+ purpose, we develop 3D dynamic rupture models
+ incorporating a shallowly dipping fault located
+ above the subducting plate. The spatial distribution
+ of dynamic stress drop over the fault has been
+ estimated for each of the events, through an
+ inversion procedure using some of the previously
+ derived kinematic fault parameters as observational
+ constraints. The results revealed quite
+ heterogeneous stress changes during these
+ earthquakes coming from medium to high dynamic
+ stress drop due to the rupture of a few patchlike
+ asperities and from stress increase in between and
+ around them. Two weak asperities located southeast
+ of the Michoacan segment were ruptured first by the
+ 1979 Petatlan event. The 1981 Playa Azul event
+ ruptured two asperities in the central zone with a
+ stress drop higher than 80 bars. The largest 1985
+ Michoacan earthquake resulted from the rupture of
+ two large-size, strong asperities located at both
+ sides of the 1981 fault zone with high stress drop
+ of XO to 100 bars and from another two asperities at
+ depth. Two days after this largest event, two
+ asperities were broken during the Zihuatanejo
+ aftershock in the southeastern adjacent zone. Many
+ aftershocks of these large events tend to be
+ distributed in the zones of stress increase outside
+ the asperities, while only small numbers of
+ aftershocks have been observed within these asperity
+ zones, It appears that several major asperities that
+ existed in this extensive segment have been ruptured
+ successively so as to fill unbroken gaps on the
+ plate interface. Thus, the stress change left over
+ from the previous earthquake has dominant effects on
+ the next event in this subduction zone.},
+}
+
@Article{Andrews:etal:2007,
author = {Andrews, D.~J. and Hanks, T.~C. and Whitney, J.~W.},
- title = {Physical limits on ground motion at Yucca Mountain},
+ title = {Physical limits on ground motion at {Yucca} {Mountain}},
journal = BSSA,
year = {2007},
volume = {97},
@@ -340,6 +423,39 @@
month = oct,
pages = {2308--2322},
doi = {10.1785/0120100075}
+}
+
+ at article{Ma:2009,
+ author = {Ma, S.},
+ title = {Distinct asymmetry in rupture-induced inelastic strain
+ across dipping faults: {An} off-fault yielding
+ model},
+ journal = GRL,
+ volume = 36,
+ number = {L20317},
+ year = 2009,
+ doi = {10.1029/2009GL040666},
+ abstract = {Off-fault stresses associated with rupture
+ propagation cause material yielding under a
+ Mohr-Coulomb condition, leaving irrecoverable
+ deformation in the fault zone. 2D dynamic rupture
+ simulations for a 30° reverse fault and a 60° normal
+ fault in a depth-dependent stress environment show
+ that the yielding (inelastic) zone off the fault
+ broadens as it nears the surface with decreasing
+ confining pressure, forming a skewed ‘flower-like’
+ structure bounded at the top by the free
+ surface. The inelastic strain in the hanging wall is
+ significantly larger and broader than the footwall
+ for both reverse and normal faults. The occurrence
+ of inelastic strain, however, significantly reduces
+ ground motion (especially on the hanging wall) and
+ gives rise to a reduced asymmetry in ground motion
+ on the hanging wall and footwall compared to elastic
+ solutions. These results provide theoretical
+ predictions for fault zone structure of dipping
+ faults that can be tested by future observational
+ experiments.},
}
@Article{Day:etal:2005,
@@ -930,6 +1046,121 @@
be far more aperiodic than has been suggested. },
}
+ at Article{Robinson:Benites:1995,
+ author = {Robinson, R. and Benites, R.},
+ title = {Synthetic seismicity models of multiple interacting faults},
+ journal = JGR-SE,
+ year = {1995},
+ volume = {100},
+ number = {B9},
+ pages = {18,229--18,238},
+ doi = {10.1029/95JB01569},
+ abstract = {A synthetic seismicity model for multiple,
+ interacting faults, of any strike and dip in a
+ three-dimensional elastic half-space has been
+ developed. The ultimate purpose is to determine how
+ hazard is modified by elastic interactions of events
+ in complex fault systems and to examine the
+ synthetic seismicity for unusual behavior before
+ large events that might be observable in the real
+ world. Each fault is subdivided into some number of
+ equal-sized patches, and a coefficient of static
+ friction, random within a specified range, is
+ assigned to each. The patch static strength is then
+ the product of friction and differential normal
+ stress. Stress (both shear and normal) accumulates
+ due to some predefined tectonic driving force until
+ one or more patches fail and slip
+ instantaneously. The slip can be in any direction
+ within the fault plane and its magnitude is
+ determined by a specified fractional stress
+ drop. This initial slip induces shear stress and
+ strength changes on all other patches and faults and
+ perhaps results in other patches failing. Patches
+ can slip more than once; once having failed, their
+ strength is reduced, resulting in an overshoot that
+ approximates the effect of dynamic stress
+ enhancement. An event is over when all patches are
+ stable. Although present computational constraints
+ do not allow models with more than about 1500
+ patches, or a large number of faults, experiments
+ with simple models of two parallel strike-slip
+ faults (25 km length, 15 km depth, 694 m × 714 m
+ patch size), driven at different rates and separated
+ by 10 or 15 km, give some indicative results. A
+ range of friction between 0.1 and 0.2, with a stress
+ drop of 10\% (about 3 MPa or 30 bars), produces a b
+ value (500,000 events) close to 1. The models
+ generate a distinct class of large “characteristicâ€
+ events that rupture the full fault plane. The
+ distribution of interevent times for the background
+ seismicity is very similar to that for a Poisson
+ process, while the characteristic events are
+ quasiperiodic. Nearly half the larger events
+ directly induce some activity on the other
+ fault. The effect of a characteristic event on the
+ occurrence of another characteristic event on the
+ opposite fault is small and depends on the
+ separation: at 10 km distance other events are
+ retarded while at 15 km distance they are
+ advanced. At 10-km separation the probability of a
+ pair of characteristic events within 1 year is
+ reduced to 0.3\% from 1.7\% for the case with no
+ interaction. Consideration of interactions between
+ slip on the faults and the driving mechanism, as
+ well as higher stress drops and more overshoot,
+ would increase the magnitude of the interaction
+ effects.},
+}
+
+ at Article{Robinson:Benites:1996,
+ author = {Robinson, R. and Benites, R.},
+ title = {Synthetic seismicity models for the {Wellington}
+ Region, {New} {Zealand}: {Implications} for the
+ temporal distribution of large events},
+ journal = JGR-SE,
+ year = {1996},
+ volume = {101},
+ number = {B12},
+ pages = {27,833--27,844},
+ doi = {10.1029/96JB02533},
+ abstract = {Our previous synthetic seismicity model for multiple
+ interacting faults in a three-dimensional half-space
+ has been extended and applied to the Wellington
+ region, New Zealand, generating long catalogs of
+ earthquakes for studying the effects of elastic
+ interactions on the temporal distribution of large
+ events. The region is one of oblique plate
+ convergence with the interface between the
+ subducting and overlying plates at an average depth
+ of 22 km. Faults included, besides the subduction
+ thrust, are segments of the four major arc-parallel,
+ strike-slip faults overlying the plate interface. We
+ have used 40 different models, with geometric and
+ mechanical parameters chosen at random within
+ reasonable ranges, to generate catalogs of 200,000
+ years duration each. For comparison, each model was
+ rerun with the elastic interactions
+ suppressed. Considering events of magnitude 7.2 or
+ more (“characteristic†events in the sense that they
+ rupture most of a fault plane), the number of short
+ (<10 years) interevent times is higher than for the
+ corresponding case with no interactions, for all
+ models but one: the ratio ranges from 0.94 to 37.21
+ (9.46 average). For longer interevent times (10 to
+ 250 years) the relative numbers are in the opposite
+ sense. For still longer interevent times (>250
+ years), the relative numbers are again mostly
+ higher. Experiments with simple models indicate that
+ this pattern requires both interfault enhancement
+ and inhibition of large events. Mutual enhancement
+ occurs most often between the subduction thrust and
+ the overlying strike-slip faults and between the two
+ segments of the Wellington fault that almost join
+ end to end. Mutual inhibition mostly occurs between
+ the subparallel strike-slip faults.},
+}
+
@article{Savage:Prescott:1978,
author = {Savage, J.~C. and Prescott, W.~H.},
title = {Asthenosphere Readjustment and the Earthquake Cycle},
@@ -1136,8 +1367,80 @@
margins (9\%-30\%).},
}
+ at article{Ida:1972,
+ author = {Ida, Yoshiaki},
+ title = {Cohesive Force across the Tip of a
+ Longitudinal-Shear Crack and {G}riffith's Specific
+ Surface Energy},
+ journal = JGR,
+ year = 1972,
+ volume = 77,
+ number = 20,
+ pages = {3796--3805},
+}
+
+ at Article{Andrews:2004,
+ author = {Andrews, D.~J.},
+ title = {Rupture calculations with dynamically-determined
+ slip-weakening friction},
+ journal = BSSA,
+ year = 2004,
+ volume = 94,
+ number = 3,
+ pages = {769--775},
+ month = jun,
+ doi = {10.1785/0120030142},
+ abstract = {The critical breakdown displacement, Dc, in which
+ friction drops to its sliding value, can be made
+ dependent on event size by specifying friction to be
+ a function of variables other than slip. Two such
+ friction laws are examined here. The first is
+ designed to achieve accuracy and smoothness in
+ discrete numerical calculations. Consistent
+ resolution throughout an evolving rupture is
+ achieved by specifying friction as a function of
+ elapsed time after peak stress is reached. Such a
+ time-weakening model produces Dc and fracture energy
+ proportional to the square root of distance rupture
+ has propagated in the case of uniform stress
+ drop. The second friction law is more physically
+ motivated. Energy loss in a damage zone outside the
+ slip zone has the effect of increasing Dc and
+ limiting peak slip velocity (Andrews, 1976). This
+ article demonstrates a converse effect, that
+ artificially limiting slip velocity on a fault in an
+ elastic medium has a toughening effect, increasing
+ fracture energy and Dc proportionally to rupture
+ propagation distance in the case of uniform stress
+ drop. Both the time-weakening and the
+ velocity-toughening models can be used in
+ calculations with heterogeneous stress drop.},
+}
+
+ at article{Dieterich:1979,
+ author = {Dieterich, J.~H.},
+ title = {Modeling of rock friction, 1. {Experimental} results
+ and constitutive equations},
+ journal = JGR-SE,
+ volume = 84,
+ number = {B5},
+ year = 1979,
+ month = may # {~10},
+ pages = {2161--2168},
+}
+
+ at article{Newmark:1959,
+ author = {Newmark, N.~M.},
+ title = {A method of computation for structural dynamics},
+ journal = {Journal of Engineering Mechanics},
+ volume = {85},
+ number = {EM3},
+ pages = {67--94},
+ year = {1959},
+}
+
@Book{Zienkiewicz:Taylor:2005,
- author = {Ziekiewicz, O.~C. and Taylor, R.~L. and Zhu, J.~Z.},
+ author = {Zienkiewicz, O.~C. and Taylor, R.~L. and Zhu, J.~Z.},
title = {The Finite-Element Method: Its Basis and Fundamentals},
publisher = {Butterworth-Heinemann},
year = {2005},
@@ -1204,6 +1507,18 @@
year = {2004}
}
+ at article{Courant:etal:1967,
+ author = {Courant, R. and Friedrichs, K. and Lewy, H.},
+ title = {On the Partial Difference Equations of Mathematical Physics},
+ journal = {IBM Journal of Research and Development},
+ volume = {11},
+ number = {2},
+ year = 1967,
+ month = mar,
+ pages = {215--234},
+ note = {English translation of the original 1928 paper published in {\it Mathematische Annalen}}
+}
+
@inproceedings{GordonBell09,
author = {Kaushik, D. and Smith, M. and Wollaber, A. and Smith, B. and Siegel, A. and Yang, W.~S.},
title = {Enabling High Fidelity Neutron Transport Simulations on Petascale Architectures},
More information about the CIG-COMMITS
mailing list