doc/user/slip.tex


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Copyright (c) 2003-2018 by The University of Queensland
% http://www.uq.edu.au
%
% Primary Business: Queensland, Australia
% Licensed under the Apache License, version 2.0
% http://www.apache.org/licenses/LICENSE-2.0
%
% Development until 2012 by Earth Systems Science Computational Center (ESSCC)
% Development 2012-2013 by School of Earth Sciences
% Development from 2014 by Centre for Geoscience Computing (GeoComp)
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

\section{Slip on a Fault}\label{Slip CHAP}
\begin{figure}[ht]
\centerline{\includegraphics{Slip1}}
\caption{Domain $\Omega=[0,1]^2$ with a vertical fault of length $0.5$.}
\label{fig:slip.1}
\end{figure}
%
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 as illustrated in \fig{fig:slip.1}.
We assume that the slip distribution $s_{i}$ on the fault is known.
We want to calculate the distribution of the displacements $u_{i}$\index{displacement}
and stress $\sigma_{ij}$\index{stress} in the domain.
Further, we assume an isotropic, linear elastic material model of the form
\begin{eqnarray} \label{Slip  stress}
\sigma_{ij} & = & \lambda u_{k,k} \delta_{ij} + \mu ( u_{i,j} + u_{j,i})
\end{eqnarray}
where $\lambda$ and $\mu$ are the Lam\'e coefficients\index{Lam\'e coefficients}
and $\delta_{ij}$ denotes the Kronecker symbol\index{Kronecker symbol}.
On the boundary the normal stress is given by
\begin{eqnarray} \label{Slip natural fault}
\sigma_{ij}n_{j}=0
\end{eqnarray}
and normal displacements are set to zero:
\begin{eqnarray} \label{Slip constraint}
u_{i}n_{i} =0
\end{eqnarray}
The stress needs to fulfill the momentum equation\index{momentum equation}
\begin{eqnarray}\label{Slip general problem}
- \sigma_{ij,j}=0
\end{eqnarray}
This problem is very similar to the elastic deformation problem presented in \Sec{ELASTIC CHAP}.
However, we need to address an additional challenge: the displacement
$u_{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 given slip $s_{i}$.
So we can split the total distribution $u_{i}$ into a component
$v_{i}$ which is continuous across the fault and the known slip $s_{i}$
\begin{eqnarray}\label{Slip Split}
u_{i} = v_{i} + \frac{1}{2} s^{\pm}_{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$ when sufficiently away from the fault.

We insert this into the stress definition in \eqn{Slip stress}
\begin{eqnarray} \label{Slip stress split}
\sigma_{ij} & = &
\sigma^c_{ij} +
\frac{1}{2} \sigma^s_{ij}
\end{eqnarray}
 with
\begin{eqnarray} \label{Slip stress split 1}
\sigma^c_{ij} = \lambda v_{k,k} \delta_{ij} + \mu ( v_{i,j} + v_{j,i})
\end{eqnarray}
and
\begin{eqnarray} \label{Slip stress split 2}
\sigma^s_{ij} = \lambda s^{\pm}_{k,k} \delta_{ij} + \mu ( s^{\pm}_{i,j} + s^{\pm}_{j,i}).
\end{eqnarray}
In fact, $\sigma^s_{ij}$ defines a stress jump across the fault.
An easy way to construct this function is to use 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_{ij} = \chi \cdot  ( \lambda s_{k,k} \delta_{ij} + \mu ( s_{i,j} + s_{j,i}) )
\end{eqnarray}
assuming that $s$ is extended by zero away from the fault.
After inserting \eqn{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_{ij,j}=\frac{1}{2} \sigma^s_{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_{ij}$ as $s$ is zero away from the fault.
In \Sec{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 \Sec{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 points 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 a 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 locations 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} define 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 of \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 the factor \var{m} which ensures that \var{s} is zero away from the fault.
It is important that the slip is zero at the ends of the faults.

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.weipa import saveVTK
  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 needs to be a 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 \file{slip.vtu} which contains the total
displacements and stress.
These values are stored as cell-centered data.
%
\begin{figure} [ht]
\centerline{\includegraphics[width=\figwidth]{Slip2}}
\caption{Total Displacement after the slip event}
\label{fig:slip.2}
\end{figure}
%
See \fig{fig:slip.2} for a visualization of the result.

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