1 
ahallam 
3003 

2 
jfenwick 
3989 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
3 
jfenwick 
6651 
% Copyright (c) 20032018 by The University of Queensland 
4 
jfenwick 
3989 
% http://www.uq.edu.au 
5 
ahallam 
3003 
% 
6 


% Primary Business: Queensland, Australia 
7 
jfenwick 
6112 
% Licensed under the Apache License, version 2.0 
8 


% http://www.apache.org/licenses/LICENSE2.0 
9 
ahallam 
3003 
% 
10 
jfenwick 
3989 
% Development until 2012 by Earth Systems Science Computational Center (ESSCC) 
11 
jfenwick 
4657 
% Development 20122013 by School of Earth Sciences 
12 


% Development from 2014 by Centre for Geoscience Computing (GeoComp) 
13 
jfenwick 
3989 
% 
14 


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
15 
ahallam 
3003 

16 


The acoustic wave equation governs the propagation of pressure waves. Wave 
17 


types that obey this law tend to travel in liquids or gases where shear waves 
18 
ahallam 
3004 
or longitudinal style wave motion is not possible. An obvious example is sound 
19 
ahallam 
3003 
waves. 
20 



21 
ahallam 
3004 
The acoustic wave equation is defined as; 
22 
ahallam 
3003 
\begin{equation} 
23 


\nabla ^2 p  \frac{1}{c^2} \frac{\partial ^2 p}{\partial t^2} = 0 
24 


\label{eqn:acswave} 
25 


\end{equation} 
26 
ahallam 
3370 
where $p$ is the pressure, $t$ is the time and $c$ is the wave velocity. In this 
27 


chapter the acoustic wave equation is demonstrated. Important steps include the 
28 


translation of the Laplacian $\nabla^2$ to the \esc general form, the stiff 
29 
ahallam 
3373 
equation stability criterion and solving for the displacement or acceleration solution. 
30 
ahallam 
3003 

31 
ahallam 
3004 
\section{The Laplacian in \esc} 
32 
ahallam 
3373 
The Laplacian operator which can be written as $\Delta$ or $\nabla^2$, is 
33 
ahallam 
3370 
calculated via the divergence of the gradient of the object, which in this 
34 


example is the scalar $p$. Thus we can write; 
35 
ahallam 
3004 
\begin{equation} 
36 
ahallam 
3029 
\nabla^2 p = \nabla \cdot \nabla p = 
37 
jfenwick 
3308 
\sum_{i}^n 
38 


\frac{\partial^2 p}{\partial x^2_{i}} 
39 
ahallam 
3004 
\label{eqn:laplacian} 
40 


\end{equation} 
41 
ahallam 
3232 
For the two dimensional case in Cartesian coordinates \autoref{eqn:laplacian} 
42 
ahallam 
3004 
becomes; 
43 


\begin{equation} 
44 


\nabla^2 p = \frac{\partial^2 p}{\partial x^2} 
45 


+ \frac{\partial^2 p}{\partial y^2} 
46 


\end{equation} 
47 
ahallam 
3003 

48 
ahallam 
3004 
In \esc the Laplacian is calculated using the divergence representation and the 
49 
ahallam 
3370 
intrinsic functions \textit{grad()} and \textit{trace()}. The function 
50 
ahallam 
3004 
\textit{grad{}} will return the spatial gradients of an object. 
51 


For a rank 0 solution, this is of the form; 
52 


\begin{equation} 
53 


\nabla p = \left[ 
54 
jfenwick 
3308 
\frac{\partial p}{\partial x _{0}}, 
55 


\frac{\partial p}{\partial x _{1}} 
56 
ahallam 
3004 
\right] 
57 


\label{eqn:grad} 
58 


\end{equation} 
59 


Larger ranked solution objects will return gradient tensors. For example, a 
60 
jfenwick 
3308 
pressure field which acts in the directions $p _{0}$ and $p 
61 


_{1}$ would return; 
62 
ahallam 
3004 
\begin{equation} 
63 


\nabla p = \begin{bmatrix} 
64 
jfenwick 
3308 
\frac{\partial p _{0}}{\partial x _{0}} & 
65 


\frac{\partial p _{1}}{\partial x _{0}} \\ 
66 


\frac{\partial p _{0}}{\partial x _{1}} & 
67 


\frac{\partial p _{1}}{\partial x _{1}} 
68 
ahallam 
3004 
\end{bmatrix} 
69 


\label{eqn:gradrank1} 
70 


\end{equation} 
71 
ahallam 
3003 

72 
ahallam 
3232 
\autoref{eqn:grad} corresponds to the Linear PDE general form value 
73 
ahallam 
3373 
$X$. Notice however, that the general form contains the term $X 
74 
jfenwick 
3308 
_{i,j}$\footnote{This is the first derivative in the $j^{th}$ 
75 
ahallam 
3029 
direction for the $i^{th}$ component of the solution.}, 
76 
ahallam 
3370 
hence for a rank 0 object there is no need to do more then calculate the 
77 
ahallam 
3004 
gradient and submit it to the solver. In the case of the rank 1 or greater 
78 
ahallam 
3373 
object, it is also necessary to calculate the trace. This is the sum of the 
79 
ahallam 
3232 
diagonal in \autoref{eqn:gradrank1}. 
80 
ahallam 
3004 

81 


Thus when solving for equations containing the Laplacian one of two things must 
82 
ahallam 
3392 
be completed. If the object \verb!p! is less than rank 1 the gradient is 
83 
ahallam 
3004 
calculated via; 
84 
ahallam 
3025 
\begin{python} 
85 
ahallam 
3063 
gradient=grad(p) 
86 
ahallam 
3025 
\end{python} 
87 
ahallam 
3370 
and if the object is greater then or equal to a rank 1 tensor, the trace of 
88 
ahallam 
3004 
the gradient is calculated. 
89 
ahallam 
3025 
\begin{python} 
90 
ahallam 
3004 
gradient=trace(grad(p)) 
91 
ahallam 
3025 
\end{python} 
92 
ahallam 
3370 
These values can then be submitted to the PDE solver via the general form term 
93 
ahallam 
3004 
$X$. The Laplacian is then computed in the solution process by taking the 
94 


divergence of $X$. 
95 



96 
ahallam 
3025 
Note, if you are unsure about the rank of your tensor, the \textit{getRank} 
97 


command will return the rank of the PDE object. 
98 


\begin{python} 
99 


rank = p.getRank() 
100 


\end{python} 
101 



102 



103 


\section{Numerical Solution Stability} \label{sec:nsstab} 
104 
ahallam 
3004 
Unfortunately, the wave equation belongs to a class of equations called 
105 


\textbf{stiff} PDEs. These types of equations can be difficult to solve 
106 
ahallam 
3373 
numerically as they tend to oscillate about the exact solution, which can 
107 
ahallam 
3370 
eventually lead to a catastrophic failure. To counter this problem, explicitly 
108 
ahallam 
3373 
stable schemes like the backwards Euler method, and correct parameterisation of 
109 
ahallam 
3370 
the problem are required. 
110 



111 


There are two variables which must be considered for 
112 


stability when numerically trying to solve the wave equation. For linear media, 
113 


the two variables are related via; 
114 
ahallam 
3004 
\begin{equation} \label{eqn:freqvel} 
115 


f=\frac{v}{\lambda} 
116 


\end{equation} 
117 
ahallam 
3025 
The velocity $v$ that a wave travels in a medium is an important variable. For 
118 
ahallam 
3370 
stability the analytical wave must not propagate faster then the numerical wave 
119 


is able to, and in general, needs to be much slower then the numerical wave. 
120 
ahallam 
3003 
For example, a line 100m long is discretised into 1m intervals or 101 nodes. If 
121 


a wave enters with a propagation velocity of 100m/s then the travel time for 
122 


the wave between each node will be 0.01 seconds. The time step, must therefore 
123 
caltinay 
4285 
be significantly less than this. Of the order $10E4$ would be appropriate. 
124 
ahallam 
3392 
This stability criterion is known as the Courant\textendash 
125 


Friedrichs\textendash Lewy condition given by 
126 
ahallam 
3384 
\begin{equation} 
127 


dt=f\cdot \frac{dx}{v} 
128 


\end{equation} 
129 


where $dx$ is the mesh size and $f$ is a safety factor. To obtain a time step of 
130 


$10E4$, a safety factor of $f=0.1$ was used. 
131 
ahallam 
3003 

132 
ahallam 
3004 
The wave frequency content also plays a part in numerical stability. The 
133 
ahallam 
3392 
Nyquistsampling theorem states that a signals bandwidth content will be 
134 
jfenwick 
3308 
accurately represented when an equispaced sampling rate $f _{n}$ is 
135 
ahallam 
3370 
equal to or greater then twice the maximum frequency of the signal 
136 
jfenwick 
3308 
$f_{s}$, or; 
137 
ahallam 
3004 
\begin{equation} \label{eqn:samptheorem} 
138 
jfenwick 
3308 
f_{n} \geqslant f_{s} 
139 
ahallam 
3004 
\end{equation} 
140 
ahallam 
3384 
For example, a 50Hz signal will require a sampling rate greater then 100Hz or 
141 
ahallam 
3004 
one sample every 0.01 seconds. The wave equation relies on a spatial frequency, 
142 
ahallam 
3373 
thus the sampling theorem in this case applies to the solution mesh spacing. 
143 


This relationship confirms that the frequency content of the input signal 
144 


directly affects the time discretisation of the problem. 
145 
ahallam 
3004 

146 


To accurately model the wave equation with high resolutions and velocities 
147 


means that very fine spatial and time discretisation is necessary for most 
148 
ahallam 
3384 
problems. This requirement makes the wave equation arduous to 
149 
ahallam 
3003 
solve numerically due to the large number of time iterations required in each 
150 
ahallam 
3004 
solution. Models with very high velocities and frequencies will be the worst 
151 
ahallam 
3029 
affected by this problem. 
152 
ahallam 
3003 

153 


\section{Displacement Solution} 
154 


\sslist{example07a.py} 
155 



156 


We begin the solution to this PDE with the centred difference formula for the 
157 


second derivative; 
158 


\begin{equation} 
159 


f''(x) \approx \frac{f(x+h  2f(x) + f(xh)}{h^2} 
160 


\label{eqn:centdiff} 
161 


\end{equation} 
162 
ahallam 
3232 
substituting \autoref{eqn:centdiff} for $\frac{\partial ^2 p }{\partial t ^2}$ 
163 


in \autoref{eqn:acswave}; 
164 
ahallam 
3003 
\begin{equation} 
165 
jfenwick 
3308 
\nabla ^2 p  \frac{1}{c^2h^2} \left[p_{(t+1)}  2p_{(t)} + 
166 


p_{(t1)} \right] 
167 
ahallam 
3003 
= 0 
168 


\label{eqn:waveu} 
169 


\end{equation} 
170 


Rearranging for $p_{(t+1)}$; 
171 


\begin{equation} 
172 
jfenwick 
3308 
p_{(t+1)} = c^2 h^2 \nabla ^2 p_{(t)} +2p_{(t)}  
173 


p_{(t1)} 
174 
ahallam 
3003 
\end{equation} 
175 


this can be compared with the general form of the \modLPDE module and it 
176 
jfenwick 
3308 
becomes clear that $D=1$, $X_{i,j}=c^2 h^2 \nabla ^2 p_{(t)}$ and 
177 
ahallam 
3004 
$Y=2p_{(t)}  p_{(t1)}$. 
178 
ahallam 
3003 

179 
ahallam 
3025 
The solution script is similar to others that we have created in previous 
180 
ahallam 
3004 
chapters. The general steps are; 
181 


\begin{enumerate} 
182 


\item The necessary libraries must be imported. 
183 


\item The domain needs to be defined. 
184 


\item The time iteration and control parameters need to be defined. 
185 


\item The PDE is initialised with source and boundary conditions. 
186 


\item The time loop is started and the PDE is solved at consecutive time steps. 
187 
ahallam 
3370 
\item All or select solutions are saved to file for visualisation later on. 
188 
ahallam 
3004 
\end{enumerate} 
189 



190 


Parts of the script which warrant more attention are the definition of the 
191 


source, visualising the source, the solution time loop and the VTK data export. 
192 



193 


\subsection{Pressure Sources} 
194 


As the pressure is a scalar, one need only define the pressure for two 
195 


time steps prior to the start of the solution loop. Two known solutions are 
196 


required because the wave equation contains a double partial derivative with 
197 


respect to time. This is often a good opportunity to introduce a source to the 
198 


solution. This model has the source located at it's centre. The source should 
199 


be smooth and cover a number of samples to satisfy the frequency stability 
200 
ahallam 
3370 
criterion. Small sources will generate high frequency signals. Here, when using 
201 


a rectangular domain, the source is defined by a cosine function. 
202 
ahallam 
3025 
\begin{python} 
203 
ahallam 
3004 
U0=0.01 # amplitude of point source 
204 


xc=[500,500] #location of point source 
205 


# define small radius around point xc 
206 


src_radius = 30 
207 


# for first two time steps 
208 
ahallam 
3025 
u=U0*(cos(length(xxc)*3.1415/src_radius)+1)*\ 
209 


whereNegative(length(xxc)src_radius) 
210 
ahallam 
3004 
u_m1=u 
211 
ahallam 
3025 
\end{python} 
212 
ahallam 
3004 

213 
ahallam 
3025 
\subsection{Visualising the Source} 
214 


There are two options for visualising the source. The first is to export the 
215 


initial conditions of the model to VTK, which can be interpreted as a scalar 
216 
caltinay 
5662 
surface in \mayavi. The second is to take a cross section of the model which 
217 
ahallam 
3370 
will require the \textit{Locator} function. 
218 
ahallam 
3063 
First \verb!Locator! must be imported; 
219 
ahallam 
3025 
\begin{python} 
220 


from esys.escript.pdetools import Locator 
221 


\end{python} 
222 


The function can then be used on the domain to locate the nearest domain node 
223 


to the point or points of interest. 
224 



225 


It is now necessary to build a list of $(x,y)$ locations that specify where are 
226 
ahallam 
3370 
model slice will go. This is easily implemented with a loop; 
227 
ahallam 
3025 
\begin{python} 
228 


cut_loc=[] 
229 


src_cut=[] 
230 


for i in range(ndx/2ndx/10,ndx/2+ndx/10): 
231 


cut_loc.append(xstep*i) 
232 


src_cut.append([xstep*i,xc[1]]) 
233 


\end{python} 
234 
ahallam 
3063 
We then submit the output to \verb!Locator! and finally return the appropriate 
235 


values using the \verb!getValue! function. 
236 
ahallam 
3025 
\begin{python} 
237 
ahallam 
3373 
src=Locator(mydomain,src_cut) 
238 
ahallam 
3025 
src_cut=src.getValue(u) 
239 


\end{python} 
240 
ahallam 
3370 
It is then a trivial task to plot and save the output using \mpl 
241 


(\autoref{fig:cxsource}). 
242 
ahallam 
3025 
\begin{python} 
243 
ahallam 
3373 
pl.plot(cut_loc,src_cut) 
244 
ahallam 
3025 
pl.axis([xc[0]src_radius*3,xc[0]+src_radius*3,0.,2*U0]) 
245 


pl.savefig(os.path.join(savepath,"source_line.png")) 
246 


\end{python} 
247 


\begin{figure}[h] 
248 
ahallam 
3029 
\centering 
249 
ahallam 
3063 
\includegraphics[width=6in]{figures/sourceline.png} 
250 
ahallam 
3025 
\caption{Cross section of the source function.} 
251 


\label{fig:cxsource} 
252 


\end{figure} 
253 



254 



255 


\subsection{Point Monitoring} 
256 


In the more general case where the solution mesh is irregular or specific 
257 


locations need to be monitored, it is simple enough to use the \textit{Locator} 
258 


function. 
259 


\begin{python} 
260 


rec=Locator(mydomain,[250.,250.]) 
261 


\end{python} 
262 
ahallam 
3029 
When the solution \verb u is updated we can extract the value at that point 
263 
ahallam 
3025 
via; 
264 


\begin{python} 
265 


u_rec=rec.getValue(u) 
266 


\end{python} 
267 
ahallam 
3063 
For consecutive time steps one can record the values from \verb!u_rec! in an 
268 


array initialised as \verb!u_rec0=[]! with; 
269 
ahallam 
3025 
\begin{python} 
270 


u_rec0.append(rec.getValue(u)) 
271 


\end{python} 
272 



273 


It can be useful to monitor the value at a single or multiple individual points 
274 


in the model during the modelling process. This is done using 
275 
ahallam 
3063 
the \verb!Locator! function. 
276 
ahallam 
3025 

277 



278 
ahallam 
3003 
\section{Acceleration Solution} 
279 


\sslist{example07b.py} 
280 



281 
ahallam 
3373 
An alternative method to the displacement solution, is to solve for the 
282 


acceleration $\frac{\partial ^2 p}{\partial t^2}$ directly. The displacement can 
283 


then be derived from the acceleration after a solution has been calculated 
284 


The acceleration is given by a modified form of \autoref{eqn:waveu}; 
285 
ahallam 
3003 
\begin{equation} 
286 


\nabla ^2 p  \frac{1}{c^2} a = 0 
287 


\label{eqn:wavea} 
288 


\end{equation} 
289 
jfenwick 
3308 
and can be solved directly with $Y=0$ and $X=c^2 \nabla ^2 p_{(t)}$. 
290 
ahallam 
3003 
After each iteration the displacement is reevaluated via; 
291 


\begin{equation} 
292 
jfenwick 
3308 
p_{(t+1)}=2p_{(t)}  p_{(t1)} + h^2a 
293 
ahallam 
3003 
\end{equation} 
294 



295 
ahallam 
3029 
\subsection{Lumping} 
296 
caltinay 
4286 
For \esc, the acceleration solution is preferred as it allows the use of matrix 
297 
ahallam 
3004 
lumping. Lumping or mass lumping as it is sometimes known, is the process of 
298 


aggressively approximating the density elements of a mass matrix into the main 
299 
caltinay 
4286 
diagonal. The use of Lumping is motivated by the simplicity of diagonal matrix 
300 
ahallam 
3004 
inversion. As a result, Lumping can significantly reduce the computational 
301 
ahallam 
3029 
requirements of a problem. Care should be taken however, as this 
302 


function can only be used when the $A$, $B$ and $C$ coefficients of the 
303 


general form are zero. 
304 
ahallam 
3003 

305 
ahallam 
3384 
More information about the lumping implementation used in \esc and its accuracy 
306 


can be found in the user guide. 
307 



308 
ahallam 
3004 
To turn lumping on in \esc one can use the command; 
309 
ahallam 
3025 
\begin{python} 
310 
sshaw 
4821 
mypde.getSolverOptions().setSolverMethod(SolverOptions.HRZ_LUMPING) 
311 
ahallam 
3025 
\end{python} 
312 
ahallam 
3004 
It is also possible to check if lumping is set using; 
313 
ahallam 
3025 
\begin{python} 
314 
jfenwick 
4853 
print(mypde.isUsingLumping()) 
315 
ahallam 
3025 
\end{python} 
316 
ahallam 
3004 

317 


\section{Stability Investigation} 
318 


It is now prudent to investigate the stability limitations of this problem. 
319 
ahallam 
3025 
First, we let the frequency content of the source be very small. If we define 
320 
ahallam 
3370 
the source as a cosine input, then the wavlength of the input is equal to the 
321 
ahallam 
3025 
radius of the source. Let this value be 5 meters. Now, if the maximum velocity 
322 
ahallam 
3373 
of the model is $c=380.0ms^{1}$, then the source 
323 
jfenwick 
3308 
frequency is $f_{r} = \frac{380.0}{5} = 76.0 Hz$. This is a worst case 
324 
ahallam 
3025 
scenario with a small source and the models maximum velocity. 
325 



326 
ahallam 
3232 
Furthermore, we know from \autoref{sec:nsstab}, that the spatial sampling 
327 
ahallam 
3025 
frequency must be at least twice this value to ensure stability. If we assume 
328 


the model mesh is a square equispaced grid, 
329 


then the sampling interval is the side length divided by the number of samples, 
330 


given by $\Delta x = \frac{1000.0m}{400} = 2.5m$ and the maximum sampling 
331 


frequency capable at this interval is 
332 
jfenwick 
3308 
$f_{s}=\frac{380.0ms^{1}}{2.5m}=152Hz$ this is just equal to the 
333 
ahallam 
3232 
required rate satisfying \autoref{eqn:samptheorem}. 
334 
ahallam 
3004 

335 
ahallam 
3232 
\autoref{fig:ex07sampth} depicts three examples where the grid has been 
336 
ahallam 
3025 
undersampled, sampled correctly, and over sampled. The grids used had 
337 


200, 400 and 800 nodes per side respectively. Obviously, the oversampled grid 
338 


retains the best resolution of the modelled wave. 
339 
ahallam 
3004 

340 
ahallam 
3025 
The time step required for each of these examples is simply calculated from 
341 


the propagation requirement. For a maximum velocity of $380.0ms^{1}$, 
342 


\begin{subequations} 
343 


\begin{equation} 
344 


\Delta t \leq \frac{1000.0m}{200} \frac{1}{380.0} = 0.013s 
345 


\end{equation} 
346 


\begin{equation} 
347 


\Delta t \leq \frac{1000.0m}{400} \frac{1}{380.0} = 0.0065s 
348 


\end{equation} 
349 


\begin{equation} 
350 


\Delta t \leq \frac{1000.0m}{800} \frac{1}{380.0} = 0.0032s 
351 


\end{equation} 
352 


\end{subequations} 
353 
ahallam 
3373 
Observe that for each doubling of the number of nodes in the mesh, we halve 
354 


the time step. To illustrate the impact this has, consider our model. If the 
355 
ahallam 
3025 
source is placed at the center, it is $500m$ from the nearest boundary. With a 
356 


velocity of $380.0ms^{1}$ it will take $\approx1.3s$ for the wavefront to 
357 


reach that boundary. In each case, this equates to $100$, $200$ and $400$ time 
358 


steps. This is again, only a best case scenario, for true stability these time 
359 
ahallam 
3373 
values may need to be halved and possibly halved again. 
360 
ahallam 
3004 

361 
ahallam 
3025 
\begin{figure}[ht] 
362 


\centering 
363 


\subfigure[Undersampled Example]{ 
364 
caltinay 
3386 
\includegraphics[width=0.45\textwidth,trim=0cm 6cm 5cm 6cm,clip]{figures/ex07usamp.png} 
365 
ahallam 
3025 
\label{fig:ex07usamp} 
366 


} 
367 


\subfigure[Just sampled Example]{ 
368 
caltinay 
3386 
\includegraphics[width=0.45\textwidth,trim=0cm 6cm 5cm 6cm,clip]{figures/ex07jsamp.png} 
369 
ahallam 
3025 
\label{fig:ex07jsamp} 
370 


} 
371 


\subfigure[Over sampled Example]{ 
372 
caltinay 
3386 
\includegraphics[width=0.45\textwidth,trim=0cm 6cm 5cm 6cm,clip]{figures/ex07nsamp.png} 
373 
ahallam 
3025 
\label{fig:ex07nsamp} 
374 


} 
375 
caltinay 
3386 
\caption{Sampling Theorem example for stability investigation} 
376 
ahallam 
3373 
\label{fig:ex07sampth} 
377 
ahallam 
3025 
\end{figure} 
378 
ahallam 
3004 
