1 
% $Id$ 
2 
% 
3 
% Copyright © 2006 by ACcESS MNRF 
4 
% \url{http://www.access.edu.au 
5 
% Primary Business: Queensland, Australia. 
6 
% Licensed under the Open Software License version 3.0 
7 
% http://www.opensource.org/license/osl3.0.php 
8 
% 
9 
\section{3D Wave Propagation} 
10 
\label{WAVE CHAP} 
11 

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

32 
At initial time $t=0$ the displacement 
33 
$u\hackscore{i}$ and the velocity $u\hackscore{i,t}$ are given: 
34 
\begin{eqnarray} \label{WAVE initial} 
35 
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 \; 
36 
\end{eqnarray} 
37 
for all $x$ in the domain. 
38 

39 
Here we are modelling a point source at the point $x\hackscore C$, in the numerical solution we 
40 
set the initial displacement to be $U\hackscore 0$ in a sphere of small radius around the point 
41 
$x\hackscore C$. 
42 

43 
We use an explicit timeintegration scheme to calculate the displacement field $u$ at 
44 
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$ 
45 
which is defined by 
46 
\begin{eqnarray} \label{WAVE dyn 2} 
47 
u^{(n)}=2u^{(n1)}u^{(n2)} + h^2 a^{(n)} \\ 
48 
\end{eqnarray} 
49 
for all $n=2,3,\ldots$. It is designed to solve a system of equations of the form 
50 
\begin{eqnarray} \label{WAVE dyn 2b} 
51 
u\hackscore{,tt}=G(u) 
52 
\end{eqnarray} 
53 
where one sets $a^{(n)}=G(u^{(n1)})$, see \Ref{Mora94}. 
54 

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

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

83 
With $u=a^{(n)}$ we can identify the values to be assigned to $D$, $X$: 
84 
\begin{equation}\label{WAVE dyn 6} 
85 
\begin{array}{ll} 
86 
D\hackscore{ij}=\rho \delta\hackscore{ij}& 
87 
X\hackscore{ij}=\sigma^{(n1)}\hackscore{ij} \\ 
88 
\end{array} 
89 
\end{equation} 
90 

91 

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

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

171 
We have seen most of the functions in previous examples but there are some new functions here: 
172 
\code{grad(u)} returns the gradient $u\hackscore{i,j}$ of $u$ (in fact \var{grad(g)[i,j]}=$=u\hackscore{i,j}$). 
173 
It is pointed out that in general there are restrictions on the argument of the \function{grad} function, for instance 
174 
for \finley the statement \code{grad(grad(u))} will raise an exception. 
175 
\code{trace(g)} returns the sum of the main diagonal elements \var{g[k,k]} of \var{g} 
176 
and \code{transpose(g)} returns the transposed of the matrix \var{g} (ie. 
177 
$\var{transpose(g)[i,j]}=\var{g[j,i]}$). 
178 

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

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

203 
The statement \code{myPDE.setLumpingOn()} 
204 
switches on the usage of an aggressive approximation of the PDE operator, in this case 
205 
formed by the coefficient \var{D} (actually the discrete 
206 
version of the PDE operator is approximated by a diagonal matrix). The approximation allows, at the costs of 
207 
an additional error, very fast 
208 
solving of the PDE when the solution is requested. In the general case, the presence of \var{A}, \var{B} or \var 
209 
{C} \code{setLumpingOn} will produce wrong results. 
210 

211 

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

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

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

261 
mydomain=Brick(ne,ne,10,l0=width,l1=width,l2=10.*width/32.) 
262 
wavePropagation(mydomain,h,tend,lam,mu,rho,U0) 
263 
\end{python} 
264 
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 
265 
working directory. 
266 
To visualize the results from the data directory: 
267 
\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. 
268 
Again use \code{Configure Data} to move backwards and forward in time, and 
269 
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. 
270 

271 
\begin{figure}{t!} 
272 
\centerline{\includegraphics[width=3.in,angle=270]{figures/WavePC.eps}} 
273 
\caption{Amplitude at Point source} 
274 
\label{WAVE FIG 1} 
275 
\end{figure} 
276 

277 
\begin{figure}{t} 
278 
\begin{center} 
279 
\includegraphics[width=2in,angle=270]{figures/Wavet0.eps} 
280 
\includegraphics[width=2in,angle=270]{figures/Wavet1.eps} 
281 
\includegraphics[width=2in,angle=270]{figures/Wavet3.eps} 
282 
\includegraphics[width=2in,angle=270]{figures/Wavet10.eps} 
283 
\includegraphics[width=2in,angle=270]{figures/Wavet30.eps} 
284 
\includegraphics[width=2in,angle=270]{figures/Wavet288.eps} 
285 
\end{center} 
286 
\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 
287 
from a point source initially at $(5\mbox{km},5\mbox{km},0)$ 
288 
with time step size $h=0.02083$. Color represents the displacement. 
289 
Here the view is oriented onto the bottom face.} 
290 
\label{WAVE FIG 2} 
291 
\end{figure} 