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\hackscore{i}}{\partial t^2}  \frac{\partial 
26 
\sigma\hackscore{ij}}{\partial x\hackscore{j}} = 0 
27 
\end{equation} 
28 
where $\sigma$ is the stress given by 
29 
\begin{equation} \label{eqn:sigw} 
30 
\sigma \hackscore{ij} = \lambda u\hackscore{k,k} \delta\hackscore{ij} + \mu ( 
31 
u\hackscore{i,j} + u\hackscore{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\hackscore{i}}{\partial t^2}$ we can derive the 
38 
displacement solution. Using $a$ we can see that \autoref{eqn:wav} becomes 
39 
\begin{equation} \label{eqn:wava} 
40 
\rho a\hackscore{i}  \frac{\partial 
41 
\sigma\hackscore{ij}}{\partial x\hackscore{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. 
45 

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

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

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

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

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

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

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

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

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

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

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

213 
\section{Time variant source} 
214 

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

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

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

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

278 
\section{Absorbing Boundary Conditions} 
279 
To mitigate the effect of the boundary on the model, absorbing boundary 
280 
conditions can be introduced. These conditions effectively dampen the wave 
281 
energy as they approach the bounday and thus prevent that energy from being 
282 
reflected. This type of approach is used typically when a model only represents 
283 
a small portion of the entire model, which in reality may have infinite bounds. 
284 
It is inpractical to calculate the solution for an infinite model and thus ABCs 
285 
allow us the create an approximate solution with small to zero boundary effects 
286 
on a model with a solvable size. 
287 

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

305 
# locate these points in the domain 
306 
left=x[0]bleft; right=x[0]bright; bottom=x[1]bbot 
307 

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

313 
gleft = calc_gamma(tgamma,bleft) 
314 
gright = calc_gamma(tgamma,bleft) 
315 
gbottom= calc_gamma(tgamma,ystep*bn) 
316 

317 
print 'gamma', gleft,gright,gbottom 
318 

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

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

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

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

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

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

389 
\section{Pycad example} 
390 
\sslist{example08c.py} 
391 
To make the problem more interesting we will now introduce an interface to the 
392 
middle of the domain. Infact we will use the same domain as we did for heat flux 
393 
in \autoref{CHAP HEAT 2}. The domain contains a syncline with two set of 
394 
material properties on either side of the interface. 
395 

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

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

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

435 
