1 

2 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
3 
% 
4 
% Copyright (c) 20032010 by University of Queensland 
5 
% Earth Systems Science Computational Center (ESSCC) 
6 
% http://www.uq.edu.au/esscc 
7 
% 
8 
% Primary Business: Queensland, Australia 
9 
% Licensed under the Open Software License version 3.0 
10 
% http://www.opensource.org/licenses/osl3.0.php 
11 
% 
12 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
13 

14 
\section{Seismic Wave Propagation in Two Dimensions} 
15 
\editor{This chapter aims to expand the readers understanding of escript by 
16 
modelling the wave equations. 
17 
Challenges will include a second order differential (multiple initial 
18 
conditions). A new PDE to fit to the general form. Movement to a 3D problem 
19 
(maybe)??? } 
20 

21 
\sslist{example08a.py} 
22 

23 
We will now expand upon the previous chapter by introducing a vector form of 
24 
the wave equation. This means that the waves will have both a scalar magnitude, 
25 
but also a direction. This type of scenario is apparent in wave forms that 
26 
exhibit compressional and transverse particle motion. A common type of wave 
27 
that obeys this principle are seismic waves. 
28 

29 
Wave propagation in the earth can be described by the elastic wave equation: 
30 
\begin{equation} \label{eqn:wav} \index{wave equation} 
31 
\rho \frac{\partial^{2}u\hackscore{i}}{\partial t^2}  \frac{\partial 
32 
\sigma\hackscore{ij}}{\partial x\hackscore{j}} = 0 
33 
\end{equation} 
34 
where $\sigma$ is the stress given by: 
35 
\begin{equation} \label{eqn:sigw} 
36 
\sigma \hackscore{ij} = \lambda u\hackscore{k,k} \delta\hackscore{ij} + \mu ( 
37 
u\hackscore{i,j} + u\hackscore{j,i}) 
38 
\end{equation} 
39 
where $\lambda$ and $\mu$ are the Lame Coefficients. Specifically $\mu$ is the 
40 
bulk modulus. The \refEq{eqn:wav} can be written with the Einstein summation 
41 
convention as: 
42 
\begin{equation} 
43 
\rho u\hackscore{i,tt} = \sigma\hackscore{ij,j} 
44 
\end{equation} 
45 

46 
In a similar process to the previous chapter, we will use the acceleration 
47 
solution to solve this PDE. By substituting $a$ directly for 
48 
$\frac{\partial^{2}u\hackscore{i}}{\partial t^2}$ we can derive the 
49 
displacement solution. Using $a$ \refEq{eqn:wav} becomes; 
50 
\begin{equation} \label{eqn:wava} 
51 
\rho a\hackscore{i}  \frac{\partial 
52 
\sigma\hackscore{ij}}{\partial x\hackscore{j}} = 0 
53 
\end{equation} 
54 

55 
\section{Vector source on the boundary} 
56 
For this particular example, we will introduce the source by applying a 
57 
displacment to the boundary during the initial time steps. The source will again 
58 
be 
59 
a radially propagating wave but due to the vector nature of the PDE used, a 
60 
direction will need to be applied to the source. 
61 

62 
The first step is to choose an amplitude and create the source as in the 
63 
previous chapter. 
64 
\begin{python} 
65 
src_length = 20; print "src_length = ",src_length 
66 
# set initial values for first two time steps with source terms 
67 
y=U0*(cos(length(xxc)*3.1415/src_length)+1)*whereNegative(length(xxc)src_leng 
68 
th) 
69 
\end{python} 
70 
where \verb xc is the source point on the boundary of the model. The source 
71 
direction is then defined as an $(x,y)$ array and multiplied by the source 
72 
function. The directional array must have a magnitude of $1$ otherwise the 
73 
amplitude of the source will become modified. For this example, the source is 
74 
directed in the $y$ direction. 
75 
\begin{python} 
76 
src_dir=numpy.array([0.,1.]) # defines direction of point source as down 
77 
y=y*src_dir 
78 
\end{python} 
79 
The function can then be applied as a boundary condition by setting it equal to 
80 
$y$ in the general form. 
81 
\begin{python} 
82 
mypde.setValue(y=y) #set the source as a function on the boundary 
83 
\end{python} 
84 
Because we are no longer using the source to define our initial condition to 
85 
the model, we must set the model state to zero for the first two time steps. 
86 
\begin{python} 
87 
# initial value of displacement at point source is constant (U0=0.01) 
88 
# for first two time steps 
89 
u=[0.0,0.0]*whereNegative(x) 
90 
u_m1=u 
91 
\end{python} 
92 

93 
If the source will introduce energy to the system over a period longer than one 
94 
or two time steps (ie the initial conditions), $y$ can be updated during the 
95 
iteration stage. 
96 

97 
\section{Time variant source} 
98 

99 
\sslist{example08b.py} 
100 

101 
Until this point, all of the wave propagation examples in this cookbook have 
102 
used impulsive sources which are smooth in space but not time. It is however, 
103 
advantageous to have a time smoothed source as it can reduce the temporal 
104 
frequency range and thus mitigate aliasing in the solution. 
105 

106 

107 
It is quite 
108 
simple to implement a source which is smooth in time. In addition to the 
109 
original source function the only extra requirement is a time function. For 
110 
this example it will be a decaying sinusoidal curve defined by; 
111 
\begin{python} 
112 
U0=0.1 # amplitude of point source 
113 
ls=100 # length of the source 
114 
source=np.zeros(ls,'float') # source array 
115 
g=np.log(0.01)/ls 
116 
for t in range(0,ls): 
117 
source[t]=np.exp(g*t)*U0*np.sin(2.*np.pi*t/(0.75*ls)) 
118 
\end{python} 
119 

120 
We then build the source and the first two time steps via; 
121 
\begin{python} 
122 
# set initial values for first two time steps with source terms 
123 
y=source[0] 
124 
*(cos(length(xxc)*3.1415/src_length)+1)*whereNegative(length(xxc)src_length) 
125 
src_dir=numpy.array([0.,1.]) # defines direction of point source as down 
126 
y=y*src_dir 
127 
mypde.setValue(y=y) #set the source as a function on the boundary 
128 
# initial value of displacement at point source is constant (U0=0.01) 
129 
# for first two time steps 
130 
u=[0.0,0.0]*whereNegative(x) 
131 
u_m1=u 
132 
\end{python} 
133 

134 
Finally, for the length of the source, we are required to update each new 
135 
solution in the itterative section of the solver. This is done via; 
136 
\begin{python} 
137 
# increment loop values 
138 
t=t+h; n=n+1 
139 
if (n < ls): 
140 
y=source[n]**(cos(length(xxc)*3.1415/src_length)+1)*\ 
141 
whereNegative(length(xxc)src_length) 
142 
y=y*src_dir; mypde.setValue(y=y) #set the source as a function on the 
143 
boundary 
144 
\end{python} 
145 

146 
\section{Absorbing Boundary Conditions} 
147 
To mitigate the effect of the boundary on the model, absorbing boundary 
148 
conditions can be introduced. These conditions effectively dampen the wave 
149 
energy as they approach the bounday and thus prevent that energy from being 
150 
reflected. This type of approach is used typically when a model only represents 
151 
a small portion of the entire model, which in reality may have infinite bounds. 
152 
It is inpractical to calculate the solution for an infinite model and thus ABCs 
153 
allow us the create an approximate solution with small to zero boundary effects 
154 
on a model with a solvable size. 
155 

156 
To dampen the waves, the method of Cerjan(1985) 
157 
\footnote{\textit{A nonreflecting boundary condition for discrete acoustic and 
158 
elastic wave equations}, 1985, Cerjan C, Geophysics 50, doi:10.1190/1.1441945} 
159 
where the solution and the stress are multiplied by a damping function defined 
160 
on $n$ nodes of the domain adjacent to the boundary, given by; 
161 
\begin{equation} 
162 
y= 
163 
\end{equation} 
164 
This is applied to the bounding 2050 pts of the model using the location 
165 
specifiers of \esc; 
166 
\begin{python} 
167 
# Define where the boundary decay will be applied. 
168 
bn=30. 
169 
bleft=xstep*bn; bright=mx(xstep*bn); bbot=my(ystep*bn) 
170 
# btop=ystep*bn # don't apply to force boundary!!! 
171 

172 
# locate these points in the domain 
173 
left=x[0]bleft; right=x[0]bright; bottom=x[1]bbot 
174 

175 
tgamma=0.98 # decay value for exponential function 
176 
def calc_gamma(G,npts): 
177 
func=np.sqrt(abs(1.*np.log(G)/(npts**2.))) 
178 
return func 
179 

180 
gleft = calc_gamma(tgamma,bleft) 
181 
gright = calc_gamma(tgamma,bleft) 
182 
gbottom= calc_gamma(tgamma,ystep*bn) 
183 

184 
print 'gamma', gleft,gright,gbottom 
185 

186 
# calculate decay functions 
187 
def abc_bfunc(gamma,loc,x,G): 
188 
func=exp(1.*(gamma*abs(locx))**2.) 
189 
return func 
190 

191 
fleft=abc_bfunc(gleft,bleft,x[0],tgamma) 
192 
fright=abc_bfunc(gright,bright,x[0],tgamma) 
193 
fbottom=abc_bfunc(gbottom,bbot,x[1],tgamma) 
194 
# apply these functions only where relevant 
195 
abcleft=fleft*whereNegative(left) 
196 
abcright=fright*wherePositive(right) 
197 
abcbottom=fbottom*wherePositive(bottom) 
198 
# make sure the inside of the abc is value 1 
199 
abcleft=abcleft+whereZero(abcleft) 
200 
abcright=abcright+whereZero(abcright) 
201 
abcbottom=abcbottom+whereZero(abcbottom) 
202 
# multiply the conditions together to get a smooth result 
203 
abc=abcleft*abcright*abcbottom 
204 
\end{python} 
205 
Note that the boundary conditions are not applied to the surface, as this is 
206 
effectively a free surface where normal reflections would be experienced. The 
207 
resulting boundary damping function can be viewed in \ref{fig:abconds} 
208 

209 
%\begin{figure}{ht} 
210 
%\centering 
211 
%\includegraphics[width=5in]{figures/abconds.png} 
212 
%\end{figure} 
213 

214 
