1 

2 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
3 
% 
4 
% Copyright (c) 20032008 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 
At initial time $t=0$ the displacement 
39 
$u\hackscore{i}$ and the velocity $u\hackscore{i,t}$ are given: 
40 
\begin{eqnarray} \label{WAVE initial} 
41 
u\hackscore{i}(0,x)= \left\{\begin{array}{cc} U\hackscore{0} & \mbox{for $x$ at point charge, $x\hackscore{C}$} \\ 0 & \mbox{elsewhere}\end{array} \right. & \mbox{ and } & u\hackscore{i,t}(0,x)=0 \; 
42 
\end{eqnarray} 
43 
for all $x$ in the domain. 
44 

45 
Here we are modelling a point source at the point $x\hackscore C$, in the numerical solution we 
46 
set the initial displacement to be $U\hackscore 0$ in a sphere of small radius around the point 
47 
$x\hackscore C$. 
48 

49 
We use an explicit time integration scheme to calculate the displacement field $u$ at 
50 
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$ 
51 
which is defined by 
52 
\begin{eqnarray} \label{WAVE dyn 2} 
53 
u^{(n)}=2u^{(n1)}u^{(n2)} + h^2 a^{(n)} \\ 
54 
\end{eqnarray} 
55 
for all $n=2,3,\ldots$. It is designed to solve a system of equations of the form 
56 
\begin{eqnarray} \label{WAVE dyn 2b} 
57 
u\hackscore{,tt}=G(u) 
58 
\end{eqnarray} 
59 
where one sets $a^{(n)}=G(u^{(n1)})$. 
60 

61 
In our case $a^{(n)}$ is given by 
62 
\begin{eqnarray}\label{WAVE dyn 3} 
63 
\rho a^{(n)}\hackscore{i}=\sigma^{(n1)}\hackscore{ij,j} 
64 
\end{eqnarray} 
65 
and boundary conditions 
66 
\begin{eqnarray} \label{WAVE natural at n} 
67 
\sigma^{(n1)}\hackscore{ij}n\hackscore{j}=0 
68 
\end{eqnarray} 
69 
derived from \eqn{WAVE natural} where 
70 
\begin{eqnarray} \label{WAVE dyn 3a} 
71 
\sigma\hackscore{ij}^{(n1)} & = & \lambda u^{(n1)}\hackscore{k,k} \delta\hackscore{ij} + \mu ( u^{(n1)}\hackscore{i,j} + u^{(n1)}\hackscore{j,i}). 
72 
\end{eqnarray} 
73 
Now we have converted our problem for displacement, $u^{(n)}$, into a problem for 
74 
acceleration, $a^(n)$, which now depends 
75 
on the solution at the previous two time steps, $u^{(n1)}$ and $u^{(n2)}$. 
76 

77 
In each time step we have to solve this problem to get the acceleration $a^{(n)}$, and we will 
78 
use the \LinearPDE class of the \linearPDEs to do so. The general form of the PDE defined through 
79 
the \LinearPDE class is discussed in \Sec{SEC LinearPDE}. The form which is relevant here is 
80 
\begin{equation}\label{WAVE dyn 100} 
81 
D\hackscore{ij} a^{(n)}\hackscore{j} =  X\hackscore{ij,j}\; . 
82 
\end{equation} 
83 
The natural boundary condition 
84 
\begin{equation}\label{WAVE dyn 101} 
85 
n\hackscore{j}X\hackscore{ij} =0 
86 
\end{equation} 
87 
is used. 
88 

89 
With $u=a^{(n)}$ we can identify the values to be assigned to $D$ and $X$: 
90 
\begin{equation}\label{WAVE dyn 6} 
91 
\begin{array}{ll} 
92 
D\hackscore{ij}=\rho \delta\hackscore{ij}& 
93 
X\hackscore{ij}=\sigma^{(n1)}\hackscore{ij} \\ 
94 
\end{array} 
95 
\end{equation} 
96 

97 

98 
%Remember that \var{y.whereZero()} returns a function which is $1$ for those points $x$ where $y$ is zero 
99 
%(up to a certain tolerance, see \Sec{SEC ESCRIPT DATA}) and $0$ for those points where $y$ is not equal to zero. 
100 

101 
The following script defines a the function \function{wavePropagation} which 
102 
implements the Verlet scheme to solve our wave propagation problem. 
103 
The argument \var{domain} which is a \Domain class object 
104 
defines the domain of the problem. \var{h} and \var{tend} are the time step size 
105 
and the end time of the simulation. \var{lam}, \var{mu} and 
106 
\var{rho} are material properties. 
107 
\begin{python} 
108 
from esys.linearPDEs import LinearPDE 
109 
from numarray import identity,zeros,ones 
110 
from esys.escript import * 
111 
from esys.escript.pdetools import Locator 
112 
def wavePropagation(domain,h,tend,lam,mu,rho,U0): 
113 
x=domain.getX() 
114 
# ... open new PDE ... 
115 
mypde=LinearPDE(domain) 
116 
mypde.setSolverMethod(LinearPDE.LUMPING) 
117 
kronecker=identity(mypde.getDim()) 
118 
# spherical source at middle of bottom face 
119 
xc=[width/2.,width/2.,0.] 
120 
# define small radius around point xc 
121 
# Lsup(x) returns the maximum value of the argument x 
122 
src_radius = 0.1*Lsup(domain.getSize()) 
123 
print "src_radius = ",src_radius 
124 
dunit=numarray.array([1.,0.,0.]) # defines direction of point source 
125 
mypde.setValue(D=kronecker*rho) 
126 
# ... set initial values .... 
127 
n=0 
128 
# initial value of displacement at point source is constant (U0=0.01) 
129 
# for first two time steps 
130 
u=U0*whereNegative(length(xxc)src_radius)*dunit 
131 
u_last=U0*whereNegative(length(xxc)src_radius)*dunit 
132 
t=0 
133 
# define the location of the point source 
134 
L=Locator(domain,xc) 
135 
# find potential at point source 
136 
u_pc=L.getValue(u) 
137 
print "u at point charge=",u_pc 
138 
u_pc_x = u_pc[0] 
139 
u_pc_y = u_pc[1] 
140 
u_pc_z = u_pc[2] 
141 
# open file to save displacement at point source 
142 
u_pc_data=open('./data/U_pc.out','w') 
143 
u_pc_data.write("%f %f %f %f\n"%(t,u_pc_x,u_pc_y,u_pc_z)) 
144 
while t<tend: 
145 
# ... get current stress .... 
146 
g=grad(u) 
147 
stress=lam*trace(g)*kronecker+mu*(g+transpose(g)) 
148 
# ... get new acceleration .... 
149 
mypde.setValue(X=stress) 
150 
a=mypde.getSolution() 
151 
# ... get new displacement ... 
152 
u_new=2*uu_last+h**2*a 
153 
# ... shift displacements .... 
154 
u_last=u 
155 
u=u_new 
156 
t+=h 
157 
n+=1 
158 
print n,"th time step t ",t 
159 
L=Locator(domain,xc) 
160 
u_pc=L.getValue(u) 
161 
print "u at point charge=",u_pc 
162 
u_pc_x=u_pc[0] 
163 
u_pc_y=u_pc[1] 
164 
u_pc_z=u_pc[2] 
165 
# save displacements at point source to file for t > 0 
166 
u_pc_data.write("%f %f %f %f\n"%(t,u_pc_x,u_pc_y,u_pc_z)) 
167 
# ... save current acceleration in units of gravity and displacements 
168 
if n==1 or n%10==0: saveVTK("./data/usoln.%i.vtu"%(n/10),acceleration=length(a)/9.81, 
169 
displacement = length(u), Ux = u[0] ) 
170 
u_pc_data.close() 
171 
\end{python} 
172 
Notice that 
173 
all coefficients of the PDE which are independent of time $t$ are set outside the \code{while} 
174 
loop. This is very efficient since it allows the \LinearPDE class to reuse information as much as possible 
175 
when iterating over time. 
176 

177 
there are a few new \escript functions in this example: 
178 
\code{grad(u)} returns the gradient $u\hackscore{i,j}$ of $u$ (in fact \var{grad(g)[i,j]}=$=u\hackscore{i,j}$). 
179 
There are restrictions on the argument of the \function{grad} function, for instance 
180 
the statement \code{grad(grad(u))} will raise an exception. 
181 
\code{trace(g)} returns the sum of the main diagonal elements \var{g[k,k]} of \var{g} 
182 
and \code{transpose(g)} returns the matrix transpose of \var{g} (ie. $\var{transpose(g)[i,j]}=\var{g[j,i]}$). 
183 

184 
We initialize the values of \code{u} and \code{u_last} to be $U\hackscore 0$ for a small 
185 
sphere of radius, \code{src_radius} around the point source located at $x\hackscore C$. 
186 
These satisfy the initial conditions defined in \eqn{WAVE initial}. 
187 

188 
The output \code{U_pc.out} and \code{vtu} files are saved in a directory called \code{data}. 
189 
The \code{U_pc.out} stores 4 columns of data: $t,u\hackscore x,u\hackscore y,u\hackscore z$ 
190 
respectively, where $u\hackscore x,u\hackscore y,u\hackscore z$ are the $x\hackscore 0,x\hackscore 1,x\hackscore 2$ components of 
191 
the displacement vector $u$ at the point source. These can be 
192 
plotted easily using any plotting package. In gnuplot the command: 
193 
\code{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, 
194 
'U_pc.out' u 1:4 title 'U_z' w l lw 2} will reproduce Figure \ref{WAVE FIG 1}. 
195 
%\begin{python} 
196 
%u=Vector(0.,ContinuousFunction(domain)) 
197 
%\end{python} 
198 
%declares \var{u} as a vectorvalued function of its coordinates 
199 
%in the domain (i.e. with two or three components in the case of 
200 
%a two or three dimensional domain, respectively). The first argument defines the value of the function which is equal 
201 
%to $0$ for all vectot components. The second argument \var{ContinuousFunction(domain)} 
202 
%specifies the `smoothness` of the function. In our case we want the displacement field to be 
203 
%a continuous function over the \Domain \var{domain}. On other option would be 
204 
%\var{Function(domain)} in which case continuity is not guaranteed. In fact the 
205 
%argument of the \function{grad} function has to be of the typ \var{ContinuousFunction(domain)} while 
206 
%returned function is of the var{Function(domain)} type. For more details see \Sec{SEC ESCRIPT DATA}. 
207 

208 
The statement \code{myPDE.setLumpingOn()} 
209 
switches on the use of an aggressive approximation of the PDE operator as a diagonal matrix 
210 
formed from the coefficient \var{D}. 
211 
The approximation allows, at the cost of 
212 
additional error, very fast 
213 
solution of the PDE. When using \code{setLumpingOn} the presence of \var{A}, \var{B} or \var 
214 
{C} will produce wrong results. 
215 

216 

217 
One of the big advantages of the Verlet scheme is the fact that the problem to be solved 
218 
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). 
219 
The problem becomes so simple because we use the stress from the last time step rather then the stress which is 
220 
actually present at the current time step. Schemes using this approach are called an explicit time integration 
221 
schemes \index{explicit scheme} \index{time integration!explicit}. The 
222 
backward Euler scheme we have used in \Chap{DIFFUSION CHAP} is 
223 
an example of an implicit scheme 
224 
\index{implicit scheme} \index{time integration!implicit}. In this case one uses the actual status of 
225 
each variable at a particular time rather then values from previous time steps. This will lead to a problem which is 
226 
more expensive to solve, in particular for nonlinear problems. 
227 
Although 
228 
explicit time integration schemes are cheap to finalize a single time step, they need significantly smaller time 
229 
steps then implicit schemes and can suffer from stability problems. Therefore they need a 
230 
very careful selection of the time step size $h$. 
231 

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

243 
The following script uses the \code{wavePropagation} function to solve the 
244 
wave equation for a point source located at the bottom face of a block. The width of the block in 
245 
each direction on the bottom face is $10\mbox{km}$ ($x\hackscore 0$ and $x\hackscore 1$ directions ie. \code{l0} and \code{l1}). 
246 
\code{ne} gives the number of elements in $x\hackscore{0}$ and $x\hackscore 1$ directions. 
247 
The depth of the block is aligned with the $x\hackscore{2}$direction. 
248 
The depth (\code{l2}) of the block in 
249 
the $x\hackscore{2}$direction is chosen so that there are 10 elements and the magnitude of the 
250 
of the depth is chosen such that the elements 
251 
become cubic. We chose 10 for the number of elements in $x\hackscore{2}$direction so that the 
252 
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 
253 
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]$. 
254 
\begin{python} 
255 
from esys.finley import Brick 
256 
ne=32 # number of cells in x_0 and x_1 directions 
257 
width=10000. # length in x_0 and x_1 directions 
258 
lam=3.462e9 
259 
mu=3.462e9 
260 
rho=1154. 
261 
tend=60 
262 
h=(1./5.)*sqrt(rho/(lam+2*mu))*(width/ne) 
263 
print "time step size = ",h 
264 
U0=0.01 # amplitude of point source 
265 

266 
mydomain=Brick(ne,ne,10,l0=width,l1=width,l2=10.*width/32.) 
267 
wavePropagation(mydomain,h,tend,lam,mu,rho,U0) 
268 
\end{python} 
269 
The script is available as \file{wave.py} in the \ExampleDirectory \index{scripts!\file{wave.py}}. Before running the code make sure you have created a directory called \code{data} in the current 
270 
working directory. 
271 
To visualize the results from the data directory: 
272 
\begin{verbatim} mayavi d usoln.1.vtu m SurfaceMap &\end{verbatim} You can rotate this figure by clicking on it with the mouse and moving it around. 
273 
Again use \code{Configure Data} to move backwards and forward in time, and 
274 
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. 
275 

276 
\begin{figure}{t!} 
277 
\centerline{\includegraphics[width=3.in,angle=270]{figures/WavePC.eps}} 
278 
\caption{Amplitude at Point source} 
279 
\label{WAVE FIG 1} 
280 
\end{figure} 
281 

282 
\begin{figure}{t} 
283 
\begin{center} 
284 
\includegraphics[width=2in,angle=270]{figures/Wavet0.eps} 
285 
\includegraphics[width=2in,angle=270]{figures/Wavet1.eps} 
286 
\includegraphics[width=2in,angle=270]{figures/Wavet3.eps} 
287 
\includegraphics[width=2in,angle=270]{figures/Wavet10.eps} 
288 
\includegraphics[width=2in,angle=270]{figures/Wavet30.eps} 
289 
\includegraphics[width=2in,angle=270]{figures/Wavet288.eps} 
290 
\end{center} 
291 
\caption{Selected time steps ($n = 0,1,30,100,300$ and 2880) of a wave propagation over a $10\mbox{km} \times 10\mbox{km} \times 3.125\mbox{km}$ block 
292 
from a point source initially at $(5\mbox{km},5\mbox{km},0)$ 
293 
with time step size $h=0.02083$. Color represents the displacement. 
294 
Here the view is oriented onto the bottom face. 
295 
\label{WAVE FIG 2}} 
296 
\end{figure} 