1 
jgs 
107 
% $Id$ 
2 
jgs 
121 
\section{3D Wave Propagation} 
3 
jgs 
107 
\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 
jgs 
110 

75 
jgs 
107 
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 \linearPDEsPack 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 
jgs 
110 
D\hackscore{ij} u\hackscore{j} =  X\hackscore{ij,j}\; . 
80 
jgs 
107 
\end{equation} 
81 
jgs 
110 
Implicitly, the natural boundary condition 
82 
jgs 
107 
\begin{equation}\label{WAVE dyn 101} 
83 


n\hackscore{j}X\hackscore{ij} =0 
84 


\end{equation} 
85 
jgs 
110 
is assumed. The \LinearPDE object allows defining constraints of the form 
86 
jgs 
107 
\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 
jgs 
110 
\var{rho} are material properties. \var{s_tt} is a function which returns $s\hackscore{,tt}$ at a given time: 
108 
jgs 
107 
\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 
jgs 
110 
kronecker=identity(mypde.getDim()) 
118 
jgs 
107 
mypde.setValue(D=kronecker*rho, \ 
119 
jgs 
110 
q=x[0].whereZero()*kronecker[1,:]) 
120 
jgs 
107 
# ... set initial values .... 
121 
jgs 
110 
n=0 
122 
jgs 
107 
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 
jgs 
110 
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 
jgs 
107 
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 
jgs 
110 
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 
jgs 
107 
\end{python} 
145 
jgs 
110 
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 
jgs 
107 
\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 
jgs 
110 
u=Vector(0.,ContinuousFunction(domain)) 
161 
jgs 
107 
\end{python} 
162 
jgs 
110 
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 
jgs 
107 

172 
jgs 
110 
The statement \code{myPDE.setLumpingOn()} 
173 
jgs 
107 
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 
jgs 
110 

180 
jgs 
107 
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 
jgs 
110 
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 
jgs 
107 

202 
jgs 
110 
An easy, heuristic way of choosing an appropriate 
203 
jgs 
107 
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 
jgs 
110 
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 
jgs 
107 
\begin{python} 
220 
jgs 
110 
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 
jgs 
107 
rho=1154. 
226 
jgs 
110 
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 
jgs 
107 
\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 
jgs 
110 
\begin{figure} 
246 
jgs 
107 
% \centerline{\includegraphics[width=\figwidth]{DiffusionRes1}} 
247 


% \centerline{\includegraphics[width=\figwidth]{DiffusionRes16}} 
248 


% \centerline{\includegraphics[width=\figwidth]{DiffusionRes32}} 
249 


%\centerline{\includegraphics[width=\figwidth]{DiffusionRes48}} 
250 
jgs 
110 
\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} 