[cig-commits] commit: Started working on solver section for dynamic simulations.
Mercurial
hg at geodynamics.org
Wed Aug 31 15:20:56 PDT 2011
changeset: 59:d61aa0f9edd9
tag: tip
user: Brad Aagaard <baagaard at usgs.gov>
date: Wed Aug 31 15:20:52 2011 -0700
files: faultRup.tex references.bib
description:
Started working on solver section for dynamic simulations.
diff -r 0a787949f1ae -r d61aa0f9edd9 faultRup.tex
--- a/faultRup.tex Wed Aug 31 08:44:33 2011 -0700
+++ b/faultRup.tex Wed Aug 31 15:20:52 2011 -0700
@@ -153,8 +153,8 @@ used in the dynamic simulations; this al
used in the dynamic simulations; this also results in using different
solvers as we will discuss later.
-\brad{Add paragraph on CIG (open-source, community code so needs to be
- robust, flexible, easy to maintain)}
+\brad{TODO: Add paragraph on CIG (open-source, community code so needs
+ to be robust, flexible, easy to maintain)}
Implementing slip on the potentially nonplanar fault surface
differentiates these types of problems from many other elasticity
@@ -168,7 +168,7 @@ and verification of its implementation u
and verification of its implementation using a few benchmarks.
% ------------------------------------------------------------------
-\section{Numerical Model of Fault Slip}\brad{Rough draft}
+\section{Numerical Model of Fault Slip}
In this section we summarize the formulation of the governing
equations using the finite-element method. We augment the conventional
@@ -380,6 +380,7 @@ displacement vector as a linear combinat
displacement vector as a linear combination of basis functions, after
some algebra we find this portion of the Jacobian is
\begin{equation}
+ \label{eqn:jacobian:implicit:stiffness}
\tensor{K} = \frac{1}{4} \int_v
(\nabla^T + \nabla) \tensor{N}_n^T \cdot
\tensorfour{C} \cdot (\nabla + \nabla^T) \tensor{N}_m \, dV.
@@ -387,6 +388,7 @@ Following a similar procedure, we find t
Following a similar procedure, we find the portion of the Jacobian
associated with equation~(\ref{eqn:quasi-static:residual:fault}) is
\begin{equation}
+ \label{eqn:jacobian:constraint}
\tensor{L} = \int_{S_f} \tensor{N}_p^T (-\tensor{N}_{m^+}+\tensor{N}_{m^-}) \, dS.
\end{equation}
Thus, the Jacobian of the entire system has the form,
@@ -454,6 +456,7 @@ Expanding the inertial term yields
\end{equation}
so that the upper portion of the Jacobian is
\begin{equation}
+ \label{eqn:jacobian:explicit:inertia}
\tensor{K} =
\frac{1}{\Delta t^2} \int_{V} \rho \tensor{N}_n^T \tensor{N}_m \, dV.
\end{equation}
@@ -552,7 +555,7 @@ corresponding to a perturbation in the L
\vec{l}_p$:
\begin{gather}
\tensor{K}_{m^+m^+} \partial \vec{u}_{m^+} = - \tensor{L}_p^T \partial \vec{l}_p, \\
- \tensor{K}_{m^-m^-} \partial \vec{u}_{m^-} = - \tensor{L}_p^T \partial \vec{l}_p, \\
+ \tensor{K}_{m^-m^-} \partial \vec{u}_{m^-} = \tensor{L}_p^T \partial \vec{l}_p, \\
\partial \vec{d} = \tensor{R} \left( \partial \vec{u}_{m^+} - \partial \vec{u}_{m^-} \right).
\end{gather}
This estimate is exact if all other degrees of freedom are constrained
@@ -595,12 +598,57 @@ for details).
\item custom preconditioner
\end{itemize}
-\subsection{Dynamic-static Simulations}\brad{Rough draft}
+\subsection{Dynamic Simulations}\brad{Rough draft}
- \begin{itemize}
- \item Lumped Jacobian
- \item Schur complement
- \end{itemize}
+% Lumped Jacobian
+
+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.
+
+The Jacobian for our system of equations involves two terms: the
+inertial term given by equation~(\ref{eqn:jacobian:explicit:inertia})
+and the fault slip constraint term given by
+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, sliminating 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.
+
+The current best available option for eliminating the off-diagonal
+terms 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 eliminates the
+off-diagonal terms in the Jacobian for both the the interial term and
+the fault slip constraint term without introducing any additional
+approximations. 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 reasonable
+approximations for many different choices of basis functions and
+quadrature points. We construct a diagonal approximation such that the
+action on rigid body motion is the same for the diagonal approximation
+as it is for the original Jacobian. 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.
+
+\brad{TODO: Add equations for intertial term and fault slip constraint
+ for lumped Jacobian}
+
+% Schur complement
+
+
+
+% Spontaneous rupture
+
% ------------------------------------------------------------------
diff -r 0a787949f1ae -r d61aa0f9edd9 references.bib
--- a/references.bib Wed Aug 31 08:44:33 2011 -0700
+++ b/references.bib Wed Aug 31 15:20:52 2011 -0700
@@ -638,7 +638,66 @@
pages = {4997--5009},
}
-
+ at article{Komatitsch:Vilotte:1998,
+ author = {Komatitsch, D. and Vilotte, J.~P.},
+ title = {The spectral element method: {An} efficient tool to
+ simulate the seismic response of 2D and 3D
+ geological structures},
+ journal = BSSA,
+ year = 1998,
+ volume = 88,
+ number 2,
+ month = apr,
+ pages = {368--392},
+ abstract = {We present the spectral element method to simulate
+ elastic-wave propagation in realistic geological
+ structures involving complieated free-surface
+ topography and material interfaces for two- and
+ three-dimensional geometries. The spectral element
+ method introduced here is a high-order variational
+ method for the spatial approximation of elastic-wave
+ equations. The mass matrix is diagonal by
+ construction in this method, which drastically
+ reduces the computational cost and allows an
+ efficient parallel implementation. Absorbing
+ boundary conditions are introduced in variational
+ form to simulate unbounded physical domains. The
+ time discretization is based on an energy-momentum
+ conserving scheme that can be put into a classical
+ explicit-implicit predictor/multi-corrector
+ format. Long-term energy conservation and stability
+ properties are illustrated as well as the efficiency
+ of the absorbing conditions. The associated Courant
+ condition behaves as {Delta}tC < O (nel–1/ndN–2),
+ with nel the number of elements, nd the spatial
+ dimension, and N the polynomial order. In practice,
+ a spatial sampling of approximately 5 points per
+ wavelength is found to be very accurate when working
+ with a polynomial degree of N = 8. The accuracy of
+ the method is shown by comparing the spectral
+ element solution to analytical solutions of the
+ classical two-dimensional (2D) problems of Lamb and
+ Garvin. The flexibility of the method is then
+ illustrated by studying more realistic 2D models
+ involving realistic geometries and complex
+ free-boundary conditions. Very accurate modeling of
+ Rayleigh-wave propagation, surface diffraction, and
+ Rayleigh-to-body-wave mode conversion associated
+ with the free-surface curvature are obtained at low
+ computational cost. The method is shown to provide
+ an efficient tool to study the diffraction of
+ elastic waves by three-dimensional (3D) surface
+ topographies and the associated local effects on
+ strong ground motion. Complex amplification
+ patterns, both in space and time, are shown to occur
+ even for a gentle hill topography. Extension to a
+ heterogeneous hill structure is considered. The
+ efficient implementation on parallel distributed
+ memory architectures will allow to perform real-time
+ visualization and interactive physical
+ investigations of 3D amplification phenomena for
+ seismic risk assessment.},
+}
@Book{Zienkiewicz:Taylor:2005,
More information about the CIG-COMMITS
mailing list