[cig-commits] commit: More work on spontaneous fault rupture.
Mercurial
hg at geodynamics.org
Tue Aug 30 15:56:44 PDT 2011
changeset: 55:d7a986bedd87
tag: tip
user: Brad Aagaard <baagaard at usgs.gov>
date: Tue Aug 30 15:56:41 2011 -0700
files: faultRup.tex
description:
More work on spontaneous fault rupture.
diff -r 7faf54660125 -r d7a986bedd87 faultRup.tex
--- a/faultRup.tex Tue Aug 30 12:20:51 2011 -0700
+++ b/faultRup.tex Tue Aug 30 15:56:41 2011 -0700
@@ -71,7 +71,7 @@ postseismic behavior often approximate d
postseismic behavior often approximate dynamic rupture behavior with
the static coseismic slip
\cite{Reilinger:etal:2000,Pollitz:etal:2001,Langbein:etal:2006,Chlieh:etal:2007}.
-Likewise, studies of radid deformation associated with earthquake
+Likewise, studies of rapid deformation associated with earthquake
rupture propagation often approximate the loading of the crust via
simplistic assumptions about the stress field at the beginning of a
rupture
@@ -132,7 +132,7 @@ behavior. Our current focus is developin
behavior. Our current focus is developing a modeling code, PyLith,
that supports both quasi-static simulations of interseismic and
coseismic deformation and dynamic simulations of earthquake rupture
-propapgation. We plan to seamlessly couple these two types of
+propagation. We plan to seamlessly couple these two types of
simulations together to resolve the earthquake cycle.
Other compelling reasons led us to develop such a code with the
@@ -145,13 +145,16 @@ over the domain. The two types of simula
over the domain. The two types of simulations tend to use different
boundary conditions, with interseismic deformation usually driven by
Dirichlet (displacement/velocity) or Neumann (traction) boundary
-conditions and rupture propgation simulations using absorbing
+conditions and rupture propagation simulations using absorbing
boundaries to truncate the domains. However, these features constitute
a small fraction of the code. The primary different between the two
types of simulations is the time integration scheme, with an implicit
scheme used in the quasi-static simulations and an explicit scheme
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)}
Implementing slip on the potentially nonplanar fault surface
differentiates these types of problems from many other elasticity
@@ -288,7 +291,7 @@ leads to
%
\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-}
+ - \tensor{N}_{m^+} \vec{u}_{m^+} + \tensor{N}_{m^-} \vec{u}_{m^-}
\right)
\, dS = \vec{0}.
\end{gather}
@@ -323,7 +326,7 @@ and the loading conditions. Considering
\label{eqn:quasi-static:residual:fault}
\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)
+ - \tensor{N}_{m^+} \vec{u}_{m^+}(\tpdt) + \tensor{N}_{m^-} \vec{u}_{m^-}(\tpdt)
\right)
\, dS = \vec{0}.
\end{gather}
@@ -331,11 +334,11 @@ equations, and add the increment in 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
+tools for solving linear systems of algebraic equations 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)$,
+\left( \begin{smallmatrix} \vec{u}_m \\ \vec{l}_m \end{smallmatrix} \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}). That is,
@@ -352,7 +355,7 @@ and~(\ref{eqn:quasi-static:residual:faul
\\
\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)
+ - \tensor{N}_{m^+} \vec{u}_{m^+}(\tpdt) + \tensor{N}_{m^-} \vec{u}_{m^-}(\tpdt)
\right)
\, dS
\end{array} \right).
@@ -383,7 +386,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}
- \tensor{L} = \int_{S_f} \tensor{N}_p (-\tensor{N}_{m+}+\tensor{N}_{m-}) \, 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}
@@ -396,7 +399,7 @@ Thus, the Jacobian of the entire system
% ------------------------------------------------------------------
\subsection{Dynamic Simulations}
-In dynamic simulations we include the interial term in order to
+In dynamic simulations we include the inertial 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
@@ -418,7 +421,7 @@ inertial term:
\\
\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)
+ - \tensor{N}_{m^+} \vec{u}_{m^+}(\tpdt) + \tensor{N}_{m^-} \vec{u}_{m^-}(\tpdt)
\right)
\, dS
\end{array} \right).
@@ -428,7 +431,7 @@ the action on the increment in the solut
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
+time using explicit time stepping via Newark's 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) =
@@ -456,7 +459,7 @@ so that the upper portion of the Jacobia
% ------------------------------------------------------------------
-% Advantages of common code base
+% Taking advantage of the commonality
The advantages of designing PyLith to handle both quasi-static and
dynamic simulations should now be apparent. There are only minor
@@ -468,38 +471,113 @@ elastic, viscoelastic, and viscoelastopl
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. This significantly reduces
+development time compared to maintaining different codes for
quasi-static and dynamic simulations.
% ------------------------------------------------------------------
\subsection{Prescribed Fault Rupture}
-Slip-time functions
+In a prescribed (kinematic) fault rupture we specify the slip-time
+history $\vec{d}(x,y,z,t)$ at every location on the fault surfaces.
+The slip-time history enters into the calculation of the residual as
+do the Lagrange multipliers, which are available as from the current
+solution. In prescribing the slip-time history we do not specify the
+initial tractions on the fault surface so the Lagrange multipliers are
+the change in the tractions on the fault surfaces corresponding to the
+slip. PyLith includes a variety of slip-time histories, including a
+step function, constant slip rate, the integral of Brune's far-field
+time function, a sine-cosine function developed by
+\citeN{Liu:etal:200?}, and a user-defined time history. These are
+discussed in detail in the PyLith manual
+\cite{PyLith:manual:1.6.2}. PyLith allows multiple earthquake ruptures
+to be superimposed on each other, thereby permitting rather complex
+slip behavior.
+
% ------------------------------------------------------------------
\subsection{Spontaneous Fault Rupture}
-Fault constitutive models
+In a spontaneous (dynamic) fault rupture a constitutive model controls
+the tractions on the fault surface. The fault slip evolves based on
+the boundary conditions and deformation. Our fault implementation
+involves specifying the fault slip and solving for the tractions on
+the fault (Lagrange multipliers). Hence, the fault constitutive model
+places bounds on the Lagrange multipliers and the systems of equations
+is nonlinear- we must adjust the slip to be consistent with the fault
+tractions. When a location on a fault is not slipping, no updates to
+the specified fault slip are necessary and the solution will converge
+in a single iteration of the nonlinear solve (i.e., the nonlinear
+problem simplifies to a linear one). In our quasi-static simulations,
+in general, the fault slips only a small fraction of the simulation
+duration, so the fault implementation yields a linear problem during
+most of the simulation. In dynamic simulations, we employ an
+alternative approach that involves a lumped Jacobian and avoids a
+nonlinear system of equations.
+
+We determine the adjustment to the fault slip for a given perturbation
+in the Lagrange multipliers assuming the rest of the solution is
+unchanged. That is, we consider only the degrees of freedom associated
+with the fault interface when computing how a perturbation in the
+Lagrange multipliers corresponds to a change in fault slip. In terms
+of the general form of a linear system of equations ($\tensor{A}
+\vec{u} = \vec{b}$), our subset of equations has the form
+\begin{equation}
+ \begin{pmatrix}
+ \tensor{K}_{m^+m^+} & 0 & \tensor{L}_p \\
+ 0 & \tensor{K}_{m^-m^-} & -\tensor{L}_p \\
+ \tensor{L}_p & -\tensor{L}_p & 0
+ \end{pmatrix}
+ \begin{pmatrix}
+ \vec{u}_{m^+} \\
+ \vec{u}_{m^-} \\
+ \vec{l}_p \\
+ \end{pmatrix}
+ =
+ \begin{pmatrix}
+ \vec{b}_{m^+} \\
+ \vec{b}_{m^-} \\
+ \vec{b}_p \\
+ \end{pmatrix},
+\end{equation}
+where $m^+$ and $m^-$ 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}_{m^+}$ and
+$\vec{b}_{m^-}$ 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$:
+\begin{gather}
+ \tensor{K}_{m^+m^+} \vec{u}_{m^+} = - \tensor{L}_p \partial \vec{l}_p, \\
+ \tensor{K}_{m^-m^-} \vec{u}_{m^-} = - \tensor{L}_p \partial \vec{l}_p, \\
+ \partial \vec{d} = \tensor{R} \left( \vec{u}_{m^+} - \vec{u}_{m^-} \right).
+\end{gather}
+
+\brad{Mention fault constitutive models currently implemented in PyLith
+(static friction, slip weakening, time weakening, Dieterich-Ruina
+rate- and state friction with aging law).}
% ------------------------------------------------------------------
\section{Finite-Element Mesh Processing}\matt{Rough draft}
- adjusting topology
+\begin{itemize}
+ \item adjusting topology
+ \item cohesive cells
+ \end{itemize}
+
% ------------------------------------------------------------------
-\section{Fault Discretization}\brad{Rough draft}
- cohesive cells
+\section{Solver Customization}
-% ------------------------------------------------------------------
-\section{Solver Customization}\matt{Rough draft}
-
-\subsection{Quasi-static Simulations}
+\subsection{Quasi-static Simulations}\matt{Rough draft}
\begin{itemize}
\item saddle point
\item custom preconditioner
\end{itemize}
-\subsection{Dynamic-static Simulations}
+\subsection{Dynamic-static Simulations}\brad{Rough draft}
\begin{itemize}
\item Lumped Jacobian
More information about the CIG-COMMITS
mailing list