[cig-commits] commit: Fixed some sign errors in the equations. Worked on custom precondtioner section.

Mercurial hg at geodynamics.org
Mon Oct 10 18:48:28 PDT 2011


changeset:   78:5b1a88528e53
tag:         tip
user:        Brad Aagaard <baagaard at usgs.gov>
date:        Mon Oct 10 18:48:21 2011 -0700
files:       faultRup.tex
description:
Fixed some sign errors in the equations. Worked on custom precondtioner section.


diff -r 7191377a3793 -r 5b1a88528e53 faultRup.tex
--- a/faultRup.tex	Mon Oct 03 09:57:29 2011 -0700
+++ b/faultRup.tex	Mon Oct 10 18:48:21 2011 -0700
@@ -197,7 +197,7 @@ We solve the elasticity equation includi
   - \tensor{\nabla} \cdot \tensor{\sigma} = \vec{0} \text{ in }V, \\
   \tensor{\sigma} \cdot \vec{n} = \vec{T} \text{ on }S_T, \\
   \vec{u} = \vec{u}_0 \text{ on }S_u, \\
-  \left( \vec{u}_{+} - \vec{u}_{-}\right) - \vec{d} = \vec{0}
+  \left( \vec{d} - (\vec{u}_{+} - \vec{u}_{-}\right)) = \vec{0}
   \text{ on }S_f, \label{eqn:fault:disp}
 \end{gather}
 where $\vec{u}$ is the displacement vector, $\rho$ is the mass
@@ -211,8 +211,8 @@ associated with both types of boundary c
 associated with both types of boundary conditions simultaneously.
 
 Following a conventional finite-element formulation (ignoring the
-fault surface for a moment), we construct the
-weak form by taking the dot product with a weighting function and
+fault surface for a moment), we construct the weak form by taking the
+dot product of the governing equation with a weighting function and
 setting the integral over the domain equal to zero,
 \begin{equation}
   \int_{V} \vec{\phi} \cdot 
@@ -222,8 +222,7 @@ setting the integral over the domain equ
 \end{equation}
 The weighting function $\vec{\phi}$ is a piecewise differential vector
 field with $\vec{\phi} = \vec{0}$ on $S_u$. After some algebra and
-use of the boundary conditions, we
-have
+use of the boundary conditions, we have
 \begin{equation}
 - \int_{V} \nabla \vec{\phi} : \tensor{\sigma} \, dV
 + \int_{S_T} \vec{\phi} \cdot \vec{T} \, dS
@@ -261,19 +260,21 @@ the Lagrange multipliers (fault traction
 \begin{equation}
 - \int_{V} \nabla\vec{\phi} : \tensor{\sigma} \, dV
 + \int_{S_T} \vec{\phi} \cdot \vec{T} \, dS
-+ \int_{S_f^{+}} \vec{\phi} \cdot \vec{l} \, dS
-- \int_{S_f^{-}} \vec{\phi} \cdot \vec{l} \, dS
+- \int_{S_f^{+}} \vec{\phi} \cdot \vec{l} \, dS
++ \int_{S_f^{-}} \vec{\phi} \cdot \vec{l} \, dS
 + \int_{V} \vec{\phi} \cdot \vec{f} \, dV 
 - \int_{V} \vec{\phi} \cdot \rho \frac{\partial^{2}\vec{u}}{\partial t^{2}} \, dV
 =0.
 \end{equation}
-We also construct the weak form for the constraint associated with
-slip on the fault by taking the dot product of the constraint equation
-with the weighting function and setting the integral over the fault
-surface to zero,
+Our sign convention for the fault normal and fault tractions (tension
+is positive) leads to the Lagrange multiplier terms having the
+opposite sign as the boundary tractions term. We also construct the
+weak form for the constraint associated with slip on the fault by
+taking the dot product of the constraint equation with the weighting
+function and setting the integral over the fault surface to zero,
 \begin{equation}
   \int_{S_f} \vec{\phi} \cdot 
-  \left( \vec{u}_{+} - \vec{u}_{-} - \vec{d}\right) \, dS = 0.
+  \left(\vec{d} - \vec{u}_{+} + \vec{u}_{-} \right) \, dS = 0.
 \end{equation}
 
 We express the weighting function $\vec{\phi}$, trial solution
@@ -300,16 +301,14 @@ product, we have
 \vec{d} = \tensor{N}_p \cdot \vec{d}_p.
 \end{gather}
 
-\brad{Add comment on collocated method somewhere near here.}
-
 The weighting function is arbitrary, so the integrands must be zero
 for all $\vec{a}_m$, which leads to
 \begin{gather}
   \begin{split}
 - \int_{V} \nabla \tensor{N}_m^T \cdot \tensor{\sigma} \, dV
 + \int_{S_T} \tensor{N}_m^T \cdot \vec{T} \, dS
-+ \int_{S_f^{+}} \tensor{N}_m^T \cdot \tensor{N}_p \cdot \vec{l}_p \, dS
-- \int_{S_f^{-}} \tensor{N}_m^T \cdot \tensor{N}_p \cdot \vec{l}_p \, dS
+- \int_{S_f^{+}} \tensor{N}_m^T \cdot \tensor{N}_p \cdot \vec{l}_p \, dS
++ \int_{S_f^{-}} \tensor{N}_m^T \cdot \tensor{N}_p \cdot \vec{l}_p \, dS
 + \int_{V} \tensor{N}_m^T \cdot \vec{f} \, dV \\
 - \int_{V} \rho \tensor{N}_m^T \cdot \tensor{N}_n \cdot \frac{\partial^2 \vec{u}_n}{\partial
   t^2} \, dV
@@ -318,10 +317,10 @@ for all $\vec{a}_m$, which leads to
 \\
 % 
   \int_{S_f} \tensor{N}_p^T \cdot 
-  \left( \tensor{N}_{n^+} \cdot \vec{u}_{n^+} 
-    - \tensor{N}_{n^-} \cdot \vec{u}_{n^-}
-    - \tensor{N}_p \cdot \vec{d}_p \right)
-  \, dS = \vec{0}.
+  \left( \tensor{N}_p \cdot \vec{d}_p
+    - \tensor{N}_{n^+} \cdot \vec{u}_{n^+} 
+    + \tensor{N}_{n^-} \cdot \vec{u}_{n^-}
+    \right) \, dS = \vec{0}.
 \end{gather}
 We want to solve these equations for the coefficients $\vec{u}_n$ and
 $\vec{l}_p$ subject to $\vec{u} = \vec{u}_0 \text{ on }S_u$.
@@ -347,8 +346,8 @@ and the loading conditions. Considering 
   \begin{split}
     - \int_{V} \nabla \tensor{N}_m^T \cdot \tensor{\sigma}(\tpdt) \, dV
     + \int_{S_T} \tensor{N}_m^T \cdot \vec{T}(\tpdt) \, dS
-    + \int_{S_f^{+}} \tensor{N}_m^T \cdot \tensor{N}_p \cdot \vec{l}_p(\tpdt) \, dS
-    - \int_{S_f^{-}} \tensor{N}_m^T \cdot \tensor{N}_p \cdot \vec{l}_p(\tpdt) \, dS
+    - \int_{S_f^{+}} \tensor{N}_m^T \cdot \tensor{N}_p \cdot \vec{l}_p(\tpdt) \, dS
+    + \int_{S_f^{-}} \tensor{N}_m^T \cdot \tensor{N}_p \cdot \vec{l}_p(\tpdt) \, dS
     \\
     + \int_{V} \tensor{N}_m^T \cdot \vec{f}(\tpdt) \, dV
     =\vec{0},
@@ -357,10 +356,10 @@ and the loading conditions. Considering 
 %
 \label{eqn:quasi-static:residual:fault}
   \int_{S_f} \tensor{N}_p^T \cdot
-  \left( \tensor{N}_{n^+} \cdot \vec{u}_{n^+}(\tpdt)
-    - \tensor{N}_{n^-} \cdot \vec{u}_{n^-}(\tpdt)
-    - \tensor{N}_p \cdot \vec{d}_p(\tpdt) \right)
-  \, dS = \vec{0}.
+  \left( \tensor{N}_p \cdot \vec{d}_p(\tpdt)
+    - \tensor{N}_{n^+} \cdot \vec{u}_{n^+}(\tpdt)
+    + \tensor{N}_{n^-} \cdot \vec{u}_{n^-}(\tpdt)
+    \right) \, dS = \vec{0}.
 \end{gather}
 In order to march forward in time, we simply increment time, solve the
 equations, and add the increment in the solution to the previous
@@ -380,17 +379,17 @@ and~(\ref{eqn:quasi-static:residual:faul
     \begin{aligned}
       - \int_{V} \nabla \tensor{N}_m^T \cdot \tensor{\sigma}(\tpdt) \, dV
       + \int_{S_T} \tensor{N}_m^T \cdot \vec{T}(\tpdt) \, dS
-      + \int_{S_f^{+}} \tensor{N}_m^T \cdot \tensor{N}_p \cdot \vec{l}_p(\tpdt) \, dS
-      &- \int_{S_f^{-}} \tensor{N}_m^T \cdot \tensor{N}_p \cdot \vec{l}_p(\tpdt) \, dS
+      - \int_{S_f^{+}} \tensor{N}_m^T \cdot \tensor{N}_p \cdot \vec{l}_p(\tpdt) \, dS
+      &+ \int_{S_f^{-}} \tensor{N}_m^T \cdot \tensor{N}_p \cdot \vec{l}_p(\tpdt) \, dS
       \\
       &+ \int_{V} \tensor{N}_m^T \cdot \vec{f}(\tpdt) \, dV
     \end{aligned}
     \\
     \int_{S_f} \tensor{N}_p^T \cdot
-    \left( \tensor{N}_{n^+} \cdot \vec{u}_{n^+}(\tpdt) -
-      \tensor{N}_{n^-} \cdot \vec{u}_{n^-}(\tpdt) - \vec{d}(\tpdt)
-    \right)
-    \, dS
+    \left( \tensor{N}_p \vec{d}(\tpdt)
+      - \tensor{N}_{n^+} \cdot \vec{u}_{n^+}(\tpdt)
+      + \tensor{N}_{n^-} \cdot \vec{u}_{n^-}(\tpdt) 
+    \right) \, dS
   \end{array} \right).
 \end{equation}
 
@@ -418,10 +417,10 @@ some algebra we find this portion of the
   \tensorfour{C} \cdot (\nabla + \nabla^T) \tensor{N}_n  \, dV.
 \end{equation}
 Following a similar procedure, we find the portion of the Jacobian
-associated with equation~(\ref{eqn:quasi-static:residual:fault}) is
+associated with the constraints, equation~(\ref{eqn:quasi-static:residual:fault}), is
 \begin{equation}
   \label{eqn:jacobian:constraint}
-  \tensor{L} = \int_{S_f} \tensor{N}_p^T \cdot (-\tensor{N}_{n^+}+\tensor{N}_{n^-}) \, dS.
+  \tensor{L} = \int_{S_f} \tensor{N}_p^T \cdot (\tensor{N}_{n^+} - \tensor{N}_{n^-}) \, dS.
 \end{equation}
 Thus, the Jacobian of the entire system has the form,
 \begin{equation}\label{eqn:saddle-point}
@@ -446,8 +445,8 @@ upper portion now includes the inertial 
     \begin{aligned}
       - \int_{V} \nabla \tensor{N}_m^T \cdot \tensor{\sigma}(\tpdt) \, dV
       + \int_{S_T} \tensor{N}_m^T \cdot \vec{T}(\tpdt) \, dS
-      &+ \int_{S_f^{+}} \tensor{N}_m^T \cdot \tensor{N}_p \cdot \vec{l}_p(\tpdt) \, dS
-      - \int_{S_f^{-}} \tensor{N}_m^T \cdot \tensor{N}_p \cdot \vec{l}_p(\tpdt) \, dS
+      &- \int_{S_f^{+}} \tensor{N}_m^T \cdot \tensor{N}_p \cdot \vec{l}_p(\tpdt) \, dS
+      + \int_{S_f^{-}} \tensor{N}_m^T \cdot \tensor{N}_p \cdot \vec{l}_p(\tpdt) \, dS
       \\
       &+ \int_{V} \tensor{N}_m^T \cdot \vec{f}(\tpdt) \, dV
       - \int_{V} \rho \tensor{N}_m^T \cdot \tensor{N}_n \cdot 
@@ -455,11 +454,10 @@ upper portion now includes the inertial 
         \end{aligned}
     \\
     \int_{S_f} \tensor{N}_p^T \cdot
-    \left( \tensor{N}_{n^+} \cdot \vec{u}_{n^+}(\tpdt)
-      - \tensor{N}_{n^-} \cdot \vec{u}_{n^-}(\tpdt)
-      - \tensor{N}_p \cdot \vec{d}_p (\tpdt)
-   \right)
-    \, dS
+    \left( \tensor{N}_p \cdot \vec{d}_p (\tpdt)
+      - \tensor{N}_{n^+} \cdot \vec{u}_{n^+}(\tpdt)
+      + \tensor{N}_{n^-} \cdot \vec{u}_{n^-}(\tpdt)
+   \right) \, dS
   \end{array} \right).
 \end{equation}
 
@@ -468,7 +466,7 @@ quasi-static simulations. In this case w
 quasi-static simulations. In this case we associate the increment in
 the solution with the temporal discretization. We march forward in
 time using explicit time stepping via Newmark's method with a central
-difference scheme, so that the acceleration and velocity are given by
+difference scheme with the acceleration and velocity given by
 \begin{gather}
   \frac{\partial^2 \vec{u}}{\partial t^2}(t) = 
   \frac{1}{\Delta t^2} \left(
@@ -579,11 +577,10 @@ where $n^+$ and $n^-$ refer to the degre
 where $n^+$ and $n^-$ refer to the degrees of freedom associated with
 the positive and negative sides of the fault,
 respectively. Furthermore, we can ignore the terms $\vec{b}_{n^+}$ and
-$\vec{b}_{n^-}$ because they do not change as we change the Lagrange
-multipliers or fault slip. This means we can solve the following
-equations to estimate the change in fault slip $\partial \vec{d}$
-corresponding to a perturbation in the Lagrange multipliers $\partial
-\vec{l}_p$:
+$\vec{b}_{n^-}$ because they remain constant as we change the Lagrange
+multipliers or fault slip. We solve the following equations to
+estimate the change in fault slip $\partial \vec{d}$ corresponding to
+a perturbation in the Lagrange multipliers $\partial \vec{l}_p$:
 \begin{gather}
   \tensor{K}_{n^+n^+} \cdot \partial \vec{u}_{n^+} = 
   - \tensor{L}_p^T \cdot \partial \vec{l}_p, \\
@@ -591,22 +588,22 @@ corresponding to a perturbation in the L
   \tensor{L}_p^T \cdot \partial \vec{l}_p, \\
   \partial \vec{d}_p =  \partial \vec{u}_{n^+} - \partial \vec{u}_{n^-}.
 \end{gather}
-This estimate is exact if all other degrees of freedom are constrained
-(which is generally only the case for simple toy problems), quite good
-if the deformation associated with the fault slip is confined close to
-the fault (which occurs in most simulations because the 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 nonlinear solve requires fewer iterations
-as the estimate of the fault slip becomes more accurate.
+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 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
 of the cohesive stress $T_c$, coefficient of friction, $\mu_f$, and
 normal traction $T_m$,
 \begin{equation}
-  T_f = T_c - \mu_f T_m.
+  T_f = T_c - \mu_f T_n.
 \end{equation}
 $T_f$ in the this equation corresponds to the magnitude of the shear
 traction vector; the shear traction vector is resolved into the
@@ -619,24 +616,25 @@ law (see the PyLith manual \cite{PyLith:
 % ------------------------------------------------------------------
 \section{Finite-Element Mesh Processing}
 
-Like all finite element engines, PyLith must be able to calculate cell
+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
+particular the Sieve package for finite-element support
 \cite{Knepley:Karpeev:2009,SieveFEMPreprint}.
 
-The Sieve API for mesh representation and manipulation is based upon a
-direct acyclic graph (DAG) representation of the \textit{covering}
-relation in a mesh. For example, faces cover cells, edges cover faces,
-etc. By focusing on the key topological relation, 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 almost no dimension or
-shape specific code. However, in order to support faults and fault
-friction models, we must expand this range of operations.
+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
+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
+able to support one, two, and three dimensional meshes, with
+simplicial, hex, and even prismatic cell shapes, while using almost no
+dimension or shape specific code. However, in order to support faults
+and fault friction models, we must expand this range of operations.
 
 In our domain decomposition approach, the finite-element mesh includes
 the fault as an interior surface. This forces alignment of the element
@@ -703,7 +701,7 @@ negative side of the fault. In cases whe
 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
+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 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
@@ -716,7 +714,7 @@ vertices. For each vertex, we examine th
 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
+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
 continue with a breadth-first search of its neighbors until a
 classified cell is found. This search must terminate because there are
@@ -741,22 +739,31 @@ result.
 
 \subsection{Quasi-static Simulations}
 
-In order to solve the large, sparse systems of linear equations arising in quasi-static simulations, we resort to
-preconditioned Krylov subspace methods~\cite{Saad03}. The Krylov iteration generates a small subspace by repeatedly
-applying the linear system operator to a vector. Since sparse matrix-vector multiplication is very scalable, 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~\cite{Saad03,Smith:etal:1996}, the
-method must often be specialized to the particular problem. Below, we describe a preconditioner specialized to the fault
-interface problem.
+In order to solve the large, sparse systems of linear equations
+arising in quasi-static simulations, we employ preconditioned Krylov
+subspace methods~\cite{Saad03}. The Krylov iteration generates a small
+subspace by repeatedly applying the linear system operator to a vector
+\matt{Need this explained in less technical jargon. Most geoscientists
+  are not familiar with ``subspace'' and ``linear system operator''.}
+Since sparse matrix-vector multiplication is very scalable, 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~\cite{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.
 
-The introduction of Lagrange multipliers to implement the fault slip constraints
-produces the saddle-point problem shown in
+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.
+preconditioners, such as the Additive Schwarz Method (ASM)
+\cite{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
+saddle point problems.
+
+\matt{Add a few sentences about algebraic multi-grid}.
 
 The PCFIELDSPLIT \cite{PETSc:manual} preconditioner in PETSc allows
 the user to define sets of unknowns which correspond to different
@@ -800,8 +807,73 @@ example preconditioners and the options 
 \label{fig:solver-options}
 \end{figure}
 
-\matt{Describe custom preconditioner, and show equation for its
-  construction (from manual I think).}
+Another option involves using using the field split preconditioner in
+PETSc in combination with a custom preconditioner matrix for the
+block associated with the Lagrange multipliers. In formulating the
+custom preconditioner, we exploit the structure of the sparse Jacobian
+matrix. Our system Jacobian has the form
+\begin{equation}
+  \tensor{A} = \left( \begin{array}{cc}
+      \tensor{K} & \tensor{L}^T \\
+      \tensor{L} & \tensor{0}
+    \end{array} \right).
+\end{equation}
+We use the Schur complement of block $\tensor{K}$ to examine the form of $\tensor{A}^{-1}$,
+\begin{equation}
+  \tensor{A}^{-1} = \left( \begin{array}{cc}
+    \tensor{K}^{-1} + \tensor{K}^{-1} \tensor{L}^{T} 
+    (-\tensor{L} \tensor{K}^{-1} \tensor{L}^{T})^{-1} \tensor{L} \tensor{K}^{-1} & 
+    -\tensor{K}^{-1} \tensor{L}^{T }(-\tensor{L} \tensor{K}^{-1} \tensor{L}^{T})^{-1} \\
+    -(-\tensor{L} \tensor{K}^{-1} \tensor{L}^{T})^{-1} \tensor{L}
+    \tensor{K}^{-1} &
+    -(\tensor{L} \tensor{K}^{-1} \tensor{L}^T)^{-1}
+  \end{array} \right),
+\end{equation}
+A suitable block diagonal approximation of $\tensor{A}^{-1}$ is
+\begin{equation}
+  \tensor{P}^{-1} = \left( \begin{array}{cc}
+      \tensor{K}^{-1} & 0 \\
+      0 & -(\tensor{L} \tensor{K}^{-1} ]tensor{L}^T)^{-1}
+    \end{array} \right),
+\end{equation}
+which leads to
+\begin{equation}
+  \tensor{P} = \left( \begin{array}{cc}
+    \tensor{P}_\mathit{elasticity} & 0 \\
+    0 & \tensor{P}_\mathit{fault}
+  \end{array} \right)
+  = \left( \begin{array}{cc}
+    \tensor{K} & 0 \\
+    0 & -\tensor{L} \tensor{K}^{-1} \tensor{L}^T
+  \end{array} \right).
+\end{equation}
+
+We focus on evaluating the lower portion of the preconditioning matrix
+associated with the Lagrange multipliers and rely on PETSc
+preconditioners to construct the upper portion. We approximate
+$\tensor{K}^{-1}$ with the inverse of the diagonal portion of
+$\tensor{K}$. $tensor{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, $\tensor{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
+cell. However, if we collocate quadrature points at the cell vertices,
+then only one basis function is nonzero at each quadrature point and
+$\tensor{L}$ becomes block diagonal; this is also true for spectral
+elements with Legendre polynomials and Gauss-Lobatto-Legendre
+quadrature points. This leads to a diagonal matrix for the lower
+portion of the conditioning matrix,
+\begin{equation}
+  \tensor{P}_f = -\tensor{L_p} (\tensor{K}_{n+n+} + \tensor{K}_{n-n-}) \tensor{K}_p^{T}.
+\end{equation}
+
+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 (CITE ML), such that we precondition the degrees of
+freedom for each global coordinate axis independently.
 
 \matt{Table of convergence for a fault example with ASM, block
   tridiagonal, block tridiagonal+custom preconditioner, LSC}



More information about the CIG-COMMITS mailing list