[cig-commits] commit: Edits to numerical model and dynamic stuff. Edits to Matt's sections. Started adding placeholders for figures.
Mercurial
hg at geodynamics.org
Tue Sep 6 13:43:07 PDT 2011
changeset: 73:20748e3ccedf
tag: tip
user: Brad Aagaard <baagaard at usgs.gov>
date: Tue Sep 06 13:43:04 2011 -0700
files: faultRup.tex references.bib
description:
Edits to numerical model and dynamic stuff. Edits to Matt's sections. Started adding placeholders for figures.
diff -r 26022948ab53 -r 20748e3ccedf faultRup.tex
--- a/faultRup.tex Thu Sep 01 20:33:15 2011 -0500
+++ b/faultRup.tex Tue Sep 06 13:43:04 2011 -0700
@@ -187,9 +187,9 @@ In this section we summarize the formula
In this section we summarize the formulation of the governing
equations using the finite-element method. We augment the conventional
finite-element formulation for elasticity with a domain decomposition
-approach \cite{Zienkiewicz:Taylor:2005,SmithBjorstadGropp96} to implement the fault slip.
-\citeN{PyLith:manual:1.6.2} provides a step-by-step description of the
-formulation.
+approach \cite{Smith:etal:1996,Zienkiewicz:Taylor:2005} to implement
+the fault slip. The PyLith manual \cite{PyLith:manual:1.6.2} provides
+a step-by-step description of the formulation.
We solve the elasticity equation including inertial terms,
\begin{gather}
@@ -201,15 +201,15 @@ We solve the elasticity equation includi
\text{ on }S_f, \label{eq:faultDisp}
\end{gather}
where $\vec{u}$ is the displacement vector, $\rho$ is the mass
-density, $f$ is the body force vector, $\tensor{\sigma}$ is the Cauchy
-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
+density, $\vec{f}$ is the body force vector, $\tensor{\sigma}$ is the
+Cauchy 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 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
-types of boundary conditions simultaneously.
+and $S_u$; however, a degree of freedom at any location cannot be
+associated with both types of boundary conditions simultaneously.
Following a conventional finite-element formulation (ignoring the
fault surface for a moment), we construct the
@@ -238,11 +238,23 @@ direction to this interior boundary and
direction to this interior boundary and ``positive'' and ``negative''
labels to the two sides of the fault, such that the fault normal is
the vector from the negative side of the fault to the positive side of
-the fault. Slip on the fault is the motion on the positive side
-relative to the motion on the negative side. Slip on the fault also
-corresponds to equal and opposite tractions on the positive and
-negative sides of the fault, which we impose using Lagrange
-multipliers with $\vec{l}_{+} - \vec{l}_{-} = 0$.
+the fault. Slip on the fault is the displacement of the positive side
+relative to the negative side. Slip on the fault also corresponds to
+equal and opposite tractions on the positive and negative sides of the
+fault, which we impose using Lagrange multipliers with $\vec{l}_{+} -
+\vec{l}_{-} = 0$.
+
+
+\begin{figure}[htbp]
+ \centering
+ \brad{TODO: ADD DIAGRAM}
+ %\includegraphics{figs/domain_decomposition}
+ \caption{Diagram of domain decomposition approach for modeling fault
+ slip. The fault slip introduces a jump in the displacement field
+ across the fault, whereas the tractions are continuous.}
+ \label{fig:domain:decomposition}
+\end{figure}
+
Recognizing that the tractions on the fault surface are analogous to
the boundary tractions, we add in the contributions from integrating
@@ -264,14 +276,13 @@ surface to zero,
\int_{S_f} \vec{\phi} \cdot
\left( \tensor{R} \cdot \vec{d} - \vec{u}_{+} + \vec{u}_{-} \right) \, dS = 0.
\end{equation}
-The rotation matrix $\tensor{R}$ is comprised of blocks of direction cosines.
We express the weighting function $\vec{\phi}$, trial solution
$\vec{u}$, and Lagrange multipliers $\vec{l}$ as linear combinations
of basis functions,
\begin{gather}
-\vec{\phi} = \sum_{n} \vec{a}_m N_m, \\
-\vec{u} = \sum_{m} \vec{u}_n N_n, \\
+\vec{\phi} = \sum_{m} \vec{a}_m N_m, \\
+\vec{u} = \sum_{n} \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
@@ -318,9 +329,9 @@ equations in terms of the increment in t
equations in terms of the increment in the solution from time $t$ to
$\tpdt$ rather than the solution at time $\tpdt$.
Consequently, rather than constructing a system with the form
-$\tensor{A} \vec{u}(\tpdt) = \vec{b}(\tpdt)$, we construct a system
-with the form $\tensor{A} d\vec{u} = \vec{b}(\tpdt) - \tensor{A}
-\vec{u}(t)$, where $\vec{u}(\tpdt) = \vec{u}(t) + d\vec{u}$.
+$\tensor{A} \cdot \vec{u}(\tpdt) = \vec{b}(\tpdt)$, we construct a system
+with the form $\tensor{A} \cdot d\vec{u} = \vec{b}(\tpdt) - \tensor{A}
+\cdot \vec{u}(t)$, where $\vec{u}(\tpdt) = \vec{u}(t) + d\vec{u}$.
% ------------------------------------------------------------------
\subsection{Quasi-static Simulations}
@@ -356,10 +367,10 @@ tools for solving linear systems of alge
tools for solving linear systems of algebraic equations with parallel
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
+- \tensor{A} \cdot \vec{u}$ and the Jacobian of the system
+($\tensor{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}
@@ -381,14 +392,14 @@ and~(\ref{eqn:quasi-static:residual:faul
\end{array} \right).
\end{equation}
-The Jacobian of the system is the action (operations) on the increment
+The Jacobian of the system is the action (operation) on the increment
in the solution. To find the portion of the Jacobian associated with
equation~(\ref{eqn:quasi-static:residual:domain}), we let
$\tensor{\sigma}(\tpdt) = \tensor{\sigma}(t) +
d\tensor{\sigma}(t)$. The action on the increment of the solution is
-associated with the term $d\tensor{\sigma}(t)$. We approximate the
-increment in the stress tensor using linear elasticity and
-infinitesimal strains,
+associated with the increment in the stress tensor
+$d\tensor{\sigma}(t)$. We approximate the increment in the stress
+tensor using linear elasticity and infinitesimal strains,
\begin{equation}
d\tensor{\sigma}(t) = \frac{1}{2} \tensorfour{C}(t) \cdot (\nabla + \nabla^T)
\vec{u}(t),
@@ -400,7 +411,7 @@ some algebra we find this portion of the
some algebra we find this portion of the Jacobian is
\begin{equation}
\label{eqn:jacobian:implicit:stiffness}
- \tensor{K} = \frac{1}{4} \int_v
+ \tensor{K} = \frac{1}{4} \int_V
(\nabla^T + \nabla) \tensor{N}_m^T \cdot
\tensorfour{C} \cdot (\nabla + \nabla^T) \tensor{N}_n \, dV.
\end{equation}
@@ -423,11 +434,11 @@ Thus, the Jacobian of the entire system
In dynamic simulations we include the inertial term in order to
resolve the propagation of seismic waves. The integral equation for
-fault slip constraint remains unchanged, so the corresponding portions
-of the Jacobian and residual are exactly the same as in the
-quasi-static simulations. Thus, the residual is the same as in the
-quasi-static simulations except the upper portion now includes the
-inertial term:
+the fault slip constraint remains unchanged, so the corresponding
+portions of the Jacobian ($\tensor{L}$) and residual ($\vec{r}_p$) are
+exactly the same as in the quasi-static simulations. Hence, the
+residual is the same as in the quasi-static simulations except the
+upper portion now includes the inertial term:
\begin{equation}
\vec{r} = \left( \begin{array}{c}
\begin{aligned}
@@ -459,7 +470,7 @@ difference scheme, so that the accelerat
\frac{\partial^2 \vec{u}}{\partial t^2}(t) =
\frac{1}{\Delta t^2} \left(
d\vec{u} - \vec{u}(t) + \vec{u}(\tmdt)
- \right) \\
+ \right), \\
%
\frac{\partial \vec{u}}{\partial t}(t) = \frac{1}{2\Delta t} \left(
d\vec{u} + \vec{u}(t) - \vec{u}(\tmdt)
@@ -484,19 +495,19 @@ so that the upper portion of the Jacobia
% ------------------------------------------------------------------
% Taking advantage of the commonality
-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 integral terms for the bulk
+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, so all of the various
-elastic, viscoelastic, and viscoelastoplastic 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
+elastic, viscoelastic, and viscoelastoplastic 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 compared to maintaining different codes for
-quasi-static and dynamic simulations.
+development time and simplifies maintenance compared to maintaining
+different codes for quasi-static and dynamic simulations.
% ------------------------------------------------------------------
@@ -505,7 +516,7 @@ In a prescribed (kinematic) fault ruptur
In a prescribed (kinematic) fault rupture we specify the slip-time
history $\vec{d}(x,y,z,t)$ at every location on the fault surfaces.
The slip-time history enters into the calculation of the residual as
-do the Lagrange multipliers, which are available as from the current
+do the Lagrange multipliers, which are available from the current
solution. In prescribing the slip-time history we do not specify the
initial tractions on the fault surface so the Lagrange multipliers are
the change in the tractions on the fault surfaces corresponding to the
@@ -527,16 +538,14 @@ the boundary conditions and deformation.
the boundary conditions and deformation. Our fault implementation
involves specifying the fault slip and solving for the tractions on
the fault (Lagrange multipliers). Hence, the fault constitutive model
-places bounds on the Lagrange multipliers and the systems of equations
-is nonlinear- we must adjust the slip to be consistent with the fault
-tractions. When a location on a fault is not slipping, no updates to
-the specified fault slip are necessary and the solution will converge
-in a single iteration of the nonlinear solve (i.e., the nonlinear
-problem simplifies to a linear one). In our quasi-static simulations,
-in general, the fault slips only a small fraction of the simulation
-duration, so the fault implementation yields a linear problem during
-most of the simulation. In dynamic simulations, we employ an
-alternative approach that involves a lumped Jacobian and avoids a
+places bounds on the Lagrange multipliers and the system of equations
+is nonlinear-- when a location on the fault is slipping, we must
+update the slip to be consistent with the fault tractions. When a
+location on a fault is not slipping, no updates to the specified fault
+slip are necessary and the solution will converge in a single
+iteration of the nonlinear solve (i.e., the nonlinear problem
+simplifies to a linear one). In dynamic simulations, we employ an
+alternative approach that involves a diagonal Jacobian and avoids a
nonlinear system of equations.
We determine the adjustment to the fault slip for a given perturbation
@@ -586,8 +595,9 @@ the fault (which occurs in most simulati
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).
+(which occurs in the rare cases in which the domain has a
+through-going fault). The nonlinear solve requires fewer iterations
+as the estimate of the fault slip becomes more accurate.
PyLith includes several commonly used fault constitutive models, all
of which specify the shear traction on the fault $T_f$ as a function
@@ -601,53 +611,110 @@ direction of the slip rate. We use the s
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).
+time-weakening, and Dieterich-Ruina rate-state friction with an aging
+law (see the PyLith manual \cite{PyLith:manual:1.6.2} 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 essential 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 FEM
-support~\cite{KnepleyKarpeev09,SieveFEMPreprint}. However, in order to support faults and fault friction models, we must
-expand this range of operations.
+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 FEM support
+\cite{Knepley:Karpeev:2009,SieveFEMPreprint}. However, in order to
+support faults and fault friction models, we must expand this range of
+operations.
-In order to impose a given fault slip, as in equation~(\ref{eq:faultDisp}), we must at least be able to represent two
-displacements for any vertex on the fault. Thus, we could designate ``fault vertices'' which would posses twice as many
-displacement degrees of freedom, as is done in many codes. However, this would mean that global variable offsets would
-have to be stored by cell, rather than vertex, significantly increasing storage costs. Moreover, this scheme would not
-be capable of supporting the cohesive cell formulation~\cite{Ortiz} which is used to enforce a fault constitutive
-model. Therefore, we choose to modify the initial finite element mesh by replacing each fault face with a zero-volume
-cohesive cell, as detailed below.
+In our domain decomposition approach, the finite-element mesh includes
+the fault as an interior surface. This forces alignment of the element
+faces along the fault. In order to impose a given fault slip, as in
+equation~(\ref{eq:faultDisp}), we must represent the displacement on
+both sides of the fault for any vertex on the fault. One option is to
+designate ``fault vertices'' which posses twice as many displacement
+degrees of freedom \cite{Aagaard:etal:BSSA:2001}. However, this
+requires storing the global variable indices by cell rather than by
+vertex or adding special fault metadata to the vertices, significantly
+increasing storage costs and/or index lookup costs.
-Given a set of oriented fault faces, we can introduce a set of cohesive cells using a step-by-step modification of the
-Sieve structure representing the mesh. First, for each vertex on the fault which take to be on $S_f^{+}$, we introduce a
-second vertex on $S_f^{-}$, and a third vertex if we are using Lagrange multiplier constraints lying on an edge between
-the first two. The fault faces are organized as a Sieve, and each face has two cells as descendents. Since they are
-consistently oriented, the first cell attached to each face is on the same side $S_f^{+}$ of the fault. We replace the
-vertices on the fault face of each second cell, on $S_f^{-}$, which the newly created vertices. Finally, we add a
-cohesive cell including the original fault face, a face with the newly created vertices, and if necessary the Lagrange
-vertices. These cohesive cells are prisms, so for example, in a tetrahedral mesh the cohesive cells are triangular prisms,
-whereas in a hexahedral meshes, they are again hexahedrons. Lastly, we update the coordinates for the newly introduced
-vertices.
+We choose another option and modify the initial finite-element mesh by
+replacing each fault face with a zero-volume cohesive cell. Many mesh
+generation tools do not support specification of faces on interior
+surfaces, so we construct the set of oriented fault faces from a set
+of vertices marked as lying on the fault. As a preprocessing step, we
+join these vertices into faces, consistently orient them (using a
+common fault normal direction), and associate them with pairs of cells
+in the original mesh.
-Unfortunately, it is uncommon for a set of faces to be specified as the fault surface. Instead, a set of vertices is
-marked as lying on the fault. Thus, as a preprocessing step, we must join these vertices into faces, consistently orient
-them, and connect them to pairs of cells in the original mesh. A more subtle problem arises for faults which do not
-stretch to the domain boundary. For all cells which touch the fault, we must replace the vertices on $S_f^{-}$ with new
-vertices. The set of cells with a face on the fault are easily classified. However, many cells may intersect the fault
-in a single vertex or an edge, which does not allow us to unambiguously determine which side of the fault their lie
-on. Our classification algorithm is given in Alg.~\ref{alg:classifyCells}.
+Given this set of oriented fault faces, we introduce a set of cohesive
+cells using a step-by-step modification of the Sieve data structure
+representing the mesh illustrated in
+Figure~\ref{fig:cohesive:cell}. First, for each vertex on the positive
+side of the fault $S_f^{+}$, we introduce a second vertex on the
+negative side of the fault $S_f^{-}$ and a third vertex corresponding
+to the Lagrange multiplier constraint. The Lagrange multiplier vertex
+lies on an edge between the vertex on $S_{f^+}$ and the vertex on
+$S_{f^-}$. The fault faces are organized as a Sieve (\matt{Add
+ description here}), and each face has the two cells it is associated
+with as descendents. Because the cells are consistently oriented, the
+first cell attached to each face is on the positive side of the fault,
+i.e., $S_f^{+}$. We replace the vertices on the fault face of each
+second cell, which is on the negative side of the fault, i.e.,
+$S_f^{-}$, with the newly created vertices. Finally, we add a cohesive
+cell including the original fault face, a face with the newly created
+vertices, and the Lagrange vertices. These cohesive cells are
+prisms. For example, in a tetrahedral mesh the cohesive cells are
+triangular prisms, whereas in a hexahedral meshes, they are again
+hexahedrons.
-In order to classify cells, we begin by iterating over the set of fault vertices. For each vertex, we examine the set of
-cells attached to that vertex, called the \emph{support} of the vertex in the Sieve API~\cite{KnepleyKarpeev09}. For
-each unclassified cell in the support, we look at all its neighbors which intersect the fault. If any is classified, we
-give the cell this classification. If not, we continue with a breadth first search of neighbors until a classified one
-is found. This search must terminate since there are a finite number of cells surrounding the vertex and at least one is
-classified (which contains a face containing this vertex). This may produce a ``wrap around'' effect on the fault edge,
-depending on the order of the iteration, however this does not impact the correctness of the result.
+
+\begin{figure}[htbp]
+ \centering
+ \brad{TODO: ADD DIAGRAM}
+ %\includegraphics{figs/cohesive_cell}
+ \caption{Construction of cohesive cells for a fault. (a) For each
+ vertex on the fault, introduce a vertex on the negative side of
+ the fault $S_{f^-}$ and a vertex corresponding to the Lagrange
+ multiplier constraint between the pair of vertices on the positive
+ and negative sides of the fault. (b) Replace the vertices in the
+ cell attached to the negative side of the fault with the newly
+ created vertices. (b) Construct cohesive cells with
+ zero volume from the vertices on the positive side of the fault,
+ negative side of the fault, and Lagrange multiplier constraints.}
+ \label{fig:cohesive:cell}
+\end{figure}
+
+
+We must also update the cells on the negative side of the fault that
+touch the fault with only an edge or single vertex. We need to replace
+the original vertices with the newly introduced vertices on the
+negative side of the fault. In cases where the the fault reaches the
+boundaries of the domain, it is relatively easy to identify these
+cells because their other vertices are shared with the cells that have
+on the negative side of the fault. However, in the case of a fault
+that does not reach the boundary of the domain, cells around the ends
+of the fault share their other vertices with cells that have a face on
+the negative side of the fault and cells that have a face on the
+positive side of the fault. We employ a breadth-first classification
+scheme to differentiate the cells so that we can replace the original
+vertices with the newly introduced vertices on the negative side of
+the fault.
+
+In order to classify cells, we iterate over the set of fault
+vertices. For each vertex, we examine the set of cells attached to
+that vertex, called the \emph{support} of the vertex in the Sieve API
+\cite{Knepley:Karpeev:2009}. For each unclassified cell in the
+support, we look at all its neighbors which touch the fault. If any is
+classified, we give the cell this same classification. If not, we
+continue with a breadth-first search of its neighbors until a
+classified cell is found. This search must terminate because there are
+a finite number of cells surrounding the vertex and at least one is
+classified (contains a fault face with this vertex). This may produce
+a ``wrap around'' effect on the fault edge, depending on the order of
+the iteration, however this does not impact the correctness of the
+result.
\matt{Show a picture of the wrap around?}
@@ -656,19 +723,6 @@ depending on the order of the iteration,
%
% We could also extend the fault with a set of ``halo faces'' to divide one side from the other
-%\begin{algorithm}\label{alg:classifyCells}
-\begin{verbatim}
-for each fault vertex:
- for each cell in the vertex support:
- if classified:
- continue
- else:
- for each nCell in a breadth-first search of neighbors:
- if nCell intersects the fault:
- if nCell is classified:
- give cell the same classification as nCell
-\end{verbatim}
-%\end{algorithm}
\matt{What about the bug we had with fault ends?}
@@ -677,17 +731,26 @@ for each fault vertex:
\subsection{Quasi-static Simulations}
-The introduction of Lagrange multipliers to implement the constraints along a fault produce the saddle-point problem
-shown in equation~(\ref{eqn:saddle-point}). Traditional black-box parallel preconditioners, such as the Additive Schwarz
-Method (ASM)~\cite{SmithBjorstadGropp96}, are not very effective for this type of problem and produce very slow
-convergence. However, PETSc provides tools to construct many variations of effective parallel preconditioners for saddle
-point problems.
+\matt{TODO: Add a paragraph on solving linear system using Krylov
+ methods. Define role of preconditioner.}
-The PCFIELDSPLIT~\cite{petsc-user-ref} preconditioner in PETSc allows the user to define sets of unknowns which
-correspond to different fields in the physical problem. This scheme is flexible enough to accomodate an arbitrary number
-of fields, mixed discretizations, fields defined over a subset of the mesh, etc. Once these fields are defined, a
-substantial range of preconditioners can be assembled using only PyLith options for
-PETSc. Figure~(\ref{fig:solver-options}) shows example preconditioners and the options necessary to construct them.
+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) \cite{Smith:etal:1996}, are not very effective for this
+type of problem and produce very slow convergence. However, PETSc
+provides tools to construct many variations of effective parallel
+preconditioners for saddle point problems.
+
+The PCFIELDSPLIT \cite{PETSc:manual} preconditioner in PETSc allows
+the user to define sets of unknowns which correspond to different
+fields in the physical problem. This scheme is flexible enough to
+accomodate an arbitrary number of fields, mixed discretizations,
+fields defined over a subset of the mesh, etc. Once these fields are
+defined, a substantial range of preconditioners can be assembled using
+only PyLith options for PETSc. Figure~(\ref{fig:solver-options}) shows
+example preconditioners and the options necessary to construct them.
\begin{figure}
\centering
@@ -722,13 +785,15 @@ PETSc. Figure~(\ref{fig:solver-options})
\label{fig:solver-options}
\end{figure}
-\matt{Describe custom preconditioner, and show equation for its construction (from manual I think).}
+\matt{Describe custom preconditioner, and show equation for its
+ construction (from manual I think).}
-\matt{Table of convergence for a fault example with ASM, block tridiagonal, block tridiagonal+custom preconditioner, LSC}
+\matt{Table of convergence for a fault example with ASM, block
+ tridiagonal, block tridiagonal+custom preconditioner, LSC}
\subsection{Dynamic Simulations}
-In the dynamic simulations the the Courant-Friderichs-Lewy condition
+In the dynamic simulations 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
@@ -744,7 +809,7 @@ finite-element basis functions in these
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,
+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
@@ -821,7 +886,7 @@ residual assuming the increment in the s
\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
+$d\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,
@@ -868,9 +933,9 @@ Letting
\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,
+and recognizing that $\tensor{S}_p$ is diagonal because $\tensor{K}$
+and $\tensor{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
@@ -896,8 +961,8 @@ in spontaneous rupture models. Because $
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
+(equation~(\ref{eqn:spontaneous:rupture:slip:update})) are exact and
+simplify to
\begin{equation}
\partial \vec{d} = - \tensor{R}^{-1} \cdot
\left( \tensor{K}_{n^+n^+}^{-1} + \tensor{K}_{n^-n^-}^{-1} \right)
diff -r 26022948ab53 -r 20748e3ccedf references.bib
--- a/references.bib Thu Sep 01 20:33:15 2011 -0500
+++ b/references.bib Tue Sep 06 13:43:04 2011 -0700
@@ -613,12 +613,12 @@
be far more aperiodic than has been suggested. },
}
- at book{SmithBjorstadGropp96,
- author = {Barry F. Smith and Petter Bj{\o}rstad and William D. Gropp},
+ at Book{Smith:etal:1996,
+ author = {Smith, B.~F. and Bj{\o}rstad, P. and Gropp, W.~D.},
title = {Domain Decomposition: Parallel Multilevel Methods for Elliptic Partial Differential Equations},
- url = {http://www.mcs.anl.gov/\~{ }bsmith/ddbook.html},
publisher = {Cambridge University Press},
year = {1996}
+ url = {http://www.mcs.anl.gov/\~{ }bsmith/ddbook.html},
}
@Article{Liu:etal:2006,
@@ -635,16 +635,6 @@
doi = {10.1785/0120060036},
}
- at 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},
- title = {{PETS}c Users Manual},
- institution = {Argonne National Laboratory},
- year = 2010,
- number = {ANL-95/11 - Revision 3.1},
- url = {http://www.mcs.anl.gov/petsc}
-}
-
@article{Brune:1970,
author = {Brune, James N.},
title = {Tectonic stress and spectra of seismic shear waves
@@ -654,14 +644,6 @@
volume = 75,
month = sep # {~10},
pages = {4997--5009},
-}
-
- at 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}
}
@article{Komatitsch:Vilotte:1998,
@@ -725,7 +707,7 @@
seismic risk assessment.},
}
- at article{KnepleyKarpeev09,
+ at article{Knepley:Karpeev:2009,
author = {Matthew G. Knepley and Dmitry A. Karpeev},
title = {Mesh Algorithms for {PDE} with {Sieve} {I}: {Mesh} Distribution},
journal = {Scientific Programming},
@@ -761,6 +743,7 @@
Number = "ANL-95/11 - Revision 3.1",
Institution = "Argonne National Laboratory",
Year = 2010,
+ url = {http://www.mcs.anl.gov/petsc}
}
@InProceedings{PETSc:efficient,
More information about the CIG-COMMITS
mailing list