[cig-commits] commit: More work on solvers section for dynamic simulations.

[Mercurial]

changeset:   62:32a982452b80
tag:         tip
user:        Brad Aagaard <baagaard at usgs.gov>
date:        Thu Sep 01 10:09:41 2011 -0700
files:       faultRup.tex references.bib
description:
More work on solvers section for dynamic simulations.


diff -r 38890ae93959 -r 32a982452b80 faultRup.tex
--- a/faultRup.tex	Wed Aug 31 16:16:36 2011 -0700
+++ b/faultRup.tex	Thu Sep 01 10:09:41 2011 -0700
@@ -603,9 +603,9 @@ In the dynamic simulations the the Coura
 In the dynamic simulations the the Courant-Friderichs-Lewy condition
 controls the stability of the time integration. In most dynamic
 problems this dictates a relatively small time step so that a typical
-simulation involves tens of thousands of time steps. Thus, we want a
-very efficient solver in order to run a simulation in a reasonable
-amount of time. 
+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
 
@@ -615,45 +615,47 @@ equation~(\ref{eqn:jacobian:constraint})
 equation~(\ref{eqn:jacobian:constraint}). Using conventional
 finite-element basis functions in these integrations results in a
 sparse matrix with off-diagonal terms. Although we can use the same
-solvers as we do for quasi-static simulations, eliminating the
-off-diagonal terms so that the Jacobian is diagonal, results in a much
-faster solver. Not only is the number of operations proportional to
-the number of degrees of freedom, but the memory requirements are
-greatly reduced by storing the diagonal of the matrix using a vector
-rather than a sparse matrix. However, the block structure of our
-Jacobian matrix, with the fault slip constraints occupying
-off-diagonal blocks requires a two step approach to being able to
-solve the linear system of equations without forming a sparse
-matrix.
+solvers as we do for quasi-static simulations to find the solution,
+eliminating the off-diagonal terms so that the Jacobian is diagonal,
+permits use of a much faster solver. With a diagonal Jacobian the
+number of operations required for the solve is proportional to the
+number of degrees of freedom, and the memory requirements are greatly
+reduced by storing the diagonal of the matrix as a vector rather
+than as a sparse matrix. However, the block structure of our Jacobian
+matrix, with the fault slip constraints occupying off-diagonal blocks,
+requires a two step approach in order to solve the linear system
+of equations without forming a sparse matrix.
 
-First, we eliminate the off-diagonal entries in each of the blocks of
-the matrix during the finite-element integrations.  The current best
+First, we eliminate the off-diagonal entries in each block of the
+matrix during the finite-element integrations.  The current best
 available option for eliminating the off-diagonal terms formed during
 the integration of the inertial term focuses on choosing a set of
 orthogonal basis functions, such as the Legendre polynomials with
 Gauss-Lobatto-Lgendre quadrature points
-\cite{Komatitsch:Vilotte:1998}, which naturally produces a diagonal
-matrix for each finite-element cell without introducing any additional
+\cite{Komatitsch:Vilotte:1998}. This discretization (often called the
+spectral element method) naturally produces a diagonal block for each
+finite-element cell without introducing any additional
 approximations. Because the fault slip constraint term also involves
 integration of the products of the basis functions over
-lower-dimension cells, orthogonal basis functions also produces a
-diagonal matrix for this integration.
+lower-dimension cells, orthogonal basis functions also produce a
+diagonal block for this integration.
 
-In contrast, the more traditional finite-element approach does
-introduce additional approximations by constructing a diagonal
-approximation of the sparse matrix. In PyLith we employ one of these
-more traditional approaches, because it produces good approximations
-for many different choices of basis functions and quadrature
-points. For each finite-element cell, we construct a diagonal
-approximation of the integral such that the action on rigid body
-motion is the same for the diagonal approximation of the integral as
-it is for the original integral. The errors associated with this
-approximation are small as long as the deformation occurs at length
-scales significantly larger than the discretization size, which is
-consistent with resolving the wave propagation. Furthermore, in
-contrast to other approaches that choose basis functions or quadrature
-points that affect the accuracy of the elasticity terms as well, this
-approach only affects the accuracy of the terms involved in the
+In contrast, traditional finite-element approaches do introduce
+additional approximations when constructing a diagonal
+approximation. In PyLith we employ one of these traditional
+approaches, because it produces good approximations for many different
+choices of basis functions and quadrature points. For each
+finite-element cell, we construct a diagonal approximation of the
+integral such that the action on rigid body motion is the same for the
+diagonal approximation of the integral as it is for the original
+integral. The errors associated with this approximation are small as
+long as the deformation occurs at length scales significantly larger
+than the discretization size, which is consistent with resolving
+seismic wave propagation accurately. Furthermore, in contrast to other
+approaches that choose basis functions or quadrature points that
+affect the accuracy of all of the finite-element integrations, such as
+choosing quadrature points coincident with the vertices of a cell,
+this approach only affects the accuracy of the terms involved in the
 Jacobian. For consistency in the formulation of the system of
 equations, these approximations are also applied to the inertial term
 and fault slip constraint term when computing the residual.
@@ -663,30 +665,98 @@ and fault slip constraint term when comp
 
 % Schur complement
 
-Second, we leverage the block structure of the off-diagonal blocks
+Second, we leverage the structure of the off-diagonal blocks
 associated with the fault slip constraint in solving the system of
-equations. We compute an initial residual assuming the increment in
-the solution is zero (i.e., $d\vec{u}_m = \vec{0}$ and $d\vec{l}_p =
-\vec{0}$,
+equations via a Schur's complement algorithm. We compute an initial
+residual assuming the increment in the solution is zero (i.e.,
+$d\vec{u}_m = \vec{0}$ and $d\vec{l}_p = \vec{0}$,
 \begin{equation}
-  \vec{r}^* = \begin{matrix} \vec{r}_m^* \\ \vec{r}_p^* \end{matrix} =
-  \begin{matrix} \vec{b}_m \\ \vec{b}_p \end{matrix}
-  - \begin{matrix}
+  \vec{r}^* = \begin{pmatrix} \vec{r}_m^* \\ \vec{r}_p^* \end{pmatrix} =
+  \begin{pmatrix} \vec{b}_m \\ \vec{b}_p \end{pmatrix}
+  - \begin{pmatrix}
     \tensor{K} & \tensor{L}^T \\ \tensor{L} & 0
-  \end{matrix}
-  \begin{matrix} \vec{u}_m \\ \vec{l}_m \end{matrix}.
+  \end{pmatrix}
+  \begin{pmatrix} \vec{u}_m \\ \vec{l}_m \end{pmatrix}.
 \end{equation}
- We compute an initial solution to the system of equations
-$\vec{u}_m^*}$ ignoring the off-diagonal blocks in the Jacobian and
+ We compute a corresponding initial solution to the system of equations
+$\vec{u}_m^*$ ignoring the off-diagonal blocks in the Jacobian and
 the increment in the Lagrange multipliers.
 \begin{equation}
 d\vec{u}_m^* = \tensor{K}^{-1} \vec{r}_m,
 \end{equation}
-taking advantage of the fact that we constructed $\tensor{K}$ so that
-it is diagonal.
+taking advantage of the fact that we construct $\tensor{K}$ so that it
+is diagonal. 
 
+We next compute the increment in the Lagrange multipliers in order to
+correct this initial solution so that the true residual is zero.
+Making use of the initial residual, the expression for the true
+residual is
+\begin{equation}
+  \label{eqn:lumped:jacobian:residual}
+  \vec{r} = \begin{pmatrix} \vec{r}_m \\ \vec{r}_p \end{pmatrix} =
+  \begin{pmatrix} \vec{r}_m^* \\ \vec{r}_p^* \end{pmatrix}
+  - \begin{pmatrix}
+    \tensor{K} & \tensor{L}^T \\ \tensor{L} & 0
+  \end{pmatrix}
+  \begin{pmatrix} d\vec{u}_m \\ d\vec{l}_m \end{pmatrix}.
+\end{equation}
+Solving the first row of equation~(\ref{eqn:lumped:jacobian:residual})
+for the increment in the solution and accounting for the structure of
+$\tensor{L}$ as we write the expressions for degrees of freedom on
+each side of the fault, we have
+\begin{gather}
+  d\vec{u}_{m^+} = 
+    d\vec{u}_{m^+}^* - \tensor{K}_{m^+m^+}^{-1} \tensor{L}_p^T d\vec{l}_p, \\
+  d\vec{u}_{m^-} = 
+    d\vec{u}_{m^-}^* + \tensor{K}_{m^-m^-}^{-1} \tensor{L}_p^T d\vec{l}_p.
+\end{gather}
+Substituting into the second row of
+equation~(\ref{eqn:lumped:jacobian:residual}) and isolating the term
+with the increment in the Lagrange multipliers yields
+\begin{equation}
+  \tensor{L}_p \left( \tensor{K}_{m^+m^+}^{-1} + \tensor{K}_{m^-m^-}^{-1} \right)
+  \tensor{L}^T_p d\vec{l}_p =
+  -\vec{r}_p^* + \tensor{L}_p \left( d\vec{u}_{m^+}^* - d\vec{u}_{m^-}^* \right).
+\end{equation}
+Letting
+\begin{equation}
+  \tensor{S}_p = \tensor{L}_p 
+  \left( \tensor{K}_{m^+m^+}^{-1} + \tensor{K}_{m^-m^-}^{-1} \right)
+  \tensor{L}^T_p,
+\end{equation}
+and recognizing that $S_p$ is diagonal because ${K}$ and ${L}_p$ are
+diagonal allows us to solve for the increment in the Lagrange
+multipliers,
+\begin{equation}
+  d\vec{l}_p = \tensor{S}_p^{-1} \left(
+  -\vec{r}_p^* + \tensor{L}_p \left( d\vec{u}_{m^+}^* - d\vec{u}_{m^-}^* \right)
+  \right).
+\end{equation}
+Now that we have the increment in the Lagrange multipliers, we can
+correct our initial solution $d\vec{u}_m^*$ so that the true residual
+is zero,
+\begin{equation}
+  d\vec{u}_m = d\vec{u}_m^* - \tensor{K}^{-1} \tensor{L}^T d\vec{l}_p.
+\end{equation}
+Because $\tensor{K}$ and $\tensor{L}$ are comprised of diagonal
+blocks, this expression for the updates to the solution are local to
+the degrees of freedom attached to the fault and the Lagrange
+multipliers.
 
-% Spontaneous rupture
+% Spontaneous rupture model and lumped Jacobian
+
+We also leverage the elimination of off-diagonal entries from the
+blocks of the Jacobian in dynamic simulations when updating the slip
+in spontaneous rupture models. Because $\tensor{K}$ is diagonal in
+this case, the expression for the change in slip for a perturbation in
+the Lagrange multipliers
+(equation~(\ref{eqn:spontaneous:rupture:slip:update})) is exact and
+simplifies to
+\begin{equation}
+  \partial \vec{d} = - \tensor{R} 
+  \left( \tensor{K}_{m^+m+}^{-1} + \tensor{K}_{m^-m^-}^{-1} \right)
+  \tensor{L}_p^T \partial \vec{l}_p.
+\end{equation}
 
 
 
diff -r 38890ae93959 -r 32a982452b80 references.bib
--- a/references.bib	Wed Aug 31 16:16:36 2011 -0700
+++ b/references.bib	Thu Sep 01 10:09:41 2011 -0700
@@ -646,7 +646,7 @@
   journal =	 BSSA,
   year =	 1998,
   volume =	 88,
-  number         2,
+  number =       2,
   month =	 apr,
   pages =	 {368--392},
   abstract =     {We present the spectral element method to simulate