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 

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

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

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

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

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

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

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

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

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

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

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

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

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

214 
\section{Time variant source} 
215 

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

438 
