[cig-commits] commit: More writing on solver section
Mercurial
hg at geodynamics.org
Thu Sep 1 18:33:53 PDT 2011
changeset: 71:4dde491c41e8
user: Matthew G. Knepley <knepley at gmail.com>
date: Thu Sep 01 13:17:57 2011 -0500
files: .hgignore faultRup.tex references.bib
description:
More writing on solver section
diff -r ec74c801eb0d -r 4dde491c41e8 .hgignore
--- a/.hgignore Thu Sep 01 06:26:09 2011 -0500
+++ b/.hgignore Thu Sep 01 13:17:57 2011 -0500
@@ -7,3 +7,4 @@ syntax: glob
*.log
*.pdf
faultRup\.out
+faultRup\.fff
diff -r ec74c801eb0d -r 4dde491c41e8 faultRup.tex
--- a/faultRup.tex Thu Sep 01 06:26:09 2011 -0500
+++ b/faultRup.tex Thu Sep 01 13:17:57 2011 -0500
@@ -16,6 +16,7 @@
\newcommand{\tensorfour}[1]{\widetilde{#1}}
\newcommand{\tpdt}{t+\Delta t}
\newcommand{\tmdt}{t-\Delta t}
+\newcommand{\code}[1]{{\tt #1}}
\widowpenalty=999999
\clubpenalty=999999
@@ -185,7 +186,7 @@ We solve the elasticity equation includi
\tensor{\sigma} \cdot \vec{n} = \vec{T} \text{ on }S_T, \\
\vec{u} = \vec{u}_0 \text{ on }S_u, \\
\vec{d} - \tensor{R} \left( \vec{u}_{+} - \vec{u}_{-}\right) = \vec{0}
- \text{ on }S_f,
+ \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
@@ -390,7 +391,7 @@ associated with equation~(\ref{eqn:quasi
\tensor{L} = \int_{S_f} \tensor{N}_p (-\tensor{N}_{m^+}+\tensor{N}_{m^-}) \, dS.
\end{equation}
Thus, the Jacobian of the entire system has the form,
-\begin{equation}
+\begin{equation}\label{eqn:saddle-point}
\tensor{A} =
\left( \begin{array}{cc}
\tensor{K} & \tensor{L}^T \\ \tensor{L} & \tensor{0}
@@ -577,22 +578,126 @@ for details).
for details).
% ------------------------------------------------------------------
-\section{Finite-Element Mesh Processing}\matt{Rough draft}
+\section{Finite-Element Mesh Processing}
-\begin{itemize}
- \item adjusting topology
- \item cohesive cells
- \end{itemize}
+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.
+
+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.
+
+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.
+
+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}.
+
+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.
+
+\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 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?}
% ------------------------------------------------------------------
\section{Solver Customization}
-\subsection{Quasi-static Simulations}\matt{Rough draft}
+\subsection{Quasi-static Simulations}
- \begin{itemize}
- \item saddle point
- \item custom preconditioner
- \end{itemize}
+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.
+
+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.
+
+\begin{figure}
+\centering
+\begin{tabular}{ll}
+ $\begin{pmatrix}\hat K & 0 \\ 0 & I\end{pmatrix}$ & $\begin{pmatrix}\hat K & L^T \\ 0 & I\end{pmatrix}$ \\
+ \code{-pc\_type fieldsplit} & \code{-pc\_type fieldsplit} \\
+ \code{-pc\_field\_split\_type additive} & \code{-pc\_field\_split\_type multiplicative} \\
+ \code{-fieldsplit\_0\_pc\_type ml} & \code{-fieldsplit\_0\_pc\_type ml} \\
+ \code{-fieldsplit\_0\_ksp\_type preonly} & \code{-fieldsplit\_0\_ksp\_type preonly} \\
+ \code{-fieldsplit\_1\_pc\_type jacobi} & \code{-fieldsplit\_1\_pc\_type jacobi} \\
+ \code{-fieldsplit\_1\_ksp\_type preonly} & \code{-fieldsplit\_1\_ksp\_type preonly} \\
+ \smallskip \\
+ $\begin{pmatrix}\hat K & 0 \\ 0 & -\hat S\end{pmatrix}$ & $\begin{pmatrix}\hat K & 0 \\ L & \hat S\end{pmatrix}$ \\
+ \code{-pc\_type fieldsplit} & \code{-pc\_type fieldsplit} \\
+ \code{-pc\_field\_split\_type schur} & \code{-pc\_field\_split\_type schur} \\
+ \code{-fieldsplit\_0\_pc\_type ml} & \code{-fieldsplit\_0\_pc\_type ml} \\
+ \code{-fieldsplit\_0\_ksp\_type preonly} & \code{-fieldsplit\_0\_ksp\_type preonly} \\
+ \code{-fieldsplit\_1\_pc\_type none} & \code{-fieldsplit\_1\_pc\_type none} \\
+ \code{-fieldsplit\_1\_ksp\_type minres} & \code{-fieldsplit\_1\_ksp\_type minres} \\
+ \code{-pc\_fieldsplit\_schur\_factorization\_type diag} & \code{-pc\_fieldsplit\_schur\_factorization\_type lower} \\
+ \smallskip \\
+ $\begin{pmatrix}\hat K & L^T \\ 0 & \hat S\end{pmatrix}$ & $\begin{pmatrix}I & 0 \\ B^T A^{-1} & I\end{pmatrix}\begin{pmatrix}\hat A & 0 \\ 0 & \hat S\end{pmatrix}\begin{pmatrix}I & A^{-1} B \\ 0 & I\end{pmatrix}$ \\
+ \code{-pc\_type fieldsplit} & \code{-pc\_type fieldsplit} \\
+ \code{-pc\_field\_split\_type schur} & \code{-pc\_field\_split\_type schur} \\
+ \code{-fieldsplit\_0\_pc\_type ml} & \code{-fieldsplit\_0\_pc\_type ml} \\
+ \code{-fieldsplit\_0\_ksp\_type preonly} & \code{-fieldsplit\_0\_ksp\_type preonly} \\
+ \code{-fieldsplit\_1\_pc\_type none} & \code{-fieldsplit\_1\_pc\_type none} \\
+ \code{-fieldsplit\_1\_ksp\_type minres} & \code{-fieldsplit\_1\_ksp\_type minres} \\
+ \code{-pc\_fieldsplit\_schur\_factorization\_type upper} & \code{-pc\_fieldsplit\_schur\_factorization\_type full}
+\end{tabular}
+\caption{Example preconditioners for the saddle point problem in equation~(\ref{eqn:saddle-point}).}
+\label{fig:solver-options}
+\end{figure}
+
+\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}
\subsection{Dynamic-static Simulations}\brad{Rough draft}
@@ -618,6 +723,8 @@ for details).
\begin{acknowledgments}
\end{acknowledgments}
+MGK acknowledges partial support from NSF grant EAR-0949446.
+
% ------------------------------------------------------------------
\bibliography{references}
diff -r ec74c801eb0d -r 4dde491c41e8 references.bib
--- a/references.bib Thu Sep 01 06:26:09 2011 -0500
+++ b/references.bib Thu Sep 01 13:17:57 2011 -0500
@@ -638,3 +638,14 @@
howpublished = {\url{http://www.mcs.anl.gov/petsc}},
year = {2011}
}
+
+ at article{KnepleyKarpeev09,
+ author = {Matthew G. Knepley and Dmitry A. Karpeev},
+ title = {Mesh Algorithms for {PDE} with {Sieve} {I}: {Mesh} Distribution},
+ journal = {Scientific Programming},
+ volume = {17},
+ number = {3},
+ pages = {215--230},
+ year = {2009},
+ doi = {10.3233/SPR-2009-0249}
+}
More information about the CIG-COMMITS
mailing list