1 

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

15 
\section{Seismic Wave Propagation in Two Dimensions} 
16 

17 
\sslist{example08a.py} 
18 
We will now expand upon the previous chapter by introducing a vector form of 
19 
the wave equation. This means that the waves will have not only a scalar 
20 
magnitude as for the pressure wave solution, but also a direction. This type of 
21 
scenario is apparent in wave types that exhibit compressional and transverse 
22 
particle motion. An example of this would be seismic waves. 
23 

24 
Wave propagation in the earth can be described by the elastic wave equation 
25 
\begin{equation} \label{eqn:wav} \index{wave equation} 
26 
\rho \frac{\partial^{2}u_{i}}{\partial t^2}  \frac{\partial 
27 
\sigma_{ij}}{\partial x_{j}} = 0 
28 
\end{equation} 
29 
where $\sigma$ is the stress given by 
30 
\begin{equation} \label{eqn:sigw} 
31 
\sigma _{ij} = \lambda u_{k,k} \delta_{ij} + \mu ( 
32 
u_{i,j} + u_{j,i}) 
33 
\end{equation} 
34 
and $\lambda$ and $\mu$ represent Lame's parameters. Specifically for seismic 
35 
waves, $\mu$ is the propagation materials shear modulus. 
36 
In a similar process to the previous chapter, we will use the acceleration 
37 
solution to solve this PDE. By substituting $a$ directly for 
38 
$\frac{\partial^{2}u_{i}}{\partial t^2}$ we can derive the 
39 
acceleration solution. Using $a$ we can see that \autoref{eqn:wav} becomes 
40 
\begin{equation} \label{eqn:wava} 
41 
\rho a_{i}  \frac{\partial 
42 
\sigma_{ij}}{\partial x_{j}} = 0 
43 
\end{equation} 
44 
Thus the problem will be solved for acceleration and then converted to 
45 
displacement using the backwards difference approximation as for the acoustic 
46 
example in the previous chapter. 
47 

48 
Consider now the stress $\sigma$. One can see that the stress consists of two 
49 
distinct terms: 
50 
\begin{subequations} 
51 
\begin{equation} \label{eqn:sigtrace} 
52 
\lambda u_{k,k} \delta_{ij} 
53 
\end{equation} 
54 
\begin{equation} \label{eqn:sigtrans} 
55 
\mu (u_{i,j} + u_{j,i}) 
56 
\end{equation} 
57 
\end{subequations} 
58 
One simply recognizes in \autoref{eqn:sigtrace} that $u_{k,k}$ is the 
59 
trace of the displacement solution and that $\delta_{ij}$ is the 
60 
kronecker delta function with dimensions equivalent to $u$. The second term 
61 
\autoref{eqn:sigtrans} is the sum of $u$ with its own transpose. Putting these 
62 
facts together we see that the spatial differential of the stress is given by the 
63 
gradient of $u$ and the aforementioned operations. This value is then submitted 
64 
to the \esc PDE as $X$. 
65 
\begin{python} 
66 
g=grad(u); stress=lam*trace(g)*kmat+mu*(g+transpose(g)) 
67 
mypde.setValue(X=stress) # set PDE values 
68 
\end{python} 
69 
The solution is then obtained via the usual method and the displacement is 
70 
calculated so that the memory variables can be updated for the next time 
71 
iteration. 
72 
\begin{python} 
73 
accel = mypde.getSolution() #get PDE solution for acceleration 
74 
u_p1=(2.*uu_m1)+h*h*accel #calculate displacement 
75 
u_m1=u; u=u_p1 # shift values by 1 
76 
\end{python} 
77 

78 
Saving the data has been handled slightly differently in this example. The VTK 
79 
files generated can be quite large and take a significant amount of time to save 
80 
to the hard disk. To avoid doing this at every iteration a test is devised which 
81 
saves only at specific time intervals. 
82 

83 
To do this there are two new parameters in our script. 
84 
\begin{python} 
85 
# data recording times 
86 
rtime=0.0 # first time to record 
87 
rtime_inc=tend/20.0 # time increment to record 
88 
\end{python} 
89 
Currently the PDE solution will be saved to file $20$ times between the start of 
90 
the modelling and the final time step. With these parameters set, an if 
91 
statement is introduced to the time loop 
92 
\begin{python} 
93 
if (t >= rtime): 
94 
saveVTK(os.path.join(savepath,"ex08a.%05d.vtu"%n),displacement=length(u),\ 
95 
acceleration=length(accel),tensor=stress) 
96 
rtime=rtime+rtime_inc #increment data save time 
97 
\end{python} 
98 
\verb!t! is the time counter. Whenever the recording time \verb!rtime! is less 
99 
then \verb!t! the solution is saved and \verb!rtime! is incremented. This 
100 
limits the number of outputs and increases the speed of the solver. 
101 

102 
\section{Multithreading} 
103 
The wave equation solution can be quite demanding on CPU time. Enhancements can 
104 
be made by accessing multiple threads or cores on your computer. This does not 
105 
require any modification to the solution script and only comes into play when 
106 
\esc is called from the shell. To use multiple threads \esc is called using 
107 
the \verb!t! option with an integer argument for the number of threads 
108 
required. For example 
109 
\begin{verbatim} 
110 
$escript t 4 example08a.py 
111 
\end{verbatim} 
112 
would call the script in this section and solve it using 4 threads. 
113 

114 
The computation times on an increasing number of cores is outlined in 
115 
\autoref{tab:wpcores}. 
116 

117 
\begin{table}[ht] 
118 
\begin{center} 
119 
\caption{Computation times for an increasing number of cores.} 
120 
\label{tab:wpcores} 
121 
\begin{tabular}{ c  c } 
122 
\hline 
123 
Number of Cores & Time (s) \\ 
124 
\hline 
125 
1 & 691.0 \\ 
126 
2 & 400.0 \\ 
127 
3 & 305.0 \\ 
128 
4 & 328.0 \\ 
129 
5 & 323.0 \\ 
130 
6 & 292.0 \\ 
131 
7 & 282.0 \\ 
132 
8 & 445.0 \\ \hline 
133 
\end{tabular} 
134 
\end{center} 
135 
\end{table} 
136 

137 
\section{Vector source on the boundary} 
138 
\sslist{example08b.py} 
139 
For this particular example, we will introduce the source by applying a 
140 
displacement to the boundary during the initial time steps. The source will 
141 
again be 
142 
a radially propagating wave but due to the vector nature of the PDE used, a 
143 
direction will need to be applied to the source. 
144 

145 
The first step is to choose an amplitude and create the source as in the 
146 
previous chapter. 
147 
\begin{python} 
148 
U0=0.01 # amplitude of point source 
149 
# will introduce a spherical source at middle left of bottom face 
150 
xc=[ndx/2,0] 
151 

152 
############################################FIRST TIME STEPS AND SOURCE 
153 
# define small radius around point xc 
154 
src_length = 40; print "src_length = ",src_length 
155 
# set initial values for first two time steps with source terms 
156 
xb=FunctionOnBoundary(domain).getX() 
157 
y=source[0]*(cos(length(xxc)*3.1415/src_length)+1)*\ 
158 
whereNegative(length(xbsrc_length)) 
159 
src_dir=numpy.array([0.,1.]) # defines direction of point source as down 
160 
y=y*src_dir 
161 
\end{python} 
162 
where \verb xc is the source point on the boundary of the model. Note that 
163 
because the source is specifically located on the boundary, we have used the 
164 
\verb!FunctionOnBoundary! call to ensure the nodes used to define the source 
165 
are also located upon the boundary. These boundary nodes are passed to 
166 
source as \verb!xb!. The source direction is then defined as an $(x,y)$ array and multiplied by the 
167 
source function. The directional array must have a magnitude of $\left 1 
168 
\right $ otherwise the amplitude of the source will become modified. For this 
169 
example, the source is directed in the $y$ direction. 
170 
\begin{python} 
171 
src_dir=numpy.array([0.,1.]) # defines direction of point source as down 
172 
y=y*src_dir 
173 
\end{python} 
174 
The function can then be applied as a boundary condition by setting it equal to 
175 
$y$ in the general form. 
176 
\begin{python} 
177 
mypde.setValue(y=y) #set the source as a function on the boundary 
178 
\end{python} 
179 
The final step is to qualify the initial conditions. Due to the fact that we are 
180 
no longer using the source to define our initial condition to the model, we 
181 
must set the model state to zero for the first two time steps. 
182 
\begin{python} 
183 
# initial value of displacement at point source is constant (U0=0.01) 
184 
# for first two time steps 
185 
u=[0.0,0.0]*wherePositive(x) 
186 
u_m1=u 
187 
\end{python} 
188 

189 
If the source is time progressive, $y$ can be updated during the 
190 
iteration stage. This is covered in the following section. 
191 

192 
\begin{figure}[htp] 
193 
\centering 
194 
\subfigure[Example 08a at 0.025s ]{ 
195 
\includegraphics[width=3in]{figures/ex08pw50.png} 
196 
\label{fig:ex08pw50} 
197 
} 
198 
\subfigure[Example 08a at 0.175s ]{ 
199 
\includegraphics[width=3in]{figures/ex08pw350.png} 
200 
\label{fig:ex08pw350} 
201 
} \\ 
202 
\subfigure[Example 08a at 0.325s ]{ 
203 
\includegraphics[width=3in]{figures/ex08pw650.png} 
204 
\label{fig:ex08pw650} 
205 
} 
206 
\subfigure[Example 08a at 0.475s ]{ 
207 
\includegraphics[width=3in]{figures/ex08pw950.png} 
208 
\label{fig:ex08pw950} 
209 
} 
210 
\label{fig:ex08pw} 
211 
\caption{Results of Example 08 at various times.} 
212 
\end{figure} 
213 
\clearpage 
214 

215 
\section{Time variant source} 
216 

217 
\sslist{example08b.py} 
218 
Until this point, all of the wave propagation examples in this cookbook have 
219 
used impulsive sources which are smooth in space but not time. It is however, 
220 
advantageous to have a time smoothed source as it can reduce the temporal 
221 
frequency range and thus mitigate aliasing in the solution. 
222 

223 
It is quite 
224 
simple to implement a source which is smooth in time. In addition to the 
225 
original source function the only extra requirement is a time function. For 
226 
this example the time variant source will be the derivative of a Gaussian curve 
227 
defined by the required dominant frequency (\autoref{fig:tvsource}). 
228 
\begin{python} 
229 
#Creating the time function of the source. 
230 
dfeq=50 #Dominant Frequency 
231 
a = 2.0 * (np.pi * dfeq)**2.0 
232 
t0 = 5.0 / (2.0 * np.pi * dfeq) 
233 
srclength = 5. * t0 
234 
ls = int(srclength/h) 
235 
print 'source length',ls 
236 
source=np.zeros(ls,'float') # source array 
237 
ampmax=0 
238 
for it in range(0,ls): 
239 
t = it*h 
240 
tt = tt0 
241 
dum1 = np.exp(a * tt * tt) 
242 
source[it] = 2. * a * tt * dum1 
243 
if (abs(source[it]) > ampmax): 
244 
ampmax = abs(source[it]) 
245 
time[t]=t*h 
246 
\end{python} 
247 
\begin{figure}[ht] 
248 
\centering 
249 
\includegraphics[width=3in]{figures/source.png} 
250 
\caption{Time variant source with a dominant frequency of 50Hz.} 
251 
\label{fig:tvsource} 
252 
\end{figure} 
253 

254 
We then build the source and the first two time steps via; 
255 
\begin{python} 
256 
# set initial values for first two time steps with source terms 
257 
y=source[0] 
258 
*(cos(length(xxc)*3.1415/src_length)+1)*whereNegative(length(xxc)src_length) 
259 
src_dir=numpy.array([0.,1.]) # defines direction of point source as down 
260 
y=y*src_dir 
261 
mypde.setValue(y=y) #set the source as a function on the boundary 
262 
# initial value of displacement at point source is constant (U0=0.01) 
263 
# for first two time steps 
264 
u=[0.0,0.0]*whereNegative(x) 
265 
u_m1=u 
266 
\end{python} 
267 

268 
Finally, for the length of the source, we are required to update each new 
269 
solution in the iterative section of the solver. This is done via; 
270 
\begin{python} 
271 
# increment loop values 
272 
t=t+h; n=n+1 
273 
if (n < ls): 
274 
y=source[n]**(cos(length(xxc)*3.1415/src_length)+1)*\ 
275 
whereNegative(length(xxc)src_length) 
276 
y=y*src_dir; mypde.setValue(y=y) #set the source as a function on the 
277 
boundary 
278 
\end{python} 
279 

280 
\section{Absorbing Boundary Conditions} 
281 
To mitigate the effect of the boundary on the model, absorbing boundary 
282 
conditions can be introduced. These conditions effectively dampen the wave energy 
283 
as they approach the boundary and thus prevent that energy from being reflected. 
284 
This type of approach is typically used when a model is shrunk to decrease 
285 
computational requirements. In practise this applies to almost all models, 
286 
especially in earth sciences where the entire planet or a large enough 
287 
portional of it cannot be modelled efficiently when considering small scale 
288 
problems. It is impractical to calculate the solution for an infinite model and thus ABCs allow 
289 
us the create an approximate solution with small to zero boundary effects on a 
290 
model with a solvable size. 
291 

292 
To dampen the waves, the method of \citet{Cerjan1985} 
293 
where the solution and the stress are multiplied by a damping function defined 
294 
on $n$ nodes of the domain adjacent to the boundary, given by; 
295 
\begin{equation} 
296 
\gamma =\sqrt{\frac{ log( \gamma _{b} ) }{n^2}} 
297 
\end{equation} 
298 
\begin{equation} 
299 
y=e^{(\gamma x)^2} 
300 
\end{equation} 
301 
This is applied to the bounding 2050 pts of the model using the location 
302 
specifiers of \esc; 
303 
\begin{python} 
304 
# Define where the boundary decay will be applied. 
305 
bn=30. 
306 
bleft=xstep*bn; bright=mx(xstep*bn); bbot=my(ystep*bn) 
307 
# btop=ystep*bn # don't apply to force boundary!!! 
308 

309 
# locate these points in the domain 
310 
left=x[0]bleft; right=x[0]bright; bottom=x[1]bbot 
311 

312 
tgamma=0.98 # decay value for exponential function 
313 
def calc_gamma(G,npts): 
314 
func=np.sqrt(abs(1.*np.log(G)/(npts**2.))) 
315 
return func 
316 

317 
gleft = calc_gamma(tgamma,bleft) 
318 
gright = calc_gamma(tgamma,bleft) 
319 
gbottom= calc_gamma(tgamma,ystep*bn) 
320 

321 
print 'gamma', gleft,gright,gbottom 
322 

323 
# calculate decay functions 
324 
def abc_bfunc(gamma,loc,x,G): 
325 
func=exp(1.*(gamma*abs(locx))**2.) 
326 
return func 
327 

328 
fleft=abc_bfunc(gleft,bleft,x[0],tgamma) 
329 
fright=abc_bfunc(gright,bright,x[0],tgamma) 
330 
fbottom=abc_bfunc(gbottom,bbot,x[1],tgamma) 
331 
# apply these functions only where relevant 
332 
abcleft=fleft*whereNegative(left) 
333 
abcright=fright*wherePositive(right) 
334 
abcbottom=fbottom*wherePositive(bottom) 
335 
# make sure the inside of the abc is value 1 
336 
abcleft=abcleft+whereZero(abcleft) 
337 
abcright=abcright+whereZero(abcright) 
338 
abcbottom=abcbottom+whereZero(abcbottom) 
339 
# multiply the conditions together to get a smooth result 
340 
abc=abcleft*abcright*abcbottom 
341 
\end{python} 
342 
Note that the boundary conditions are not applied to the surface, as this is 
343 
effectively a free surface where normal reflections would be experienced. 
344 
Special conditions can be introduced at this surface if they are known. The 
345 
resulting boundary damping function can be viewed in 
346 
\autoref{fig:abconds}. 
347 

348 
\section{Second order Meshing} 
349 
For stiff problems like the wave equation it is often prudent to implement 
350 
second order meshing. This creates a more accurate mesh approximation with some 
351 
increased processing cost. To turn second order meshing on, the \verb!rectangle! 
352 
function accepts an \verb!order! keyword argument. 
353 
\begin{python} 
354 
domain=Rectangle(l0=mx,l1=my,n0=ndx, n1=ndy,order=2) # create the domain 
355 
\end{python} 
356 
Other pycad functions and objects have similar keyword arguments for higher 
357 
order meshing. 
358 

359 
Note that when implementing second order meshing, a smaller timestep is required 
360 
then for first order meshes as the second order essentially reduces the size of 
361 
the mesh by half. 
362 

363 
\begin{figure}[ht] 
364 
\centering 
365 
\includegraphics[width=5in]{figures/ex08babc.png} 
366 
\label{fig:abconds} 
367 
\caption{Absorbing boundary conditions for example08b.py} 
368 
\end{figure} 
369 

370 
\begin{figure}[htp] 
371 
\centering 
372 
\subfigure[Example 08b at 0.03s ]{ 
373 
\includegraphics[width=3in]{figures/ex08sw060.png} 
374 
\label{fig:ex08pw060} 
375 
} 
376 
\subfigure[Example 08b at 0.16s ]{ 
377 
\includegraphics[width=3in]{figures/ex08sw320.png} 
378 
\label{fig:ex08pw320} 
379 
} \\ 
380 
\subfigure[Example 08b at 0.33s ]{ 
381 
\includegraphics[width=3in]{figures/ex08sw660.png} 
382 
\label{fig:ex08pw660} 
383 
} 
384 
\subfigure[Example 08b at 0.44s ]{ 
385 
\includegraphics[width=3in]{figures/ex08sw880.png} 
386 
\label{fig:ex08pw880} 
387 
} 
388 
\label{fig:ex08pw} 
389 
\caption{Results of Example 08b at various times.} 
390 
\end{figure} 
391 
\clearpage 
392 

393 
\section{Pycad example} 
394 
\sslist{example08c.py} 
395 
To make the problem more interesting we will now introduce an interface to the 
396 
middle of the domain. In fact we will use the same domain as we did fora 
397 
different set of material properties on either side of the interface. 
398 

399 
\begin{figure}[ht] 
400 
\begin{center} 
401 
\includegraphics[width=5in]{figures/gmshexample08c.png} 
402 
\caption{Domain geometry for example08c.py showing line tangents.} 
403 
\label{fig:ex08cgeo} 
404 
\end{center} 
405 
\end{figure} 
406 

407 
It is simple enough to slightly modify the scripts of the previous sections to 
408 
accept this domain. Multiple material parameters must now be deined and assigned 
409 
to specific tagged areas. Again this is done via 
410 
\begin{python} 
411 
lam=Scalar(0,Function(domain)) 
412 
lam.setTaggedValue("top",lam1) 
413 
lam.setTaggedValue("bottom",lam2) 
414 
mu=Scalar(0,Function(domain)) 
415 
mu.setTaggedValue("top",mu1) 
416 
mu.setTaggedValue("bottom",mu2) 
417 
rho=Scalar(0,Function(domain)) 
418 
rho.setTaggedValue("top",rho1) 
419 
rho.setTaggedValue("bottom",rho2) 
420 
\end{python} 
421 
Don't forget that the source boundary must also be tagged and added so it can 
422 
be referenced 
423 
\begin{python} 
424 
# Add the subdomains and flux boundaries. 
425 
d.addItems(PropertySet("top",tblock),PropertySet("bottom",bblock),\ 
426 
PropertySet("linetop",l30)) 
427 
\end{python} 
428 
It is now possible to solve the script as in the previous examples 
429 
(\autoref{fig:ex08cres}). 
430 

431 
\begin{figure}[ht] 
432 
\centering 
433 
\includegraphics[width=4in]{figures/ex08c2601.png} 
434 
\caption{Modelling results of example08c.py at 0.2601 seconds. Notice the 
435 
refraction of the wave front about the boundary between the two layers.} 
436 
\label{fig:ex08cres} 
437 
\end{figure} 
438 

439 
