[cig-commits] commit: More work on solver stuff. Added outline for benchmarks.
Mercurial
hg at geodynamics.org
Tue Mar 27 11:24:37 PDT 2012
changeset: 82:197e088cf409
tag: tip
user: Brad Aagaard <baagaard at usgs.gov>
date: Tue Mar 27 11:24:23 2012 -0700
files: faultRup.tex
description:
More work on solver stuff. Added outline for benchmarks.
diff -r 020cfc86ea0a -r 197e088cf409 faultRup.tex
--- a/faultRup.tex Mon Mar 26 17:34:10 2012 -0700
+++ b/faultRup.tex Tue Mar 27 11:24:23 2012 -0700
@@ -277,8 +277,8 @@ the Lagrange multipliers (fault traction
\begin{split}
- \int_{V} \nabla\mathbf{\phi} : \mathsf{\sigma} \, dV
+ \int_{S_T} \mathbf{\phi} \cdot \mathbf{T} \, dS
- - \int_{S_f^{+}} \mathbf{\phi} \cdot \mathbf{l} \, dS
- + \int_{S_f^{-}} \mathbf{\phi} \cdot \mathbf{l} \, dS \\
+ - \int_{S_{f^+}} \mathbf{\phi} \cdot \mathbf{l} \, dS
+ + \int_{S_{f^-}} \mathbf{\phi} \cdot \mathbf{l} \, dS \\
+ \int_{V} \mathbf{\phi} \cdot \mathbf{f} \, dV
- \int_{V} \mathbf{\phi} \cdot \rho \frac{\partial^{2}\mathbf{u}}{\partial t^{2}} \, dV
=0.
@@ -325,8 +325,8 @@ for all $\mathbf{a}_m$, which leads to
\begin{split}
- \int_{V} \nabla \mathbf{N}_m^T \cdot \mathsf{\sigma} \, dV
+ \int_{S_T} \mathbf{N}_m^T \cdot \mathbf{T} \, dS
-- \int_{S_f^{+}} \mathbf{N}_m^T \cdot \mathbf{N}_p \cdot \mathbf{l}_p \, dS \\
-+ \int_{S_f^{-}} \mathbf{N}_m^T \cdot \mathbf{N}_p \cdot \mathbf{l}_p \, dS
+- \int_{S_{f^+}} \mathbf{N}_m^T \cdot \mathbf{N}_p \cdot \mathbf{l}_p \, dS \\
++ \int_{S_{f^-}} \mathbf{N}_m^T \cdot \mathbf{N}_p \cdot \mathbf{l}_p \, dS
+ \int_{V} \mathbf{N}_m^T \cdot \mathbf{f} \, dV \\
- \int_{V} \rho \mathbf{N}_m^T \cdot \mathbf{N}_n \cdot \frac{\partial^2 \mathbf{u}_n}{\partial
t^2} \, dV
@@ -352,9 +352,12 @@ equations in terms of the increment in t
equations in terms of the increment in the solution from time $t$ to
$t+\Delta t$ rather than the solution at time $t+\Delta t$.
Consequently, rather than constructing a system with the form
-$\mathbf{A} \cdot \mathbf{u}(t+\Delta t) = \mathbf{b}(t+\Delta t)$, we construct a system
-with the form $\mathbf{A} \cdot \mathbf{du} = \mathbf{b}(t+\Delta t) - \mathbf{A}
-\cdot \mathbf{u}(t)$, where $\mathbf{u}(t+\Delta t) = \mathbf{u}(t) + \mathbf{du}$.
+$\mathbf{A} \cdot \mathbf{u}(t+\Delta t) = \mathbf{b}(t+\Delta t)$, we
+construct a system with the form $\mathbf{A} \cdot \mathbf{du} =
+\mathbf{b}(t+\Delta t) - \mathbf{A} \cdot \mathbf{u}(t)$, where
+$\mathbf{u}(t+\Delta t) = \mathbf{u}(t) + \mathbf{du}$.
+
+\brad{Add description of nondimensionalization?}
% ------------------------------------------------------------------
\subsection{Quasi-static Simulations}
@@ -368,8 +371,8 @@ and the loading conditions. Considering
\begin{split}
- \int_{V} \nabla \mathbf{N}_m^T \cdot \mathsf{\sigma}(t+\Delta t) \, dV
+ \int_{S_T} \mathbf{N}_m^T \cdot \mathbf{T}(t+\Delta t) \, dS \\
- - \int_{S_f^{+}} \mathbf{N}_m^T \cdot \mathbf{N}_p \cdot \mathbf{l}_p(t+\Delta t) \, dS \\
- + \int_{S_f^{-}} \mathbf{N}_m^T \cdot \mathbf{N}_p \cdot \mathbf{l}_p(t+\Delta t) \, dS
+ - \int_{S_{f^+}} \mathbf{N}_m^T \cdot \mathbf{N}_p \cdot \mathbf{l}_p(t+\Delta t) \, dS \\
+ + \int_{S_{f^-}} \mathbf{N}_m^T \cdot \mathbf{N}_p \cdot \mathbf{l}_p(t+\Delta t) \, dS
\\
+ \int_{V} \mathbf{N}_m^T \cdot \mathbf{f}(t+\Delta t) \, dV
=\mathbf{0},
@@ -432,12 +435,37 @@ associated with the constraints, equatio
\mathbf{L} = \int_{S_f} \mathbf{N}_p^T \cdot (\mathbf{N}_{n^+} - \mathbf{N}_{n^-}) \, dS.
\end{equation}
Thus, the Jacobian of the entire system has the form,
-\begin{equation}\label{eqn:saddle-point}
+\begin{equation}\label{eqn:saddle:point}
\mathbf{A} =
\left( \begin{array}{cc}
\mathbf{K} & \mathbf{L}^T \\ \mathbf{L} & \mathbf{0}
\end{array} \right).
\end{equation}
+Note that the terms in $\mathbf{N_{n^+}}$ and $\mathbf{N_{n^-}}$ are
+identical, but they refer to degrees of freedom on the positive and
+negative sides of the fault, respectively. Consequently, in practice
+we compute the terms for the positive side of the fault and assemble
+the terms into the appropriate degrees of freedom for both sides of
+the fault. Hence, we compute
+\begin{equation}
+ \label{eqn:jacobian:constraint:code}
+ \mathbf{L_p} = \int_{S_f} \mathbf{N}_p^T \cdot \mathbf{N}_{n^+} \, dS,
+\end{equation}
+with the Jacobian of the entire system taking the form,
+\begin{equation}\label{eqn:saddle:point:code}
+ \mathbf{A} =
+ \left( \begin{array}{cccc}
+ \mathbf{K}_{nn} & \mathbf{K}_{nn^+} & \mathbf{K}_{nn^-} & \mathbf{0} \\
+ \mathbf{K}_{n^+n} & \mathbf{K}_{n^+n^+} & \mathbf{0} & \mathbf{L}_p^T \\
+ \mathbf{K}_{n^-n} & \mathbf{0} & \mathbf{K}_{n^-n^-} & -\mathbf{L}_p^T \\
+ \mathbf{0} & \mathbf{L}_p & -\mathbf{L}_p & \mathbf{0}
+ \end{array} \right),
+\end{equation}
+where $n$ denotes degrees of freedom not associated with the fault,
+$n^-$ denotes degrees of freedom associated with the negative side of
+the fault, $n^+$ denotes degrees of freedom associated with the
+positive side of the fault, and $p$ denotes degrees of freedom
+associated with the Lagrange multipliers.
% ------------------------------------------------------------------
\subsection{Dynamic Simulations}
@@ -453,8 +481,8 @@ inertial term in equation~\ref{eqn:quasi
\begin{split}
- \int_{V} \nabla \mathbf{N}_m^T \cdot \mathsf{\sigma}(t+\Delta t) \, dV
+ \int_{S_T} \mathbf{N}_m^T \cdot \mathbf{T}(t+\Delta t) \, dS \\
- - \int_{S_f^{+}} \mathbf{N}_m^T \cdot \mathbf{N}_p \cdot \mathbf{l}_p(t+\Delta t) \, dS \\
- + \int_{S_f^{-}} \mathbf{N}_m^T \cdot \mathbf{N}_p \cdot \mathbf{l}_p(t+\Delta t) \, dS
+ - \int_{S_{f^+}} \mathbf{N}_m^T \cdot \mathbf{N}_p \cdot \mathbf{l}_p(t+\Delta t) \, dS \\
+ + \int_{S_{f^-}} \mathbf{N}_m^T \cdot \mathbf{N}_p \cdot \mathbf{l}_p(t+\Delta t) \, dS
\\
+ \int_{V} \mathbf{N}_m^T \cdot \mathbf{f}(t+\Delta t) \, dV \\
- \int_{V} \rho \mathbf{N}_m^T \cdot \mathbf{N}_n \cdot
@@ -524,12 +552,13 @@ do the Lagrange multipliers, which are a
do the Lagrange multipliers, which are available from the current
trial 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 slip. PyLith includes a variety of slip-time histories, including
-a step function, a linear ramp (constant slip rate), the integral of
-Brune's far-field time function \citep{Brune:1970}, a sine-cosine
-function developed by \cite{Liu:etal:2006}, and a user-defined time
-history. These are discussed in detail in the PyLith manual
+are the \textit{change} in the tractions on the fault surfaces
+corresponding to the slip. PyLith includes a variety of slip-time
+histories, including a step function, a linear ramp (constant slip
+rate), the integral of Brune's far-field time function
+\citep{Brune:1970}, a sine-cosine function developed by
+\cite{Liu:etal:2006}, and a user-defined time history. These are
+discussed in detail in the PyLith manual
\citep{PyLith:manual:1.6.2}. PyLith allows specification of the slip
initiation time independently at each location as well as
superposition of multiple earthquake ruptures with different origin
@@ -541,24 +570,25 @@ times, thereby permitting complex slip b
In a spontaneous (dynamic) fault rupture a constitutive model controls
the tractions on the fault surface. The fault slip evolves based on
-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
+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 find
-both the fault slip and Lagrange multipliers that are consistent with
-the fault constitutive model.
+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
-tractions required to keep the fault locked. If the Lagrange
-multipliers exceed the fault tractions allowed by the fault
-constitutive model, then we iterate to find the increment in slip that
-yields Lagrange multipliers that satisfy the fault constitutive model.
-On the other hand, if the Lagrange multipliers do not exceed the fault
-tractions allowed by the fault constitutive model, then the increment
-in fault slip remains zero, and no adjustements to the solution are
-necessary.
+tractions required to lock the fault. If the Lagrange multipliers
+exceed the fault tractions allowed by the fault constitutive model,
+then we iterate to find the increment in slip that yields Lagrange
+multipliers that satisfy the fault constitutive model. On the other
+hand, if the Lagrange multipliers do not exceed the fault tractions
+allowed by the fault constitutive model, then the increment in fault
+slip remains zero, and no adjustements to the solution are necessary.
In iterating to find the fault slip and Lagrange multipliers that
satisfy the fault constitutive model, we employ the following
@@ -566,12 +596,13 @@ multipliers necessary to satisfy the fau
multipliers necessary to satisfy the fault constitutive model for the
current estimate of slip. We then compute the increment in fault slip
corresponding to this perturbation in the Lagrange multipliers
-assuming the rest of the solution is unchanged. That is, we consider
-only the degrees of freedom associated with the fault interface when
-computing how a perturbation in the Lagrange multipliers corresponds
-to a change in fault slip. In terms of the general form of a linear
-system of equations ($\mathbf{A} \mathbf{u} = \mathbf{b}$), our subset
-of equations has the form
+assuming deformation is limited to vertices on the fault. That is, we
+consider only the degrees of freedom associated with the fault
+interface when computing how a perturbation in the Lagrange
+multipliers corresponds to a change in fault slip. In terms of the
+general form of a linear system of equations ($\mathbf{A} \mathbf{u} =
+\mathbf{b}$), our subset of equations based on
+equation~(\ref{eqn:saddle:point:code}) has the form
\begin{equation}
\begin{pmatrix}
\mathbf{K}_{n^+n^+} & 0 & \mathbf{L}_p^T \\
@@ -605,18 +636,17 @@ multipliers $\partial \mathbf{l}_p$:
\mathbf{L}_p^T \cdot \partial \mathbf{l}_p, \\
\partial \mathbf{d}_p = \partial \mathbf{u}_{n^+} - \partial \mathbf{u}_{n^-}.
\end{gather}
-This estimate is exact if all degrees of freedom not adjacent to the
-fault are constrained (which is generally only the case for simple toy
-problems), quite good if the deformation associated with the fault
-slip is confined close to the fault (which occurs in most simulations
-because the 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). The number of iterations required
-to satisfy the fault constitutive model decreases as the estimate of
-the fault slip becomes more accurate. We will discuss the procedures
-for solving these equations for both quasi-static and dynamic
-simulations in section~\ref{sec:??}.
+This estimate is exact for constant friction if all degrees of freedom
+not adjacent to the fault are constrained (which is generally only the
+case for simple toy problems), quite good if the deformation
+associated with the fault slip is confined close to the fault (which
+occurs in most simulations due to a finite-extent of slip in the
+middle of the domain), and poor if the deformation associated with
+fault slip extends to the boundaries (which is relatively rare). The
+number of iterations required to satisfy the fault constitutive model
+decreases as the estimate of the fault slip becomes more accurate. We
+will discuss the procedures for solving these equations for both
+quasi-static and dynamic simulations in section~\ref{sec:??}.
PyLith includes several commonly used fault constitutive models, all
of which specify the shear traction on the fault $T_f$ as a function
@@ -647,7 +677,7 @@ particular the Sieve package for finite-
The Sieve application programming interface (API) for mesh
representation and manipulation is based upon a direct acyclic graph
-(DAG) representation of the \textit{covering} relation in a mesh. For
+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 relation, the interface can
be both concise and quite general. Using this generic API, PyLith is
@@ -681,58 +711,62 @@ Given this set of oriented fault faces,
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
+Figure~\ref{fig:cohesive:cell}. First, for each vertex on the negative
+side of the fault $S_{f^-}$, we introduce a second vertex on the
+positive 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, 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\brad{CHECK THIS} side of the fault, i.e., $S_f^{+}$. We
+the negative 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
+on the positive 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 hexahedrons.
-We must also update all cells on the negative side of the fault that
+We must also update all cells on the positive 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
+positive 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
-faces on the negative side of the fault. However, in the case of a
-fault that does not reach the boundary of the domain, cells near 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 between these two cases, so
-that we can replace the original vertices with the newly introduced
-vertices on the negative side of the fault.
+cells because these vertices are shared with the cells that have faces
+on the positive side of the fault. However, in the case of a fault
+that does not reach the boundary of the domain, cells near the ends of
+the fault share vertices with cells that have a face on the positive
+side of the fault \textit{and} cells that have a face on the negative
+side of the fault. We use a breadth-first classification scheme to
+classify all cells with vertices on the fault into those having vertices
+on the positive side of the fault and those having vertices on the
+negative side of the fault, so that we can replace the original
+vertices with the newly introduced vertices on the positive side of
+the fault.
-In classifying the cells, we iterate over the set of fault
+In classifying the 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
\citep{Knepley:Karpeev:2009}. For each unclassified cell in the
-support, we look at all its neighbors that touch the fault. If any is
-classified, we give the cell this same classification. If not, we
+support, we look at all of its neighbors that 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.
+classified (contains a face on the fault with this vertex). Depending
+on the order of the iteration, this can produce a ``wrap around''
+effect at the ends of the fault, but it does not affect the numerical
+solution because the fault slip should be zero at the edges of the
+fault.
-\matt{Show a picture of the wrap around?}
-% We could instead do a breadth-first classification of cells starting from the cells with a face on the fault, and
-% proceedings in sets of neighbors. This would ``meet halfway'' around fault edges.
+% We could instead do a breadth-first classification of cells starting
+% from the cells with a face on the fault, and proceedings in sets of
+% neighbors. This would ``meet halfway'' around fault edges.
%
-% We could also extend the fault with a set of ``halo faces'' to divide one side from the other
+% We could also extend the fault with a set of ``halo faces'' to
+% divide one side from the other
% ------------------------------------------------------------------
@@ -741,23 +775,24 @@ result.
\subsection{Quasi-static Simulations}
In order to solve the large, sparse systems of linear equations
-arising in quasi-static simulations, we employ preconditioned Krylov
-subspace methods~\citep{Saad03}. We create a sequence of vectors by
-repeatedly applying the system matrix to the rhs vector, $\{ A^k b\}$,
-and they form a basis for a subspace, termed the Krylov space. We can
-efficiently find an approximate solution in this subspace.
-Since sparse matrix-vector multiplication is very scalable, this is
-the method of choice for parallel simulation. However, for most
+arising in our quasi-static simulations, we employ preconditioned
+Krylov subspace methods~\citep{Saad03}. We create a sequence of
+vectors by repeatedly applying the system matrix to the
+right-hand-side vector, $\mathbf{A^k} \cdot \mathbf{b}$, and they form
+a basis for a subspace, termed the Krylov space. We can efficiently
+find an approximate solution in this subspace. Because sparse
+matrix-vector multiplication is scalable via parallel processing, this
+is the method of choice for parallel simulation. However, for most
physically relevant problems, the Krylov solver requires a
preconditioner in order to accelerate convergence. While generic
preconditioners exist~\citep{Saad03,Smith:etal:1996}, the method must
-often be specialized to a particular problem. In this section we describe a
-preconditioner specialized to the fault interface problem and Lagrange
-multipliers.
+often be specialized to a particular problem. In this section we
+describe a preconditioner specialized to our formulation for fault
+slip with 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
+equation~(\ref{eqn:saddle:point}). Traditional black-box parallel
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
@@ -770,7 +805,7 @@ accomodate an arbitrary number of fields
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
+only PyLith options for PETSc. Table~\ref{tab:solver:options} shows
example preconditioners and the options necessary to construct them.
Another option involves using using the field split preconditioner in
@@ -784,17 +819,19 @@ matrix. Our system Jacobian has the form
\mathbf{L} & \mathbf{0}
\end{array} \right).
\end{equation}
-We use the Schur complement of block $\mathbf{K}$ to examine the form of $\mathbf{A}^{-1}$,
-\begin{equation}
+We use the Schur complement of block $\mathbf{K}$ to examine the form
+of $\mathbf{A}^{-1}$,
+\begin{gather}
\mathbf{A}^{-1} = \left( \begin{array}{cc}
- \mathbf{K}^{-1} + \mathbf{K}^{-1} \mathbf{L}^{T}
- (-\mathbf{L} \mathbf{K}^{-1} \mathbf{L}^{T})^{-1} \mathbf{L} \mathbf{K}^{-1} &
- -\mathbf{K}^{-1} \mathbf{L}^{T }(-\mathbf{L} \mathbf{K}^{-1} \mathbf{L}^{T})^{-1} \\
- -(-\mathbf{L} \mathbf{K}^{-1} \mathbf{L}^{T})^{-1} \mathbf{L}
- \mathbf{K}^{-1} &
- -(\mathbf{L} \mathbf{K}^{-1} \mathbf{L}^T)^{-1}
- \end{array} \right),
-\end{equation}
+ \mathbf{A}^{-1}_{nn} & \mathbf{A}^{-1}_{np} \\
+ \mathbf{A}^{-1}_{pn} & \mathbf{A}^{-1}_{pp}
+ \end{array} \right), \\
+ \mathbf{A}^{-1}_{nn} = \mathbf{K}^{-1} + \mathbf{K}^{-1} \mathbf{L}^{T}
+ (-\mathbf{L} \mathbf{K}^{-1} \mathbf{L}^{T})^{-1} \mathbf{L} \mathbf{K}^{-1}, \\
+ \mathbf{A}^{-1}_{np} = -\mathbf{K}^{-1} \mathbf{L}^{T }(-\mathbf{L} \mathbf{K}^{-1} \mathbf{L}^{T})^{-1}, \\
+ \mathbf{A}^{-1}_{pn} = -(-\mathbf{L} \mathbf{K}^{-1} \mathbf{L}^{T})^{-1} \mathbf{L} \mathbf{K}^{-1}, \\
+ \mathbf{A}^{-1}_{pp} = -(\mathbf{L} \mathbf{K}^{-1} \mathbf{L}^T)^{-1}.
+\end{gather}
A suitable block diagonal approximation of $\mathbf{A}^{-1}$ is
\begin{equation}
\mathbf{P}^{-1} = \left( \begin{array}{cc}
@@ -814,23 +851,30 @@ which leads to
\end{array} \right).
\end{equation}
-The elastic block $K$ can be further split into blocks associated with displacements along the $x$, $y$, and $z$
-axes. It is known that the vector Laplacian is spectrally equivalent to this operator~\citep{AskMarkAdams}, and each
-component is efficiently preconditioned by Algebraic Multigrid (AMG), such as ML~\citep{ML}. Algebraic Multigrid mimics
-the action of traditional geometric multgrid, but it generates coarse level operators and interpolation matrices using
-only the system matrix, treated as a weighted graph, rather than a separate description of the problem geometry, such as
-a mesh. We use PCFIELDSPLIT to split the elastic block and separately apply AMG to each component.
+The elastic block $K$ can be further split into blocks associated with
+displacements along the $x$, $y$, and $z$ axes. It is known that the
+vector Laplacian is spectrally equivalent to this
+operator~\citep{AskMarkAdams}\matt{Ask Mark Adams}, and each component
+is efficiently preconditioned by Algebraic Multigrid (AMG), such as
+the ML library \citep{ML}. AMG mimics the action of traditional
+geometric multgrid, but it generates coarse level operators and
+interpolation matrices using only the system matrix, treated as a
+weighted graph, rather than a separate description of the problem
+geometry, such as a mesh. We use PCFIELDSPLIT to split the elastic
+block and separately apply AMG to each component.
-We will now focus on evaluating the lower portion of the preconditioning matrix
-associated with the Lagrange multipliers since stock PETSc
-preconditioners can handle the upper portion. We approximate
-$\mathbf{K}^{-1}$ with the inverse of the diagonal portion of
-$\mathbf{K}$. $\mathbf{L}$ consists of integrating the products of
-basis functions over the fault faces. Its structure depends on the
-quadrature scheme and the choice of basis functions. For conventional
-low order finite-elements and Gauss quadrature, $\mathbf{L}$ contains
-nonzero terms coupling the degree of freedom for each coordinate axes
-with the corresponding degree of freedom of the other vertices in a
+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
+discussed in the previous paragraph. In computing
+$\mathbf{P_\mathit{fault}}$ we we approximate $\mathbf{K}^{-1}$ with
+the inverse of the diagonal portion of $\mathbf{K}$. $\mathbf{L}$
+consists of integrating the products of basis functions over the fault
+faces. Its structure depends on the quadrature scheme and the choice
+of basis functions. For conventional low order finite-elements and
+Gauss quadrature, $\mathbf{L}$ contains nonzero terms coupling the
+degree of freedom for each coordinate axes of a vertex with the
+corresponding degree of freedom of the other vertices in a
cell. However, if we collocate quadrature points at the cell vertices,
then only one basis function is nonzero at each quadrature point and
$\mathbf{L}$ becomes block diagonal; this is also true for spectral
@@ -838,46 +882,35 @@ quadrature points. This leads to a diago
quadrature points. This leads to a diagonal matrix for the lower
portion of the conditioning matrix,
\begin{equation}
- \mathbf{P}_f = -\mathbf{L_p} (\mathbf{K}_{n+n+} + \mathbf{K}_{n-n-}) \mathbf{K}_p^{T}.
+ \mathbf{P}_\mathit{fault} = -\mathbf{L}_p (\mathbf{K}_{n+n+} + \mathbf{K}_{n-n-}) \mathbf{L}_p^{T},
\end{equation}
+where $\mathbf{L}_p$ is given in equation~(\ref{eqn:jacobian:constraint:code}) and
+$\mathbf{K}_{n+n+}$ and $\mathbf{K}_{n-n-}$ are the diagonal terms
+from equation~(\ref{eqn:saddle:point:code}).
-\brad{Matt: I conjecture that collocation, because it makes $\mathbf{L}$ block diagonal, is more tolerant of the
- diagonal approximation for $\mathbf{K}$.}
+% Matt conjectures that collocation, because it makes $\mathbf{L}$
+% block diagonal, is more tolerant of the diagonal approximation for
+% $\mathbf{K}$.}
For the upper portion of the preconditioning matrix associated with
-elasticity, we have found algrebraic multigrid preconditioners
-provide substantially faster convergence that the Additive Schwarz
-method. We combine the field split preconditioner with the algebraic multigrid
-preconditioner, such that we precondition the degrees of
-freedom for each global coordinate axis independently. Table~\ref{tab:iterates}
-shows the number of iterates required to solve a problem with three faults to a
-relative tolerance of $10^{-8}$. It clearly shows the surperiority of our custom
-fault preconditioner. It also reveals that while is scales fairly well, the custom
-preconditioner does have a dependence on problem size.
-
-\begin{table}
-\centering
-\begin{tabular}{llrr}
-KSP & PC & P1 Iterates & P2 Iterates \\ % Use refine.cfg for P2
-\hline
-GMRES & ASM/Shifted ICC & 77 & $>500$ \\ % asm.cfg
-GMRES & FS/ML & $>500$ & $>500$ \\ % fieldsplit.cfg
-GMRES & FS/Schur/ML & 62 & 135 \\ % schur.cfg
-GMRES & FS/LSC/ML & 53 & 108 \\ % lsc.cfg
-GMRES & FS/ML/Custom & 37 & 49 \\ % fieldsplit.cfg custompc.cfg
-\end{tabular}
-\caption{Performance of Krylov solvers on a problem with 3 faults. Problem 1 (P1) has 25000 unknowns and 645
- constraints, whereas the larger Problem 2 (P2) has 191,000 unknowns and 2241 constraints. The relative
- solver tolerance is $10^{-8}$, and the preconditioners used were Additive Schwarz Method (ASM), FieldSplit (FS),
- Incomplete Cholesky (ICC), algebraic multigrid (ML), Schur complement, least-square commutator (LSC), and a custom
- preconditioner for the fault problem.}
-\label{tab:iterates}
-\end{table}
+elasticity, we have found AMG preconditioners provide substantially
+faster convergence that the Additive Schwarz method. We combine the
+field split preconditioner with the AMG preconditioner, such that we
+precondition the degrees of freedom for each global coordinate axis
+independently. Table~\ref{tab:iterates} shows the number of iterates
+required to solve a problem with prescribed slip on three faults
+\brad{Setup solver test for weak scaling computations or use
+ strike-slip quasi-static benchmark? Add figure describing
+ problem. Add problm description to PyLith manual and cite manual for
+ details?} to a relative tolerance of $10^{-8}$. It clearly shows
+the surperiority of our custom fault preconditioner. It also reveals
+that while is scales fairly well, the custom preconditioner does have
+a dependence on problem size.
\subsection{Dynamic Simulations}
-In the dynamic simulations the Courant-Friderichs-Lewy condition
-controls the stability of the time integration. In most dynamic
+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
@@ -959,14 +992,14 @@ associated with the fault slip constrain
associated with the fault slip constraint in solving the system of
equations via a Schur's complement algorithm. We compute an initial
residual assuming the increment in the solution is zero (i.e.,
-$\mathbf{du}_n = \mathbf{0}$ and $d\mathbf{l}_p = \mathbf{0}$,
+$\mathbf{du}_n = \mathbf{0}$ and $\mathbf{dl}_p = \mathbf{0}$,
\begin{equation}
\mathbf{r}^* = \begin{pmatrix} \mathbf{r}_n^* \\ \mathbf{r}_p^* \end{pmatrix} =
\begin{pmatrix} \mathbf{b}_n \\ \mathbf{b}_p \end{pmatrix}
- \begin{pmatrix}
\mathbf{K} & \mathbf{L}^T \\ \mathbf{L} & 0
\end{pmatrix}
- \begin{pmatrix} \mathbf{u}_n \\ \mathbf{l}_n \end{pmatrix}.
+ \begin{pmatrix} \mathbf{u}_n \\ \mathbf{l}_p \end{pmatrix}.
\end{equation}
We compute a corresponding initial solution to the system of equations
$\mathbf{du}_n^*$ ignoring the off-diagonal blocks in the Jacobian and
@@ -988,7 +1021,7 @@ residual is
- \begin{pmatrix}
\mathbf{K} & \mathbf{L}^T \\ \mathbf{L} & 0
\end{pmatrix}
- \begin{pmatrix} \mathbf{du}_n \\ d\mathbf{l}_n \end{pmatrix}.
+ \begin{pmatrix} \mathbf{du}_n \\ \mathbf{dl}_p \end{pmatrix}.
\end{equation}
Solving the first row of equation~(\ref{eqn:lumped:jacobian:residual})
for the increment in the solution and accounting for the structure of
@@ -996,9 +1029,9 @@ each side of the fault, we have
each side of the fault, we have
\begin{gather}
\mathbf{du}_{n^+} =
- \mathbf{du}_{n^+}^* - \mathbf{K}_{n^+n^+}^{-1} \cdot \mathbf{L}_p^T \cdot d\mathbf{l}_p, \\
+ \mathbf{du}_{n^+}^* - \mathbf{K}_{n^+n^+}^{-1} \cdot \mathbf{L}_p^T \cdot \mathbf{dl}_p, \\
\mathbf{du}_{n^-} =
- \mathbf{du}_{n^-}^* + \mathbf{K}_{n^-n^-}^{-1} \cdot \mathbf{L}_p^T \cdot d\mathbf{l}_p.
+ \mathbf{du}_{n^-}^* + \mathbf{K}_{n^-n^-}^{-1} \cdot \mathbf{L}_p^T \cdot \mathbf{dl}_p.
\end{gather}
Substituting into the second row of
equation~(\ref{eqn:lumped:jacobian:residual}) and isolating the term
@@ -1006,7 +1039,7 @@ with the increment in the Lagrange multi
\begin{equation}
\mathbf{L}_p \cdot
\left( \mathbf{K}_{n^+n^+}^{-1} + \mathbf{K}_{n^-n^-}^{-1} \right) \cdot
- \mathbf{L}^T_p \cdot d\mathbf{l}_p =
+ \mathbf{L}^T_p \cdot \mathbf{dl}_p =
-\mathbf{r}_p^* + \mathbf{L}_p \cdot
\left( \mathbf{du}_{n^+}^* - \mathbf{du}_{n^-}^* \right).
\end{equation}
@@ -1020,17 +1053,17 @@ and $\mathbf{L}_p$ are diagonal allows u
and $\mathbf{L}_p$ are diagonal allows us to solve for the increment
in the Lagrange multipliers,
\begin{equation}
- d\mathbf{l}_p = \mathbf{S}_p^{-1} \cdot \left(
+ \mathbf{dl}_p = \mathbf{S}_p^{-1} \cdot \left[
-\mathbf{r}_p^* + \mathbf{L}_p \cdot
\left( \mathbf{du}_{n^+}^* - \mathbf{du}_{n^-}^* \right)
- \right).
+ \right].
\end{equation}
Now that we have the increment in the Lagrange multipliers, we can
correct our initial solution $\mathbf{du}_n^*$ so that the true residual
is zero,
\begin{equation}
\mathbf{du}_n =
- \mathbf{du}_n^* - \mathbf{K}^{-1} \cdot \mathbf{L}^T \cdot d\mathbf{l}_p.
+ \mathbf{du}_n^* - \mathbf{K}^{-1} \cdot \mathbf{L}^T \cdot \mathbf{dl}_p.
\end{equation}
Because $\mathbf{K}$ and $\mathbf{L}$ are comprised of diagonal
blocks, this expression for the updates to the solution are local to
@@ -1054,11 +1087,34 @@ simplify to
% ------------------------------------------------------------------
+\section{Performance Benchmark}
+\begin{itemize}
+\item \brad{Solver test or quasi-static strike-slip benchmark?}
+\item Comparison of preconditioner performance
+\item Weak scaling performance for field split with custom preconditioner
+\end{itemize}
+
+% ------------------------------------------------------------------
\section{Code Verification Benchmarks}
\begin{itemize}
- \item Savage and Prescott\charles{Rough draft}
- \item Spontaneous rupture 2-D\brad{Rough draft}
- \item Spontaneous rupture 3-D\brad{Rough draft}
+\item Savage and Prescott
+ \begin{itemize}
+ \item Compare performance of hex8 and tet4
+ \end{itemize}
+\item Spontaneous rupture benchmark: TPV12, TPV13
+ \begin{itemize}
+ \item dipping fault, depth dependent stresses, super-shear rupture,
+ Drucker-Prage elastoplastic model
+ \item Not ideal due to discontinuities in spatial variation of parameters
+ \item 2-D, compare quad4 and tri3
+ \item 3-D, only tet4 (complex geometry)
+ \end{itemize}
+\item Spontaneous rupture 3-D: TPV16
+ \begin{itemize}
+ \item Heterogeneous initial stresses, more consistent nucleation
+ \item No discontinuities in spatial variation of parameters
+ \item compare hex8 and tet4
+ \end{itemize}
\end{itemize}
@@ -1075,6 +1131,9 @@ simplify to
elasticity equation.\\
$\mu_f$ & coefficient of friction.\\
$\mathbf{n}$ & normal vector.\\
+ $\mathbf{P}$ & preconditioning matrix.\\
+ $\mathbf{P}_\mathit{elastic}$ & preconditioning matrix associated with elasticity.\\
+ $\mathbf{P}_\mathit{fault}$ & preconditioning matrix associated with fault slip constraints (Lagrange multipliers).\\
$S_f$ & fault surface.\\
$S_T$ & surface with Neumann boundary conditions.\\
$S_u$ & surface with Dirichlet boundary conditions.\\
@@ -1083,13 +1142,12 @@ simplify to
$T_c$ & scalar shear traction associated with cohesion.\\
$T_f$ & scalar shear traction associated with friction.\\
$T_n$ & scalar normal traction.\\
+ $\mathbf{u}$ & displacement vector.\\
$V$ & spatial domain of model.\\
- $\mathbf{T}$ & traction vector.\\
- $\mathbf{u}$ & displacement vector.\\
+ $\Delta t$ & Time step.\\
$\mathbf{\phi}$ & weighting function.\\
$\rho$ & mass density.\\
- $\mathsf{\sigma}$ & Cauchy stress tensor.\\
-
+ $\mathsf{\sigma}$ & Cauchy stress tensor.
\end{notation}
@@ -1149,11 +1207,11 @@ MGK acknowledges partial support from NS
\begin{table}
\centering
\caption{Example preconditioners for the saddle point problem in
- equation~(\ref{eqn:saddle-point}). All of these field split
+ equation~(\ref{eqn:saddle:point}). All of these field split
preconditioners require the use of the parameters
\texttt{split\_fields = True} and \texttt{matrix\_type = aij} for
\texttt{pylithapp.problem.formulation}.}
-\label{tab:solver-options}
+\label{tab:solver:options}
\begin{tabular}{ll}
$\begin{pmatrix}\hat K & 0 \\ 0 & I\end{pmatrix}$ & $\begin{pmatrix}\hat K & L^T \\ 0 & I\end{pmatrix}$ \\
\texttt{[pylithapp.problem.formulation]} & \texttt{[pylithapp.problem.formulation]} \\
@@ -1200,6 +1258,25 @@ MGK acknowledges partial support from NS
\end{tabular}
\end{table}
+\begin{table}
+\centering
+\begin{tabular}{llrr}
+KSP & PC & P1 Iterates & P2 Iterates \\ % Use refine.cfg for P2
+\hline
+GMRES & ASM/Shifted ICC & 77 & $>500$ \\ % asm.cfg
+GMRES & FS/ML & $>500$ & $>500$ \\ % fieldsplit.cfg
+GMRES & FS/Schur/ML & 62 & 135 \\ % schur.cfg
+GMRES & FS/LSC/ML & 53 & 108 \\ % lsc.cfg
+GMRES & FS/ML/Custom & 37 & 49 \\ % fieldsplit.cfg custompc.cfg
+\end{tabular}
+\caption{Performance of Krylov solvers on a problem with 3 faults. Problem 1 (P1) has 25000 unknowns and 645
+ constraints, whereas the larger Problem 2 (P2) has 191,000 unknowns and 2241 constraints. The relative
+ solver tolerance is $10^{-8}$, and the preconditioners used were Additive Schwarz Method (ASM), FieldSplit (FS),
+ Incomplete Cholesky (ICC), algebraic multigrid (ML), Schur complement, least-square commutator (LSC), and a custom
+ preconditioner for the fault problem.}
+\label{tab:iterates}
+\end{table}
+
% ------------------------------------------------------------------
More information about the CIG-COMMITS
mailing list