[cig-commits] commit: Minor cleanup of quasi-static simulations section. Worked on dynamic simulations section.
Mercurial
hg at geodynamics.org
Tue Aug 30 12:20:57 PDT 2011
changeset: 54:7faf54660125
tag: tip
user: Brad Aagaard <baagaard at usgs.gov>
date: Tue Aug 30 12:20:51 2011 -0700
files: faultRup.tex
description:
Minor cleanup of quasi-static simulations section. Worked on dynamic simulations section.
diff -r 924b986f337e -r 7faf54660125 faultRup.tex
--- a/faultRup.tex Mon Aug 29 17:54:05 2011 -0700
+++ b/faultRup.tex Tue Aug 30 12:20:51 2011 -0700
@@ -11,7 +11,8 @@
\newcommand\matt[1]{{\color{blue}\bf [MATT: #1]}}
\newcommand\charles[1]{{\color{green}\bf [CHARLES: #1]}}
-\newcommand{\tensor}[1]{\underline{#1}}
+\newcommand{\tensor}[1]{\overline{#1}}
+\newcommand{\tensorfour}[1]{\widetilde{#1}}
\newcommand{\tpdt}{t+\Delta t}
\newcommand{\tmdt}{t-\Delta t}
@@ -176,10 +177,10 @@ We solve the elasticity equation includi
We solve the elasticity equation including inertial terms,
\begin{gather}
\rho \frac{\partial^2\vec{u}}{\partial t^2} - \vec{f}
- - \nabla \cdot \tensor{\sigma} = \vec{0} \text{ in }V, \\
+ - \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, \\
- \tensor{R} \left( \vec{u}^{+} - \vec{u}^{-}\right) = \vec{d}
+ \vec{d} - \tensor{R} \left( \vec{u}^{+} - \vec{u}^{-}\right) = \vec{0}
\text{ on }S_f,
\end{gather}
where $\vec{u}$ is the displacement vector, $\rho$ is the mass
@@ -198,20 +199,20 @@ weak form by taking the dot product with
weak form by taking the dot product with a weighting function and
setting the integral over the domain equal to zero,
\begin{equation}
- \int_{V} \left( \nabla \cdot \tensor{\sigma} + \vec{f} -
+ \int_{V} \vec{\phi} \cdot
+ \left( \tensor{\nabla} \cdot \tensor{\sigma} + \vec{f} -
\rho\frac{\partial^{2}\vec{u}}{\partial t^{2}} \right)
- \cdot \vec{\phi}\, dV=0.
+ \, dV=0.
\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
\begin{equation}
-- \int_{V} \tensor{\sigma} : \nabla\vec{\phi} \, dV
-+ \int_{S_T} \vec{T} \cdot \vec{\phi}\, dS
-+ \int_{V} \vec{f} \cdot \vec{\phi} \, dV
-- \int_{V} \rho \frac{\partial^{2}\vec{u}}{\partial t^{2}} \cdot
-\vec{\phi}\, dV
+- \int_{V} \nabla \vec{\phi} : \tensor{\sigma} \, dV
++ \int_{S_T} \vec{\phi} \cdot \vec{T} \, 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}
@@ -230,13 +231,12 @@ the boundary tractions, we add in the co
the boundary tractions, we add in the contributions from integrating
the Lagrange multipliers (fault tractions) over the fault surface,
\begin{equation}
-- \int_{V} \tensor{\sigma} : \nabla\vec{\phi} \, dV
-+ \int_{S_T} \vec{T} \cdot \vec{\phi}\, dS
-+ \int_{S_f^{+}} \vec{l} \cdot \vec{\phi}\, dS
-- \int_{S_f^{-}} \vec{l} \cdot \vec{\phi}\, dS
-+ \int_{V} \vec{f} \cdot \vec{\phi} \, dV
-- \int_{V} \rho \frac{\partial^{2}\vec{u}}{\partial t^{2}} \cdot
-\vec{\phi}\, dV
+- \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_{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
@@ -244,8 +244,8 @@ with the weighting function and setting
with the weighting function and setting the integral over the fault
surface to zero,
\begin{equation}
- \int_{S_f} \left( \vec{u}^{+} - \vec{u}^{-} -
- \tensor{R}_T \vec{d} \right) \cdot \vec{\phi} \, dS = 0.
+ \int_{S_f} \vec{\phi} \cdot
+ \left( \tensor{R}_T \vec{d} - \vec{u}^{+} + \vec{u}^{-} \right) \, dS = 0.
\end{equation}
The rotation matrix $\tensor{R}_T$ is the inverse of
$\tensor{R}$ and is composed of blocks that are the transpose of the
@@ -255,33 +255,45 @@ We express the trial solution $\vec{u}$,
$\vec{\phi}$, and Lagrange multipliers $\vec{l}$ as linear combinations
of basis functions,
\begin{gather}
-\vec{u} = \sum_{m} \vec{u^{m}} N^{m}, \\
-\vec{\phi} = \sum_{n} \vec{a^{n}} N^{n}, \\
-\vec{l} = \sum_{p} \vec{l^{p}} N^{p}.
+\vec{u} = \sum_{m} \vec{u}_m N_m, \\
+\vec{\phi} = \sum_{n} \vec{a}_n N_n, \\
+\vec{l} = \sum_{p} \vec{l}_p N_p.
\end{gather}
Because the weighting function is zero on $S_u$, the number of basis
functions for the trial solution $\vec{u}$ is generally greater than
the number of basis functions for the weighting function $\vec{\phi}$,
i.e., $m > n$. The basis functions for the Lagrange multipliers are
associated with the fault surface, which is a lower dimension than the
-domain, so $p \ll n$ in most cases. The weighting function is
-arbitrary, so the integrands must be zero for all $\vec{a^n}$, which
+domain, so $p \ll n$ in most cases. If we express the linear
+combination of basis functions in terms of a matrix-vector product, we
+have
+\begin{gather}
+\vec{u} = \tensor{N}_m \vec{u}_m, \\
+\vec{\phi} = \tensor{N}_n \vec{a}_n, \\
+\vec{l} = \tensor{N}_p \vec{l}_p.
+\end{gather}
+
+The weighting function is
+arbitrary, so the integrands must be zero for all $\vec{a}_n$, which
leads to
\begin{gather}
-- \int_{V} \tensor{\sigma} : \nabla N^{n} \, dV
-+ \int_{S_T} \vec{T} N^n \, dS
-+ \int_{S_f^{+}} \sum_{p} \vec{l^p} N^p N^n \, dS
-- \int_{S_f^{-}} \sum_{p} \vec{l^p} N^p N^n \, dS
-+ \int_{V} \vec{f} N^n \, dV
-- \int_{V} \rho\sum_{m} \frac{\partial^2 \vec{u^{m}}}{\partial
- t^{2}}N^{m}N^{n}\, dV
+- \int_{V} \nabla \tensor{N}_n^T \cdot \tensor{\sigma} \, dV
++ \int_{S_T} \tensor{N}_n^T \vec{T} \, dS
++ \int_{S_f^{+}} \tensor{N}_n^T \tensor{N}_p \vec{l}_p \, dS
+- \int_{S_f^{-}} \tensor{N}_n^T \tensor{N}_p \vec{l}_p \, dS
++ \int_{V} \tensor{N}_n^T \vec{f} \, dV
+- \int_{V} \rho \tensor{N}_n^T \tensor{N}_m \frac{\partial^2 \vec{u}_m}{\partial
+ t^2} \, dV
=\vec{0}, \\
%
- \int_{S_f} \left( \sum_m \vec{u^{m+}} N^m - \sum_m \vec{u^{m-}} N^m -
- \tensor{R}_T \vec{d} \right) N^p \, dS = \vec{0}.
+ \int_{S_f} \tensor{N}_p^T
+ \left( \tensor{R}_T \vec{d}
+ - \tensor{N}_{m+} \vec{u}_{m+} + \tensor{N}_{m-} \vec{u}_{m-}
+ \right)
+ \, dS = \vec{0}.
\end{gather}
-We want to solve these equations for the coefficients $\vec{u^m}$ and
-$\vec{l^p}$ subject to $\vec{u} = \vec{u}_0 \text{ on }S_u$.
+We want to solve these equations for the coefficients $\vec{u}_m$ and
+$\vec{l}_p$ subject to $\vec{u} = \vec{u}_0 \text{ on }S_u$.
For nonlinear bulk rheologies it is convenient to work with the
increment in stress and strain, so we formulate the solution of the
@@ -289,9 +301,10 @@ equations in terms of the increment in t
$\tpdt$ rather than the solution at time $\tpdt$.
Consequently, rather than constructing a system with the form
$\tensor{A} \vec{u}(\tpdt) = \vec{b}(\tpdt)$, we construct a system
-with the form $\tensor{A} \vec{du} = \vec{b}(\tpdt) - \tensor{A}
-\vec{u}(t)$.
+with the form $\tensor{A} d\vec{u} = \vec{b}(\tpdt) - \tensor{A}
+\vec{u}(t)$, where $\vec{u}(\tpdt) = \vec{u}(t) + d\vec{u}$.
+% ------------------------------------------------------------------
\subsection{Quasi-static Simulations}
For quasi-static simulations we ignore the inertial term and
@@ -300,50 +313,77 @@ and the loading conditions. Considering
$\tpdt$,
\begin{gather}
\label{eqn:quasi-static:residual:domain}
-- \int_{V} \tensor{\sigma}(\tpdt) : \nabla N^{n} \, dV
-+ \int_{S_T} \vec{T}(\tpdt) N^n \, dS
-+ \int_{S_f^{+}} \sum_{p} \vec{l^p}(\tpdt) N^p N^n \, dS
-- \int_{S_f^{-}} \sum_{p} \vec{l^p}(\tpdt) N^p N^n \, dS
-+ \int_{V} \vec{f}(\tpdt) N^n \, dV
+- \int_{V} \nabla \tensor{N}_n^T \cdot \tensor{\sigma}(\tpdt) \, dV
++ \int_{S_T} \tensor{N}_n^T \vec{T}(\tpdt) \, dS
++ \int_{S_f^{+}} \tensor{N}_n^T \tensor{N}_p \vec{l}_p(\tpdt) \, dS
+- \int_{S_f^{-}} \tensor{N}_n^T \tensor{N}_p \vec{l}_p(\tpdt) \, dS
++ \int_{V} \tensor{N}_n^T \vec{f}(\tpdt) \, dV
=\vec{0}, \\
%
\label{eqn:quasi-static:residual:fault}
- \int_{S_f} \left( \sum_m \vec{u^{m+}}(\tpdt) N^m - \sum_m \vec{u^{m-}}(\tpdt) N^m -
- \tensor{R}_T \vec{d}(\tpdt) \right) N^p \, dS = \vec{0}.
+ \int_{S_f} \tensor{N}_p^T
+ \left( \tensor{R}_T \vec{d}(\tpdt)
+ - \tensor{N}_{m+} \vec{u}_{m+}(\tpdt) + \tensor{N}_{m-} \vec{u}_{m-}(\tpdt)
+ \right)
+ \, dS = \vec{0}.
\end{gather}
-We solve these equations using the Portable, Extensible Toolkit for
-Scientific Computation (PETSc), which provides a suite of tools for
-solving linear systems of algebraic equstions with parallel processing
-\cite{PETSc}. In solving the system, we compute the residual (i.e.,
-$\vec{r} = \vec{b} - \tensor{A} \vec{u}$ and the Jacobian of the
-system ($A$). The residual is simply the left hand sides of
+In order to march forward in time, we simply increment time, solve the
+equations, and add the increment in the solution to the previous
+solution. We solve these equations using the Portable, Extensible
+Toolkit for Scientific Computation (PETSc), which provides a suite of
+tools for solving linear systems of algebraic equstions with parallel
+processing \cite{PETSc}. In solving the system, we compute the
+residual (i.e., $\vec{r} = \vec{b} - \tensor{A} \vec{u}$ and the
+Jacobian of the system ($A$). In our case the solution is $\vec{u} =
+\left( \begin{array}{c} \vec{u}_m \\ \vec{l}_m \end{array} \right)$,
+and the residual is simply the left hand sides of
equations~(\ref{eqn:quasi-static:residual:domain})
-and~(\ref{eqn:quasi-static:residual:fault}).
+and~(\ref{eqn:quasi-static:residual:fault}). That is,
+\begin{equation}
+\vec{r} = \left( \begin{array}{c}
+ \begin{aligned}
+ - \int_{V} \nabla \tensor{N}_n^T \cdot \tensor{\sigma}(\tpdt) \, dV
+ + \int_{S_T} \tensor{N}_n^T \vec{T}(\tpdt) \, dS
+ + \int_{S_f^{+}} \tensor{N}_n^T \tensor{N}_p \vec{l}_p(\tpdt) \, dS
+ &- \int_{S_f^{-}} \tensor{N}_n^T \tensor{N}_p \vec{l}_p(\tpdt) \, dS
+ \\
+ &+ \int_{V} \tensor{N}_n^T \vec{f}(\tpdt) \, dV
+ \end{aligned}
+ \\
+ \int_{S_f} \tensor{N}_p^T
+ \left( \tensor{R}_T \vec{d}(\tpdt)
+ - \tensor{N}_{m+} \vec{u}_{m+}(\tpdt) + \tensor{N}_{m-} \vec{u}_{m-}(\tpdt)
+ \right)
+ \, dS
+ \end{array} \right).
+\end{equation}
The Jacobian of the system is the action (operations) on the increment
in the solution. To find the portion of the Jacobian associated with
equation~(\ref{eqn:quasi-static:residual:domain}), we let
$\tensor{\sigma}(\tpdt) = \tensor{\sigma}(t) +
-\tensor{d\sigma}(t)$. The action on the increment of the solution is
-associated with the term $\tensor{d\sigma}(t)$. We approximate the
+d\tensor{\sigma}(t)$. The action on the increment of the solution is
+associated with the term $d\tensor{\sigma}(t)$. We approximate the
increment in the stress tensor using linear elasticity and
infinitesimal strains,
\begin{equation}
- \tensor{d\sigma}(t) = \frac{1}{2} \tensor{C}(t) (\nabla + \nabla^T)
+ d\tensor{\sigma}(t) = \frac{1}{2} \tensorfour{C}(t) \cdot (\nabla + \nabla^T)
\vec{u}(t),
\end{equation}
-where $\tensor{C}$ is the fourth order tensor of elastic
+where $\tensorfour{C}$ is the fourth order tensor of elastic
constants. Substituting into the first term in
-equation~(\ref{eqn:quasi-static:residual:domain}), after some algebra,
-we find this portion of the Jacobian is
+equation~(\ref{eqn:quasi-static:residual:domain}) and expressing the
+displacement vector as a linear combination of basis functions, after
+some algebra we find this portion of the Jacobian is
\begin{equation}
- \tensor{K} = \frac{1}{2} \int_v \tensor{C} (\nabla +
- \nabla^T) : \nabla N^n \, dV.
+ \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.
\end{equation}
-The portion of the Jacobian associated with
-equation~(\ref{eqn:quasi-static:residual:fault}) is
+Following a similar procedure, we find the portion of the Jacobian
+associated with equation~(\ref{eqn:quasi-static:residual:fault}) is
\begin{equation}
- \tensor{L} = \int_{S_f} (N^{m+}-N^{m-})N^p \, dS.
+ \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}
@@ -353,12 +393,90 @@ Thus, the Jacobian of the entire system
\end{array} \right).
\end{equation}
+% ------------------------------------------------------------------
\subsection{Dynamic Simulations}
+In dynamic simulations we include the interial term in order to
+resolve the propagation of seismic waves. The integral equation for
+fault slip constraint remains unchanged, so the corresponding portions
+of the Jacobian and residual are exactly the same as in the
+quasi-static simulations. Thus, the residual is the same as in the
+quasi-static simulations except the upper portion now includes the
+inertial term:
+\begin{equation}
+\vec{r} = \left( \begin{array}{c}
+ \begin{aligned}
+ - \int_{V} \nabla \tensor{N}_n^T \cdot \tensor{\sigma}(\tpdt) \, dV
+ + \int_{S_T} \tensor{N}_n^T \vec{T}(\tpdt) \, dS
+ &+ \int_{S_f^{+}} \tensor{N}_n^T \tensor{N}_p \vec{l}_p(\tpdt) \, dS
+ - \int_{S_f^{-}} \tensor{N}_n^T \tensor{N}_p \vec{l}_p(\tpdt) \, dS
+ \\
+ &+ \int_{V} \tensor{N}_n^T \vec{f}(\tpdt) \, dV
+ - \int_{V} \rho \tensor{N}_n^T \tensor{N}_m
+ \frac{\partial^2 \vec{u}_m}{\partial t^2} \, dV
+ \end{aligned}
+ \\
+ \int_{S_f} \tensor{N}_p^T
+ \left( \tensor{R}_T \vec{d}(\tpdt)
+ - \tensor{N}_{m+} \vec{u}_{m+}(\tpdt) + \tensor{N}_{m-} \vec{u}_{m-}(\tpdt)
+ \right)
+ \, dS
+ \end{array} \right).
+\end{equation}
+
+We find the upper portion of the Jacobian of the system by considering
+the action on the increment in the solution, just as we did for the
+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 Newmarks method with a central
+difference scheme, so that the acceleration and velocity are given by
+\begin{gather}
+ \frac{\partial^2 \vec{u}}{\partial t^2}(t) =
+ \frac{1}{\Delta t^2} \left(
+ d\vec{u} - \vec{u}(t) + \vec{u}(\tmdt)
+ \right) \\
+%
+ \frac{\partial \vec{u}}{\partial t}(t) = \frac{1}{2\Delta t} \left(
+ d\vec{u} + \vec{u}(t) - \vec{u}(\tmdt)
+ \right).
+\end{gather}
+Expanding the inertial term yields
+\begin{equation}
+ - \int_{V} \rho \tensor{N}_n^T \tensor{N}_m \frac{\partial^2 \vec{u}_m}{\partial
+ t^2} \, dV =
+ - \frac{1}{\Delta t^2} \int_{V} \rho \tensor{N}_n^T \tensor{N}_m \left(
+ d\vec{u}_m(t) - \vec{u}_m(t) + \vec{u}_m(\tmdt)
+ \right) \, dV,
+\end{equation}
+so that the upper portion of the Jacobian is
+\begin{equation}
+ \tensor{K} =
+ \frac{1}{\Delta t^2} \int_{V} \rho \tensor{N}_n^T \tensor{N}_m \, dV.
+\end{equation}
+
+
+% ------------------------------------------------------------------
+% Advantages of common code base
+
+The advantages of designing PyLith to handle both quasi-static and
+dynamic simulations should now be apparent. There are only minor
+differences between the formulations for the quasi-static and
+dynamic simulations. As a result both types of simulations can use the
+same data structures and many of the same routines for assembling the
+system of equations. In particular, the integral terms for the bulk
+rheologies and fault behavior are identical, so all of the various
+elastic, viscoelastic, and viscoelastoplastic constitutive models and fault
+rupture behavior via slip-time functions and fault constitutive models
+need only a single implementation in order to be used for both
+quasi-static and dynamic simulations.
+
+
+% ------------------------------------------------------------------
\subsection{Prescribed Fault Rupture}
Slip-time functions
+% ------------------------------------------------------------------
\subsection{Spontaneous Fault Rupture}
Fault constitutive models
@@ -373,10 +491,21 @@ Fault constitutive models
% ------------------------------------------------------------------
\section{Solver Customization}\matt{Rough draft}
+
+\subsection{Quasi-static Simulations}
+
\begin{itemize}
\item saddle point
\item custom preconditioner
\end{itemize}
+
+\subsection{Dynamic-static Simulations}
+
+ \begin{itemize}
+ \item Lumped Jacobian
+ \item Schur complement
+ \end{itemize}
+
% ------------------------------------------------------------------
\section{Code Verification Benchmarks}
More information about the CIG-COMMITS
mailing list