1 

2 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
3 
% 
4 
% Copyright (c) 20032009 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 

15 
\section{3D Wave Propagation} 
16 
\label{WAVE CHAP} 
17 

18 
In this next example we want to calculate the displacement field $u\hackscore{i}$ for any time $t>0$ by solving the wave equation: 
19 
\index{wave equation} 
20 
\begin{eqnarray}\label{WAVE general problem} 
21 
\rho u\hackscore{i,tt}  \sigma\hackscore{ij,j}=0 
22 
\end{eqnarray} 
23 
in a three dimensional block of length $L$ in $x\hackscore{0}$ 
24 
and $x\hackscore{1}$ direction and height $H$ 
25 
in $x\hackscore{2}$ direction. $\rho$ is the known density which may be a function of its location. 
26 
$\sigma\hackscore{ij}$ is the stress field \index{stress} which in case of an isotropic, linear elastic material is given by 
27 
\begin{eqnarray} \label{WAVE stress} 
28 
\sigma\hackscore{ij} & = & \lambda u\hackscore{k,k} \delta\hackscore{ij} + \mu ( u\hackscore{i,j} + u\hackscore{j,i}) 
29 
\end{eqnarray} 
30 
where $\lambda$ and $\mu$ are the Lame coefficients 
31 
\index{Lame coefficients} and $\delta\hackscore{ij}$ denotes the Kronecker symbol\index{Kronecker symbol}. 
32 
On the boundary the normal stress is given by 
33 
\begin{eqnarray} \label{WAVE natural} 
34 
\sigma\hackscore{ij}n\hackscore{j}=0 
35 
\end{eqnarray} 
36 
for all time $t>0$. 
37 

38 
\begin{figure}[t!] 
39 
\centerline{\includegraphics[angle=90,width=4.in]{figures/waveu}} 
40 
\caption{Input Displacement at Source Point ($\alpha=0.7$, $t\hackscore{0}=3$, $U\hackscore{0}=1$).} 
41 
\label{WAVE FIG 1.2} 
42 
\end{figure} 
43 

44 
\begin{figure}[t!] 
45 
\centerline{\includegraphics[angle=90,width=4.in]{figures/wavea}} 
46 
\caption{Input Acceleration at Source Point ($\alpha=0.7$, $t\hackscore{0}=3$, $U\hackscore{0}=1$).} 
47 
\label{WAVE FIG 1.1} 
48 
\end{figure} 
49 

50 
Here we are modelling a point source at the point $x\hackscore C$ in the $x\hackscore{0}$direction 
51 
which is a negative pulse of amplitude $U\hackscore{0}$ followed by the same 
52 
positive pulse. In mathematical terms we use 
53 
\begin{eqnarray} \label{WAVE source} 
54 
u\hackscore{0}(x\hackscore C,t)= U\hackscore{0} \; \sqrt{2} \frac{(tt\hackscore{0})}{\alpha} e^{\frac{1}{2}\frac{(tt\hackscore{0})^2}{\alpha^2}} \ 
55 
\end{eqnarray} 
56 
for all $t\ge0$ where $\alpha$ is the width of the pulse and $t\hackscore{0}$ is the time when 
57 
the pulse changes from negative to positive. In the simulations we will choose $\alpha=0.3$ and $t\hackscore{0}=2$, see Figure~\ref{WAVE FIG 1.2} 
58 
and will apply the source as a constraint in a in a sphere of small radius around the point 
59 
$x\hackscore C$. 
60 

61 
We use an explicit time integration scheme to calculate the displacement field $u$ at 
62 
certain time marks $t^{(n)}$ where $t^{(n)}=t^{(n1)}+h$ with time step size $h>0$. In the following the upper index ${(n)}$ refers to values at time $t^{(n)}$. We use the Verlet scheme \index{Verlet scheme} with constant time step size $h$ 
63 
which is defined by 
64 
\begin{eqnarray} \label{WAVE dyn 2} 
65 
u^{(n)}=2u^{(n1)}u^{(n2)} + h^2 a^{(n)} \\ 
66 
\end{eqnarray} 
67 
for all $n=2,3,\ldots$. It is designed to solve a system of equations of the form 
68 
\begin{eqnarray} \label{WAVE dyn 2b} 
69 
u\hackscore{,tt}=G(u) 
70 
\end{eqnarray} 
71 
where one sets $a^{(n)}=G(u^{(n1)})$. 
72 

73 
In our case $a^{(n)}$ is given by 
74 
\begin{eqnarray}\label{WAVE dyn 3} 
75 
\rho a^{(n)}\hackscore{i}=\sigma^{(n1)}\hackscore{ij,j} 
76 
\end{eqnarray} 
77 
and boundary conditions 
78 
\begin{eqnarray} \label{WAVE natural at n} 
79 
\sigma^{(n1)}\hackscore{ij}n\hackscore{j}=0 
80 
\end{eqnarray} 
81 
derived from \eqn{WAVE natural} where 
82 
\begin{eqnarray} \label{WAVE dyn 3a} 
83 
\sigma\hackscore{ij}^{(n1)} & = & \lambda u^{(n1)}\hackscore{k,k} \delta\hackscore{ij} + \mu ( u^{(n1)}\hackscore{i,j} + u^{(n1)}\hackscore{j,i}). 
84 
\end{eqnarray} 
85 
We also need to apply the constraint 
86 
\begin{eqnarray} \label{WAVE dyn 3aa} 
87 
a^{(n)}\hackscore{0}(x\hackscore C,t)= U\hackscore{0} 
88 
\frac{\sqrt(2.)}{\alpha^2} (4\frac{(tt\hackscore{0})^3}{\alpha^3}6\frac{tt\hackscore{0}}{\alpha})e^{\frac{1}{2}\frac{(tt\hackscore{0})^2}{\alpha^2}} 
89 
\end{eqnarray} 
90 
which is derived from equation~\ref{WAVE source} by calculating the second order time derivative, 
91 
see Figure~\ref{WAVE FIG 1.1}. Now we have converted our problem for displacement, $u^{(n)}$, into a problem for 
92 
acceleration, $a^{(n)}$, which now depends 
93 
on the solution at the previous two time steps, $u^{(n1)}$ and $u^{(n2)}$. 
94 

95 
In each time step we have to solve this problem to get the acceleration $a^{(n)}$, and we will 
96 
use the \LinearPDE class of the \linearPDEs to do so. The general form of the PDE defined through 
97 
the \LinearPDE class is discussed in \Sec{SEC LinearPDE}. The form which is relevant here is 
98 
\begin{equation}\label{WAVE dyn 100} 
99 
D\hackscore{ij} a^{(n)}\hackscore{j} =  X\hackscore{ij,j}\; . 
100 
\end{equation} 
101 
The natural boundary condition 
102 
\begin{equation}\label{WAVE dyn 101} 
103 
n\hackscore{j}X\hackscore{ij} =0 
104 
\end{equation} 
105 
is used. 
106 
With $u=a^{(n)}$ we can identify the values to be assigned to $D$ and $X$: 
107 
\begin{equation}\label{WAVE dyn 6} 
108 
\begin{array}{ll} 
109 
D\hackscore{ij}=\rho \delta\hackscore{ij}& 
110 
X\hackscore{ij}=\sigma^{(n1)}\hackscore{ij} \\ 
111 
\end{array} 
112 
\end{equation} 
113 
Moreover we need to define the location $r$ where the constraint~\ref{WAVE dyn 3aa} is applied. We will apply 
114 
the constraint on a small sphere of radius $R$ around $x\hackscore C$ (we will use 3p.c. of the width of the domain): 
115 
\begin{equation}\label{WAVE dyn 6 1} 
116 
q\hackscore{i}(x) = 
117 
\left\{ 
118 
\begin{array}{rc} 
119 
1 & \xx_c\\le R \\ 
120 
0 & \mbox{otherwise} 
121 
\end{array} 
122 
\right. 
123 
\end{equation} 
124 
The following script defines a the function \function{wavePropagation} which 
125 
implements the Verlet scheme to solve our wave propagation problem. 
126 
The argument \var{domain} which is a \Domain class object 
127 
defines the domain of the problem. \var{h} and \var{tend} are the time step size 
128 
and the end time of the simulation. \var{lam}, \var{mu} and 
129 
\var{rho} are material properties. 
130 
\begin{python} 
131 
def wavePropagation(domain,h,tend,lam,mu,rho, x_c, src_radius, U0): 
132 
# lists to collect displacement at point source which is returned to the caller 
133 
ts, u_pc0,u_pc1,u_pc2=[], [], [], [] 
134 

135 
x=domain.getX() 
136 
# ... open new PDE ... 
137 
mypde=LinearPDE(domain) 
138 
mypde.getSolverOptions().setSolverMethod(mypde.getSolverOptions().LUMPING) 
139 
kronecker=identity(mypde.getDim()) 
140 
dunit=numpy.array([1.,0.,0.]) # defines direction of point source 
141 
mypde.setValue(D=kronecker*rho, q=whereNegative(length(xxc)src_radius)*dunit) 
142 
# ... set initial values .... 
143 
n=0 
144 
# for first two time steps 
145 
u=Vector(0.,Solution(domain)) 
146 
u_last=Vector(0.,Solution(domain)) 
147 
t=0 
148 
# define the location of the point source 
149 
L=Locator(domain,xc) 
150 
# find potential at point source 
151 
u_pc=L.getValue(u) 
152 
print "u at point charge=",u_pc 
153 
# open file to save displacement at point source 
154 
u_pc_data=FileWriter('./data/U_pc.out') 
155 
ts.append(t); u_pc0.append(u_pc[0]), u_pc1.append(u_pc[1]), u_pc2.append(u_pc[2]) 
156 

157 
while t<tend: 
158 
t+=h 
159 
# ... get current stress .... 
160 
g=grad(u) 
161 
stress=lam*trace(g)*kronecker+mu*(g+transpose(g)) 
162 
# ... get new acceleration .... 
163 
amplitude=U0*(4*(tt0)**3/alpha**36*(tt0)/alpha)*sqrt(2.)/alpha**2 \ 
164 
*exp(1./2.(tt0)**2/alpha**2) 
165 
mypde.setValue(X=stress, r=dunit*amplitude) 
166 
a=mypde.getSolution() 
167 
# ... get new displacement ... 
168 
u_new=2*uu_last+h**2*a 
169 
# ... shift displacements .... 
170 
u_last=u 
171 
u=u_new 
172 
n+=1 
173 
print n,"th time step t ",t 
174 
u_pc=L.getValue(u) 
175 
print "u at point charge=",u_pc 
176 
# save displacements at point source to file for t > 0 
177 
ts.append(t); u_pc0.append(u_pc[0]), u_pc1.append(u_pc[1]), \ 
178 
u_pc2.append(u_pc[2]) 
179 

180 
# ... save current acceleration in units of gravity and displacements 
181 
if n==1 or n%10==0: saveVTK("./data/usoln.%i.vtu"%(n/10), \ 
182 
acceleration=length(a)/9.81, 
183 
displacement = length(u), \ 
184 
tensor = stress, Ux = u[0] ) 
185 

186 
return ts, u_pc0, u_pc1, u_pc2 
187 
\end{python} 
188 
Notice that 
189 
all coefficients of the PDE which are independent of time $t$ are set outside the \code{while} 
190 
loop. This is very efficient since it allows the \LinearPDE class to reuse information as much as possible 
191 
when iterating over time. 
192 

193 
The statement 
194 
\begin{python} 
195 
mypde.getSolverOptions().setSolverMethod(mypde.getSolverOptions().LUMPING) 
196 
\end{python} 
197 
switches on the use of an aggressive approximation of the PDE operator as a diagonal matrix 
198 
formed from the coefficient \var{D}. 
199 
The approximation allows, at the cost of 
200 
additional error, very fast 
201 
solution of the PDE. When using lumping the presence of \var{A}, \var{B} or \var{C} will produce wrong results. 
202 

203 
There are a few new \escript functions in this example: 
204 
\code{grad(u)} returns the gradient $u\hackscore{i,j}$ of $u$ (in fact \var{grad(g)[i,j]}=$=u\hackscore{i,j}$). 
205 
There are restrictions on the argument of the \function{grad} function, for instance 
206 
the statement \code{grad(grad(u))} will raise an exception. 
207 
\code{trace(g)} returns the sum of the main diagonal elements \var{g[k,k]} of \var{g} 
208 
and \code{transpose(g)} returns the matrix transpose of \var{g} (ie. $\var{transpose(g)[i,j]}=\var{g[j,i]}$). 
209 

210 
We initialize the values of \code{u} and \code{u_last} to be zero. It is important 
211 
to initialize both with the \SolutionFS \FunctionSpace as they have to be seen as solutions of PDEs from previous time steps. In fact, the \function{grad} does not accept arguments with a certain \FunctionSpace, for more details see \Sec{SEC ESCRIPT DATA}. 
212 

213 
The \class{Locator} is designed to extract values at a given location (in this case $x\hackscore C$) from functions such as the displacement vector \code{u}. Typically the \class{Locator} is used in the following form: 
214 
\begin{python} 
215 
L=Locator(domain,xc) 
216 
u=... 
217 
u_pc=L.getValue(u) 
218 
\end{python} 
219 
The return value \code{u_pc} is the value of \code{u} at the location \code{xc}\footnote{In fact the finite element node which is closest to the given position. The usage of \class{Locator} is MPI save.}. The values 
220 
are collected in the lists \var{u_pc0}, \var{u_pc1} and \var{u_pc2} together with the 
221 
corresponding time marker in \var{ts}. The values are handed back to the caller. Later we will show to ways to 
222 
access these data. 
223 

224 
One of the big advantages of the Verlet scheme is the fact that the problem to be solved 
225 
in each time step is very simple and does not involve any spatial derivatives (which is what allows us to use lumping in this simulation). 
226 
The problem becomes so simple because we use the stress from the last time step rather then the stress which is 
227 
actually present at the current time step. Schemes using this approach are called an explicit time integration 
228 
schemes \index{explicit scheme} \index{time integration!explicit}. The 
229 
backward Euler scheme we have used in \Chap{DIFFUSION CHAP} is 
230 
an example of an implicit scheme 
231 
\index{implicit scheme} \index{time integration!implicit}. In this case one uses the actual status of 
232 
each variable at a particular time rather then values from previous time steps. This will lead to a problem which is 
233 
more expensive to solve, in particular for nonlinear problems. 
234 
Although 
235 
explicit time integration schemes are cheap to finalize a single time step, they need significantly smaller time 
236 
steps then implicit schemes and can suffer from stability problems. Therefore they need a 
237 
very careful selection of the time step size $h$. 
238 

239 
An easy, heuristic way of choosing an appropriate 
240 
time step size is the Courant condition \index{Courant condition} \index{explicit scheme!Courant condition} 
241 
which says that within a time step a information should not travel further than a cell used in the 
242 
discretization scheme. In the case of the wave equation the velocity of a (p) wave is given as 
243 
$\sqrt{\frac{\lambda+2\mu}{\rho}}$ so one should choose $h$ from 
244 
\begin{eqnarray}\label{WAVE dyn 66} 
245 
h= \frac{1}{5} \sqrt{\frac{\rho}{\lambda+2\mu}} \Delta x 
246 
\end{eqnarray} 
247 
where $\Delta x$ is the cell diameter. The factor $\frac{1}{5}$ is a safety factor considering the heuristics of 
248 
the formula. 
249 

250 
\begin{figure}[t] 
251 
\begin{center} 
252 
\includegraphics[width=2in]{figures/Wave11} 
253 
\includegraphics[width=2in]{figures/Wave22} 
254 
\includegraphics[width=2in]{figures/Wave28} 
255 
\includegraphics[width=2in]{figures/Wave32} 
256 
\includegraphics[width=2in]{figures/Wave36} 
257 
\end{center} 
258 
\caption{Selected time steps ($n = 11, 22, 32, 36$) of a wave propagation over a $10\mbox{km} \times 10\mbox{km} \times 3.125\mbox{km}$ block 
259 
from a point source initially at $(5\mbox{km},5\mbox{km},0)$ 
260 
with time step size $h=0.02083$. Color represents the displacement. 
261 
Here the view is oriented onto the bottom face. 
262 
\label{WAVE FIG 2}} 
263 
\end{figure} 
264 

265 
The following script uses the \code{wavePropagation} function to solve the 
266 
wave equation for a point source located at the bottom face of a block. The width of the block in 
267 
each direction on the bottom face is $10\mbox{km}$ ($x\hackscore 0$ and $x\hackscore 1$ directions ie. \code{l0} and \code{l1}). 
268 
The \code{ne} gives the number of elements in $x\hackscore{0}$ and $x\hackscore 1$ directions. 
269 
The depth of the block is aligned with the $x\hackscore{2}$direction. 
270 
The depth (\code{l2}) of the block in 
271 
the $x\hackscore{2}$direction is chosen so that there are 10 elements and the magnitude of the 
272 
of the depth is chosen such that the elements 
273 
become cubic. We chose 10 for the number of elements in $x\hackscore{2}$direction so that the 
274 
computation would be faster. \code{Brick($n\hackscore 0, n\hackscore 1, n\hackscore 2,l\hackscore 0,l\hackscore 1,l\hackscore 2$)} is an \finley function which creates a rectangular mesh 
275 
with $n\hackscore 0 \times n\hackscore 1 \times n\hackscore2$ elements over the brick $[0,l\hackscore 0] \times [0,l\hackscore 1] \times [0,l\hackscore 2]$. 
276 
\begin{python} 
277 
from esys.finley import Brick 
278 
ne=32 # number of cells in x_0 and x_1 directions 
279 
width=10000. # length in x_0 and x_1 directions 
280 
lam=3.462e9 
281 
mu=3.462e9 
282 
rho=1154. 
283 
tend=60 
284 
U0=1. # amplitude of point source 
285 
# spherical source at middle of bottom face 
286 
xc=[width/2.,width/2.,0.] 
287 
# define small radius around point xc 
288 
src_radius = 0.03*width 
289 
print "src_radius = ",src_radius 
290 
mydomain=Brick(ne,ne,10,l0=width,l1=width,l2=10.*width/32.) 
291 
h=(1./5.)*inf(sqrt(rho/(lam+2*mu))*inf(domain.getSize()) 
292 
print "time step size = ",h 
293 
ts, u_pc0, u_pc1, u_pc2 = \ 
294 
wavePropagation(mydomain,h,tend,lam,mu,rho,xc, src_radius, U0) 
295 
\end{python} 
296 
The \function{domain.getSize()} return the local element size $\Delta x$. The 
297 
\function{inf} makes sure that the Courant condition~\ref{WAVE dyn 66} olds everywhere in the domain. 
298 

299 
The script is available as \file{wave.py} in the \ExampleDirectory \index{scripts!\file{wave.py}}. 
300 
To visualize the results from the data directory: 
301 
\begin{verbatim} 
302 
mayavi2 d usoln.1.vtu m SurfaceMap & 
303 
\end{verbatim} 
304 
You can rotate this figure by clicking on it with the mouse and moving it around. 
305 
Again use \code{Configure Data} to move backwards and forward in time, and 
306 
also to choose the results (acceleration, displacement or $u\hackscore x$) by using \code{Select Scalar}. Figure \ref{WAVE FIG 2} shows the results for the displacement at various time steps. 
307 

308 
\begin{figure}[t!] 
309 
\centerline{\includegraphics[width=4.in, angle=90]{figures/WavePC}} 
310 
\caption{Amplitude at Point source from the Simulation.} 
311 
\label{WAVE FIG 1} 
312 
\end{figure} 
313 

314 
It remains to show some possibilities to inspect the collected data \var{u_pc0}, \var{u_pc1} and \var{u_pc2}. 
315 
One way is to write the data to a file and then use an external package such as \gnuplot, excel or OpenOffice to read the data for further analysis. The following code shows one possible form to write the data to the 
316 
file \file{./data/U_pc.out}: 
317 
\begin{python} 
318 
u_pc_data=FileWriter('./data/U_pc.out') 
319 
for i in xrange(len(ts)) : 
320 
u_pc_data.write("%f %f %f %f\n"%(ts[i],u_pc0[i],u_pc1[i],u_pc2[i])) 
321 
u_pc_data.close() 
322 
\end{python} 
323 
The \code{U_pc.out} stores 4 columns of data: $t,u\hackscore x,u\hackscore y,u\hackscore z$ 
324 
respectively, where $u\hackscore x,u\hackscore y,u\hackscore z$ are the $x\hackscore 0,x\hackscore 1,x\hackscore 2$ components of the displacement vector $u$ at the point source. These can be 
325 
plotted easily using any plotting package. In \gnuplot the command: 
326 
\begin{verbatim} 
327 
plot 'U_pc.out' u 1:2 title 'U_x' w l lw 2, 'U_pc.out' u 1:3 title 'U_y' w l lw 2, 
328 
'U_pc.out' u 1:4 title 'U_z' w l lw 2 
329 
\end{verbatim} 
330 
will reproduce Figure~\ref{WAVE FIG 1} (As expected this is identical to the input signal shown in Figure~\ref{WAVE FIG 1.2})2. It is pointed out that we are not using the 
331 
standard \PYTHON \function{open} to write to the file \code{U_pc.out} as it is not safe 
332 
when running \escript under MPI, see chapter~\ref{EXECUTION} for more details. 
333 

334 
Alternatively, one can implement plotting the results at run time rather than in a postprocessing step. This avoids 
335 
the generation of an intermediate data file. In {\it escript} the preferred way of creating 2D plots of 
336 
time dependent data is \MATPLOTLIB. The following script creates the plot and writes it into the 
337 
file \file{u_pc.png} in the PNG image format: 
338 
\begin{python} 
339 
import matplotlib.pyplot as plt 
340 
if getMPIRankWorld() == 0: 
341 
plt.title("Displacement at Point Source") 
342 
plt.plot(ts, u_pc0, '', label="x_0", linewidth=1) 
343 
plt.plot(ts, u_pc1, '', label="x_1", linewidth=1) 
344 
plt.plot(ts, u_pc2, '', label="x_2", linewidth=1) 
345 
plt.xlabel('time') 
346 
plt.ylabel('displacement') 
347 
plt.legend() 
348 
plt.savefig('u_pc.png', format='png') 
349 
\end{python} 
350 
You can add the \function{plt.show()} to create a interactive browser window. Please not that 
351 
through the \code{getMPIRankWorld() == 0} statement the plot is generated on one processor only (in this case the rank 0 processor) when run under MPI. 
352 

353 
Both options for processing the point source data are include in the example file \file{wave.py}. There other options available to process these data in particular through the \SCIPY package , eg Fourier transformations. It is beyond the scope of this users guide to document the usage of \SCIPY for time series analysis but is highly recommended that users use \SCIPY to any further data processing. 
354 
