1 
% $Id$ 
2 
\section{3D Wave Propagation} 
3 
\label{WAVE CHAP} 
4 

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

30 
The points on the left face of the domain, ie. points $x=(x\hackscore{i})$ with $x\hackscore{0}=0$, 
31 
are moved into the $x\hackscore{2}$ direction by $u^{max}$. This is expressed by the constraint 
32 
\begin{eqnarray} \label{WAVE impact} 
33 
u\hackscore{2}(x,t)=s(t) \mbox{ for } x=(x\hackscore{i}) with x\hackscore{0}=0 
34 
\end{eqnarray} 
35 
where 
36 
\begin{eqnarray} \label{WAVE impact s} 
37 
s(t)=u^{max} (1e^{(\tau^{1}t)^3}) 
38 
\end{eqnarray} 
39 
$\tau$ measures the time needed to reach the maximum displacement $u^{max}$ ( 
40 
actually $\tau (ln 2)^{\frac{1}{3}}) \approx 0.9 \tau$ is the time until $\frac{u^{max}}{2}$ is reached). 
41 
The smaller $\tau$ the faster $u^{max}$ is reached. 
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^{(n)} + h^2 a^{(n)} \\ 
48 
\end{eqnarray} 
49 
for all $n=1,2,\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} 
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 
Additionally $a^{(n)}$ has to fullfill the constraint 
68 
\begin{eqnarray}\label{WAVE dyn 4} 
69 
a^{(n)}\hackscore{2}(x)= s\hackscore{,tt}(t^{(n)})\mbox{ where } s\hackscore{,tt}(t)=\frac{u^{max}}{\tau^2} 
70 
(6\tau^{1}t 9(\tau^{1}t)^4) e^{(\tau^{1}t)^3} 
71 
\end{eqnarray} 
72 
derived from \eqn{WAVE impact}. 
73 

74 

75 
In each time step we have to solve this problem to get the acceleration $a^{(n)}$ and we will 
76 
use the \LinearPDE class of the \linearPDEs to do so. The general form of the PDE defined through 
77 
the \LinearPDE class is discussed in \Sec{SEC LinearPDE}. The form which are relevant here are 
78 
\begin{equation}\label{WAVE dyn 100} 
79 
D\hackscore{ij} u\hackscore{j} =  X\hackscore{ij,j}\; . 
80 
\end{equation} 
81 
Implicitly, the natural boundary condition 
82 
\begin{equation}\label{WAVE dyn 101} 
83 
n\hackscore{j}X\hackscore{ij} =0 
84 
\end{equation} 
85 
is assumed. The \LinearPDE object allows defining constraints of the form 
86 
\begin{equation}\label{WAVE dyn 101} 
87 
u\hackscore{i}(x)=r\hackscore{i}(x) \mbox{ for } x \mbox{ with } q\hackscore{i}>0 
88 
\end{equation} 
89 
where $r$ and $q$ are given functions. 
90 

91 
With $u=a^{(n)}$ we can identify the values to be assigned to $D$, $X$, $r$ and $q$: 
92 
\begin{equation}\label{WAVE dyn 6} 
93 
\begin{array}{ll} 
94 
D\hackscore{ij}=\rho \delta\hackscore{ij}& 
95 
X\hackscore{ij}=\sigma^{(n1)}\hackscore{ij} \\ 
96 
r\hackscore{i}=\delta\hackscore{i2} \cdot s\hackscore{,tt}(t^{(n)}) & q\hackscore{i}=\delta\hackscore{i2} \cdot x\hackscore{0}.\texttt{whereZero()} 
97 
\end{array} 
98 
\end{equation} 
99 
Remember that \var{y.whereZero()} returns a function which is $1$ for those points $x$ where $y$ is zero 
100 
(up to a certain tolerance, see \Sec{SEC ESCRIPT DATA}) and $0$ for those points where $y$ is not equal to zero. 
101 

102 
In the following script defines a the function \function{wavePropagation} which 
103 
implements the Verlet scheme to solve our wave propagation problem. 
104 
The argument \var{domain} which is a \Domain class object 
105 
defines the domain of the problem. \var{h} and \var{tend} is the step size 
106 
and the time mark after which the simulation will be terminated. \var{lam}, \var{mu} and 
107 
\var{rho} are material properties. \var{s_tt} is a function which returns $s\hackscore{,tt}$ at a given time: 
108 
\begin{python} 
109 
from esys.linearPDEs import LinearPDE 
110 
from numarray import identity 
111 
from esys.escript import * 
112 
def wavePropagation(domain,h,tend,lam,mu,rho,s_tt): 
113 
x=domain.getX() 
114 
# ... open new PDE ... 
115 
mypde=LinearPDE(domain) 
116 
mypde.setLumpingOn() 
117 
kronecker=identity(mypde.getDim()) 
118 
mypde.setValue(D=kronecker*rho, \ 
119 
q=x[0].whereZero()*kronecker[1,:]) 
120 
# ... set initial values .... 
121 
n=0 
122 
u=Vector(0,ContinuousFunction(domain)) 
123 
u_last=Vector(0,ContinuousFunction(domain)) 
124 
t=0 
125 
while t<tend: 
126 
# ... get current stress .... 
127 
g=grad(u) 
128 
stress=lam*trace(g)*kronecker+mu*(g+transpose(g)) 
129 
# ... get new acceleration .... 
130 
mypde.setValue(X=stress,r=s_tt(t+h)*kronecker[1,:]) 
131 
a=mypde.getSolution() 
132 
# ... get new displacement ... 
133 
u_new=2*uu_last+h**2*a 
134 
# ... shift displacements .... 
135 
u_last=u 
136 
u=u_new 
137 
t+=h 
138 
n+=1 
139 
print n,"th time step t ",t 
140 
print "a=",inf(a),sup(a) 
141 
print "u=",inf(u),sup(u) 
142 
# ... save current acceleration in units of gravity 
143 
if n%10==0 : (length(a)/9.81).saveDX("/tmp/res/u.%i.dx"%(n/10)) 
144 
\end{python} 
145 
Notice that 
146 
all coefficients of the PDE which are independent from time $t$ are set outside the \code{while} 
147 
loop. This allows the \LinearPDE class to ruse information as much as possible 
148 
when iterating over time. 
149 

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

158 
The initialization of \var{u} and \var{u_last} needs some explanation. The statement 
159 
\begin{python} 
160 
u=Vector(0.,ContinuousFunction(domain)) 
161 
\end{python} 
162 
declares \var{u} as a vectorvalued function of its coordinates 
163 
in the domain (i.e. with two or three components in the case of 
164 
a two or three dimensional domain, repectively). The first argument defines the value of the function which is equal 
165 
to $0$ for all vectot components. The second argument \var{ContinuousFunction(domain)} 
166 
specifies the `smoothness` of the function. In our case we want the displacement field to be 
167 
a continuous function over the \Domain \var{domain}. On other option would be 
168 
\var{Function(domain)} in which case continuity is not garanteed. In fact the 
169 
argument of the \function{grad} function has to be of the typ \var{ContinuousFunction(domain)} while 
170 
returned function is of the var{Function(domain)} type. For more details see \Sec{SEC ESCRIPT DATA}. 
171 

172 
The statement \code{myPDE.setLumpingOn()} 
173 
switches on the usage of an aggressive approximation of the PDE operator, in this case 
174 
formed by the coefficient \var{D} (actually the discrete 
175 
version of the PDE operator is approximated by a diagonal matrix). The approximation allows, at the costs of 
176 
an additional error, very fast 
177 
solving of the PDE when the solution is requested. In the general case of the presents of \var{A}, \var{B} or \var 
178 
{C} \code{setLumpingOn} will produce wrong results. 
179 

180 
It is pointed out that the value of the constraint \var{r} for the acceleration \var{a} at time $t=0$ is consistent 
181 
with the initial value of the acceleration which is identical zero. Otherwise, a discontinuity of \var{a} in time 
182 
direction would lead to heavy oscillations in the displacement. It is also pointed out, that at time $t=0$ the 
183 
constraint defined by \eqn{WAVE impact} fulfills 
184 
the initial conditions in \eqn{WAVE initial}. 
185 

186 

187 
One of the big advantages of the Verlet scheme is the fact that the problem to be solved 
188 
in each time step is very simple and does not involve any spatial derivatives (which actually allows us to use lumping). 
189 
The problem becomes so simple because we use the stress from the last time step rather then the stress which 
190 
actually is present at the current time step. Schemes using this approach are called an explicit time integration 
191 
scheme \index{explicit scheme} \index{time integration!explicit}. The 
192 
backward Euler scheme we have used in \Chap{DIFFUSION CHAP} is 
193 
an example of an implicit schemes 
194 
\index{implicit scheme} \index{time integration!implicit}. In this case one uses the actual status of 
195 
each variable at a particular time rather then values from previous time steps. This will lead to problem which is 
196 
more expensive to solve, in particular for nonlinear problems. 
197 
Although 
198 
explicit time integration schemes are cheep to finalize a single time step, they need a significant smaller time 
199 
step then implicit schemes and can suffer from stability problems. Therefor they need a 
200 
very careful selection of the time step size $h$. 
201 

202 
An easy, heuristic way of choosing an appropriate 
203 
time step size is the Courant condition \index{Courant condition} \index{explicit scheme!Courant condition} 
204 
which says that within a time step a information should not travel further than a cell used in the 
205 
discretization scheme. In the case of the wave equation the velocity of a (p) wave is given as 
206 
$\sqrt{\frac{\lambda+2\mu}{\rho}}$ so one should choose $h$ from 
207 
\begin{eqnarray}\label{WAVE dyn 66} 
208 
h= \frac{1}{5} \sqrt{\frac{\rho}{\lambda+2\mu}} \Delta x 
209 
\end{eqnarray} 
210 
where $\Delta x$ is the cell diameter. The factor $\frac{1}{5}$ is a safety factor considering the heuristic of 
211 
the formula. 
212 

213 
The following script uses the \code{wavePropagtion} function to solve the 
214 
wave equation over a block in crust of the earth. The depth of the block is 
215 
$10km$ and the width in each direction is $100km$. The depth of the block is alligned 
216 
with the $x\hackscore{0}$direction. \var{ne} gives the number of elements in $x\hackscore{0}$direction. The number 
217 
of elements in $x\hackscore{1}$ and $x\hackscore{2}$direction is chosen such that the elements 
218 
become cubic. The impact is $2m$ acting within a time of about $10sec$: 
219 
\begin{python} 
220 
ne=10 # number of cells in x_0direction 
221 
depth=10000. # length in x_0direction 
222 
width=100000. # length in x_1 and x_2 direction 
223 
lam=3.462e9 
224 
mu=3.462e9 
225 
rho=1154. 
226 
tau=10. 
227 
umax=2. 
228 
tend=60 
229 
h=1./5.*sqrt(rho/(lam+2*mu))*(depth/ne) 
230 
def s_tt(t): return umax/tau**2*(6*t/tau9*(t/tau)**4)*exp((t/tau)**3) 
231 

232 
print "time step size = ",h 
233 
mydomain=Brick(ne,int(width/depth)*ne,int(width/depth)*ne,l0=depth,l1=width,l2=width) 
234 
wavePropagation(mydomain,h,tend,lam,mu,rho,s_tt) 
235 
\end{python} 
236 
The script is available as \file{wave.py} in the \ExampleDirectory \index{scripts!\file{wave.py}}. 
237 

238 

239 
% \begin{figure} 
240 
% \centerline{\includegraphics[width=\figwidth]{DiffusionDomain}} 
241 
% \caption{Temperature Diffusion Problem with Circular Heat Source} 
242 
% \label{WAVE FIG 1} 
243 
% \end{figure} 
244 

245 
\begin{figure} 
246 
% \centerline{\includegraphics[width=\figwidth]{DiffusionRes1}} 
247 
% \centerline{\includegraphics[width=\figwidth]{DiffusionRes16}} 
248 
% \centerline{\includegraphics[width=\figwidth]{DiffusionRes32}} 
249 
%\centerline{\includegraphics[width=\figwidth]{DiffusionRes48}} 
250 
\caption{Selected time steps of a wave propagation over a $10km \times 100km \times 100km$ 
251 
with time step size $h=0.0666$. Color represents the acceleration.} 
252 
\label{WAVE FIG 2} 
253 
\end{figure} 