doc/user/slip.tex

\begin{figure} [ht]
\centerline{\includegraphics[width=\figwidth]{figures/Slip1}}
\caption{Domain $\Omega=[0,1]^2$ with a vertical fault of length $0.5$.}
\label{fig:slip.1}
\end{figure}

\section{Slip on a Fault}
\label{Slip CHAP}

In this example we illustrate how to calculate the stress distribution around a fault \index{fault} in the 
Earth's crust caused by a slip \index{slip} through an earthquake.

To simplify the presentation we assume a simple domain $\Omega=[0,1]^2$ with
a vertical fault in its center. We assume that the slip distribution $s\hackscore{i}$ on the fault is known. We want to calculate the distribution of the displacements $u\hackscore{i}$\index{displacement} and stress $\sigma_{ij}$\index{stress} in the domain. We assume an isotropic, linear elastic material model of the form
\begin{eqnarray} \label{Slip  stress}
\sigma\hackscore{ij} & = & \lambda u\hackscore{k,k} \delta\hackscore{ij} + \mu ( u\hackscore{i,j} + u\hackscore{j,i})
\end{eqnarray}
where $\lambda$ and $\mu$ are the Lame coefficients 
\index{Lame coefficients} and $\delta\hackscore{ij}$ denotes the Kronecker symbol\index{Kronecker symbol}.
On the boundary the normal stress is given by
\begin{eqnarray} \label{Slip natural fault}
\sigma\hackscore{ij}n\hackscore{j}=0
\end{eqnarray}
and normal displacements are set to zero:
\begin{eqnarray} \label{Slip constraint}
u\hackscore{i}n\hackscore{i} =0 
\end{eqnarray}
The stress needs to fulfill the momentum equation
\index{momentum equation}
\begin{eqnarray}\label{Slip general problem}
- \sigma\hackscore{ij,j}=0
\end{eqnarray}
This problem is very similar to the elastic deformation problem presented in Section~\ref{ELASTIC CHAP}.

However, we need to address an additional challenge: The displacement $u\hackscore{i}$ is in fact discontinuous across the fault
but we are in the lucky situation that we know the jump of the displacements across the fault. This is
in fact the the given slip $s\hackscore{i}$. So we can split the total distribution $u\hackscore{i}$ into a component $v\hackscore{i}$ 
which is continuous across the fault and the known slip $s\hackscore{i}$
\begin{eqnarray}\label{Slip Split}
u\hackscore{i} = v\hackscore{i} + \frac{1}{2} s^{\pm}\hackscore{i}
\end{eqnarray}
where $s^{\pm}=s$ when right of the fault and $s^{\pm}=-s$ 
when left of the fault. We assume that $s^{\pm}=0$ sufficiently away from the fault.

We insert this into the stress definition~\ref{Slip  stress}
\begin{eqnarray} \label{Slip  stress split}
\sigma\hackscore{ij} & = &
\sigma^c\hackscore{ij} + 
\frac{1}{2} \sigma^s\hackscore{ij}
\end{eqnarray}
 with 
\begin{eqnarray} \label{Slip  stress split 1 }
\sigma^c\hackscore{ij} = \lambda v\hackscore{k,k} \delta\hackscore{ij} + \mu ( v\hackscore{i,j} + v\hackscore{j,i})
\end{eqnarray}
and 
\begin{eqnarray} \label{Slip  stress split 2 }
\sigma^s\hackscore{ij} = \lambda s^{\pm}\hackscore{k,k} \delta\hackscore{ij} + \mu ( s^{\pm}\hackscore{i,j} + s^{\pm}\hackscore{j,i})
\end{eqnarray}
In fact $\sigma^s\hackscore{ij}$ defines a stress jump across the fault. In easy way to construct this 
function is to us a function $\chi$ which is $1$ on the right and $-1$ on the left side from the fault. One can then
set 
\begin{eqnarray} \label{Slip  stress split 23 }
\sigma^s\hackscore{ij} = \chi \cdot  ( \lambda s\hackscore{k,k} \delta\hackscore{ij} + \mu ( s\hackscore{i,j} + s\hackscore{j,i}) )
\end{eqnarray}
assuming that $s$ is extended by zero away from the fault. After inserting~\ref{Slip  stress split} into~\ref{Slip general problem} we get the differential equation
\index{momentum equation}
\begin{eqnarray}\label{Slip general problem 2 }
- \sigma^c\hackscore{ij,j}=\frac{1}{2} \sigma^s\hackscore{ij,j}
\end{eqnarray}
Together with the definition~\ref{Slip  stress split 1 } we have a differential equation for the
continuous function $v_i$. Notice that the boundary condition~\ref{Slip constraint}
and~\ref{Slip natural fault} transfer to $v_i$ and $\sigma^c\hackscore{ij}$ as 
$s$ is zero away from the fault. In Section~\ref{ELASTIC CHAP} we have discussed how this problem
is solved using the \LinearPDE class. We refer to this section for further details. 

To define the fault we use the \class{FaultSystem} class introduced in Section~\ref{Fault System}. 
The following statements define a fault system \var{fs} and add the fault \var{1} 
to the system.
\begin{python}
fs=FaultSystem(dim=2)
fs.addFault(fs.addFault(V0=[0.5,0.25], strikes=90*DEG, ls=0.5, tag=1)
\end{python}
The fault added starts at point $(0.5,0.25)$ has length $0.5$ and is pointing north.
The main purpose of the \class{FaultSystem} class is to define a parameterization
of the fault using a local coordinate system. One can inquire the class
to get the range used to parameterize a fault.
\begin{python}
p0,p1= fs.getW0Range(tag=1)
\end{python}
Typically \var{p0} is equal to zero while \var{p1} is equal to the length of the fault. The parameterization
is given as a mapping from a set of local coordinates onto parameter range (in our case 
the range \var{p0} to \var{p1}). For instance to map the entire domain \var{mydomain} onto the fault 
one can use
\begin{python}
x=mydomain.getX()
p, m=fs.getParametrization(x,tag=1)
\end{python}
Of course there is the problem that not all location are on the fault. For those locations which are on the 
fault \var{m} is set to $1$ otherwise $0$ is used. So on return the values of \var{p} are defining
the value of the fault parameterization (typically the distance from the starting point of the fault along the fault) 
where \var{m} is positive. On all other locations the value if \var{p} is undefined. Now \var{p} can
be used to define a slip distribution on the fault via
\begin{python}
s=m*(p-p0)*(p1-p)/((p1-p0)/2)**2*slip_max*[0.,1.]
\end{python}
Notice that the factor \var{m} which makes sure that \var{s} is zero away from the fault. It is important
that the slip is zero at the ends of the faults.

\begin{figure} [ht]
\centerline{\includegraphics[width=\figwidth]{figures/Slip2}}
\caption{Total Displacement after slip event.}
\label{fig:slip.2}
\end{figure}

We can now put all components together to get the script:
\begin{python}
from esys.escript import *
from esys.escript.linearPDEs import LinearPDE
from esys.escript.models import FaultSystem
from esys.finley import Rectangle
from esys.escript.unitsSI import DEG

#... set some parameters ...
lam=1.
mu=1
slip_max=1.
mydomain = Rectangle(l0=1.,l1=1.,n0=16, n1=16)  # n1 need to be multiple of 4!!!
# .. create the fault system
fs=FaultSystem(dim=2)
fs.addFault(V0=[0.5,0.25], strikes=90*DEG, ls=0.5, tag=1)
# ... create a slip distribution on the fault:
p, m=fs.getParametrization(mydomain.getX(),tag=1)
p0,p1= fs.getW0Range(tag=1)
s=m*(p-p0)*(p1-p)/((p1-p0)/2)**2*slip_max*[0.,1.]
# ... calculate stress according to slip:
D=symmetric(grad(s))
chi, d=fs.getSideAndDistance(D.getFunctionSpace().getX(),tag=1)
sigma_s=(mu*D+lam*trace(D)*kronecker(mydomain))*chi
#... open symmetric PDE ...
mypde=LinearPDE(mydomain)
mypde.setSymmetryOn()
#... set coefficients ...
C=Tensor4(0.,Function(mydomain))
for i in range(mydomain.getDim()):
  for j in range(mydomain.getDim()):
     C[i,i,j,j]+=lam
     C[j,i,j,i]+=mu
     C[j,i,i,j]+=mu
# ... fix displacement in normal direction 
x=mydomain.getX()
msk=whereZero(x[0])*[1.,0.] + whereZero(x[0]-1.)*[1.,0.] \
   +whereZero(x[1])*[0.,1.] + whereZero(x[1]-1.)*[0.,1.]
mypde.setValue(A=C,X=-0.5*sigma_s,q=msk)
#... solve pde ...
mypde.getSolverOptions().setVerbosityOn()
v=mypde.getSolution()
# .. write the displacement to file:
D=symmetric(grad(v))
sigma=(mu*D+lam*trace(D)*kronecker(mydomain))+0.5*sigma_s
saveVTK("slip.vtu",disp=v+0.5*chi*s, stress= sigma)
\end{python}
The script creates the \file{slip.vtu} giving the total displacements and stress. These values are stored as cell 
centered data. See Figure~\ref{fig:slip.2} for the result.
1	gross	2654	\begin{figure} [ht]
2			\centerline{\includegraphics[width=\figwidth]{figures/Slip1}}
3			\caption{Domain $\Omega=[0,1]^2$ with a vertical fault of length $0.5$.}
4			\label{fig:slip.1}
5			\end{figure}
6
7	gross	2647	\section{Slip on a Fault}
8			\label{Slip CHAP}
9
10	gross	2654	In this example we illustrate how to calculate the stress distribution around a fault \index{fault} in the
11			Earth's crust caused by a slip \index{slip} through an earthquake.
12
13			To simplify the presentation we assume a simple domain $\Omega=[0,1]^2$ with
14			a vertical fault in its center. We assume that the slip distribution $s\hackscore{i}$ on the fault is known. We want to calculate the distribution of the displacements $u\hackscore{i}$\index{displacement} and stress $\sigma_{ij}$\index{stress} in the domain. We assume an isotropic, linear elastic material model of the form
15			\begin{eqnarray} \label{Slip stress}
16	gross	2647	\sigma\hackscore{ij} & = & \lambda u\hackscore{k,k} \delta\hackscore{ij} + \mu ( u\hackscore{i,j} + u\hackscore{j,i})
17			\end{eqnarray}
18			where $\lambda$ and $\mu$ are the Lame coefficients
19			\index{Lame coefficients} and $\delta\hackscore{ij}$ denotes the Kronecker symbol\index{Kronecker symbol}.
20			On the boundary the normal stress is given by
21	gross	2654	\begin{eqnarray} \label{Slip natural fault}
22	gross	2647	\sigma\hackscore{ij}n\hackscore{j}=0
23			\end{eqnarray}
24	gross	2654	and normal displacements are set to zero:
25			\begin{eqnarray} \label{Slip constraint}
26			u\hackscore{i}n\hackscore{i} =0
27			\end{eqnarray}
28			The stress needs to fulfill the momentum equation
29			\index{momentum equation}
30			\begin{eqnarray}\label{Slip general problem}
31			- \sigma\hackscore{ij,j}=0
32			\end{eqnarray}
33			This problem is very similar to the elastic deformation problem presented in Section~\ref{ELASTIC CHAP}.
34
35			However, we need to address an additional challenge: The displacement $u\hackscore{i}$ is in fact discontinuous across the fault
36			but we are in the lucky situation that we know the jump of the displacements across the fault. This is
37			in fact the the given slip $s\hackscore{i}$. So we can split the total distribution $u\hackscore{i}$ into a component $v\hackscore{i}$
38			which is continuous across the fault and the known slip $s\hackscore{i}$
39			\begin{eqnarray}\label{Slip Split}
40			u\hackscore{i} = v\hackscore{i} + \frac{1}{2} s^{\pm}\hackscore{i}
41			\end{eqnarray}
42			where $s^{\pm}=s$ when right of the fault and $s^{\pm}=-s$
43			when left of the fault. We assume that $s^{\pm}=0$ sufficiently away from the fault.
44
45			We insert this into the stress definition~\ref{Slip stress}
46			\begin{eqnarray} \label{Slip stress split}
47			\sigma\hackscore{ij} & = &
48			\sigma^c\hackscore{ij} +
49			\frac{1}{2} \sigma^s\hackscore{ij}
50			\end{eqnarray}
51			with
52			\begin{eqnarray} \label{Slip stress split 1 }
53			\sigma^c\hackscore{ij} = \lambda v\hackscore{k,k} \delta\hackscore{ij} + \mu ( v\hackscore{i,j} + v\hackscore{j,i})
54			\end{eqnarray}
55			and
56			\begin{eqnarray} \label{Slip stress split 2 }
57			\sigma^s\hackscore{ij} = \lambda s^{\pm}\hackscore{k,k} \delta\hackscore{ij} + \mu ( s^{\pm}\hackscore{i,j} + s^{\pm}\hackscore{j,i})
58			\end{eqnarray}
59			In fact $\sigma^s\hackscore{ij}$ defines a stress jump across the fault. In easy way to construct this
60			function is to us a function $\chi$ which is $1$ on the right and $-1$ on the left side from the fault. One can then
61			set
62			\begin{eqnarray} \label{Slip stress split 23 }
63			\sigma^s\hackscore{ij} = \chi \cdot ( \lambda s\hackscore{k,k} \delta\hackscore{ij} + \mu ( s\hackscore{i,j} + s\hackscore{j,i}) )
64			\end{eqnarray}
65			assuming that $s$ is extended by zero away from the fault. After inserting~\ref{Slip stress split} into~\ref{Slip general problem} we get the differential equation
66			\index{momentum equation}
67			\begin{eqnarray}\label{Slip general problem 2 }
68			- \sigma^c\hackscore{ij,j}=\frac{1}{2} \sigma^s\hackscore{ij,j}
69			\end{eqnarray}
70	gross	2690	Together with the definition~\ref{Slip stress split 1 } we have a differential equation for the
71	gross	2654	continuous function $v_i$. Notice that the boundary condition~\ref{Slip constraint}
72			and~\ref{Slip natural fault} transfer to $v_i$ and $\sigma^c\hackscore{ij}$ as
73			$s$ is zero away from the fault. In Section~\ref{ELASTIC CHAP} we have discussed how this problem
74			is solved using the \LinearPDE class. We refer to this section for further details.
75
76			To define the fault we use the \class{FaultSystem} class introduced in Section~\ref{Fault System}.
77			The following statements define a fault system \var{fs} and add the fault \var{1}
78			to the system.
79			\begin{python}
80			fs=FaultSystem(dim=2)
81	gross	2663	fs.addFault(fs.addFault(V0=[0.5,0.25], strikes=90*DEG, ls=0.5, tag=1)
82	gross	2654	\end{python}
83	gross	2663	The fault added starts at point $(0.5,0.25)$ has length $0.5$ and is pointing north.
84	gross	2654	The main purpose of the \class{FaultSystem} class is to define a parameterization
85	gross	2690	of the fault using a local coordinate system. One can inquire the class
86	gross	2654	to get the range used to parameterize a fault.
87			\begin{python}
88			p0,p1= fs.getW0Range(tag=1)
89			\end{python}
90			Typically \var{p0} is equal to zero while \var{p1} is equal to the length of the fault. The parameterization
91			is given as a mapping from a set of local coordinates onto parameter range (in our case
92	gross	2690	the range \var{p0} to \var{p1}). For instance to map the entire domain \var{mydomain} onto the fault
93	gross	2654	one can use
94			\begin{python}
95			x=mydomain.getX()
96			p, m=fs.getParametrization(x,tag=1)
97			\end{python}
98			Of course there is the problem that not all location are on the fault. For those locations which are on the
99			fault \var{m} is set to $1$ otherwise $0$ is used. So on return the values of \var{p} are defining
100			the value of the fault parameterization (typically the distance from the starting point of the fault along the fault)
101			where \var{m} is positive. On all other locations the value if \var{p} is undefined. Now \var{p} can
102			be used to define a slip distribution on the fault via
103			\begin{python}
104			s=m(p-p0)(p1-p)/((p1-p0)/2)*2slip_max*[0.,1.]
105			\end{python}
106			Notice that the factor \var{m} which makes sure that \var{s} is zero away from the fault. It is important
107			that the slip is zero at the ends of the faults.
108
109			\begin{figure} [ht]
110			\centerline{\includegraphics[width=\figwidth]{figures/Slip2}}
111			\caption{Total Displacement after slip event.}
112			\label{fig:slip.2}
113			\end{figure}
114
115			We can now put all components together to get the script:
116			\begin{python}
117			from esys.escript import *
118			from esys.escript.linearPDEs import LinearPDE
119			from esys.escript.models import FaultSystem
120			from esys.finley import Rectangle
121	gross	2663	from esys.escript.unitsSI import DEG
122
123	gross	2654	#... set some parameters ...
124			lam=1.
125			mu=1
126			slip_max=1.
127			mydomain = Rectangle(l0=1.,l1=1.,n0=16, n1=16) # n1 need to be multiple of 4!!!
128			# .. create the fault system
129			fs=FaultSystem(dim=2)
130	gross	2663	fs.addFault(V0=[0.5,0.25], strikes=90*DEG, ls=0.5, tag=1)
131	gross	2654	# ... create a slip distribution on the fault:
132			p, m=fs.getParametrization(mydomain.getX(),tag=1)
133			p0,p1= fs.getW0Range(tag=1)
134			s=m(p-p0)(p1-p)/((p1-p0)/2)*2slip_max*[0.,1.]
135			# ... calculate stress according to slip:
136			D=symmetric(grad(s))
137			chi, d=fs.getSideAndDistance(D.getFunctionSpace().getX(),tag=1)
138			sigma_s=(muD+lamtrace(D)kronecker(mydomain))chi
139			#... open symmetric PDE ...
140			mypde=LinearPDE(mydomain)
141			mypde.setSymmetryOn()
142			#... set coefficients ...
143			C=Tensor4(0.,Function(mydomain))
144			for i in range(mydomain.getDim()):
145			for j in range(mydomain.getDim()):
146			C[i,i,j,j]+=lam
147			C[j,i,j,i]+=mu
148			C[j,i,i,j]+=mu
149			# ... fix displacement in normal direction
150			x=mydomain.getX()
151			msk=whereZero(x[0])[1.,0.] + whereZero(x[0]-1.)[1.,0.] \
152			+whereZero(x[1])[0.,1.] + whereZero(x[1]-1.)[0.,1.]
153			mypde.setValue(A=C,X=-0.5*sigma_s,q=msk)
154			#... solve pde ...
155			mypde.getSolverOptions().setVerbosityOn()
156			v=mypde.getSolution()
157			# .. write the displacement to file:
158			D=symmetric(grad(v))
159			sigma=(muD+lamtrace(D)kronecker(mydomain))+0.5sigma_s
160			saveVTK("slip.vtu",disp=v+0.5chis, stress= sigma)
161			\end{python}
162			The script creates the \file{slip.vtu} giving the total displacements and stress. These values are stored as cell
163			centered data. See Figure~\ref{fig:slip.2} for the result.