23 
\nabla ^2 p  \frac{1}{c^2} \frac{\partial ^2 p}{\partial t^2} = 0 
\nabla ^2 p  \frac{1}{c^2} \frac{\partial ^2 p}{\partial t^2} = 0 
24 
\label{eqn:acswave} 
\label{eqn:acswave} 
25 
\end{equation} 
\end{equation} 
26 
where $p$ is the pressure, $t$ is the time and $c$ is the wave velocity. 
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 

equation stability criterion and solving for the displacement of acceleration solution. 
30 


31 
\section{The Laplacian in \esc} 
\section{The Laplacian in \esc} 
32 
The Laplacian opperator which can be written as $\Delta$ or $\nabla^2$ is 
The Laplacian operator which can be written as $\Delta$ or $\nabla^2$ is 
33 
calculated via the divergence of the gradient of the object, which is in this 
calculated via the divergence of the gradient of the object, which in this 
34 
example $p$. Thus we can write; 
example is the scalar $p$. Thus we can write; 
35 
\begin{equation} 
\begin{equation} 
36 
\nabla^2 p = \nabla \cdot \nabla p = 
\nabla^2 p = \nabla \cdot \nabla p = 
37 
\sum_{i}^n 
\sum_{i}^n 
46 
\end{equation} 
\end{equation} 
47 


48 
In \esc the Laplacian is calculated using the divergence representation and the 
In \esc the Laplacian is calculated using the divergence representation and the 
49 
intrinsic functions \textit{grad()} and \textit{trace()}. The fucntion 
intrinsic functions \textit{grad()} and \textit{trace()}. The function 
50 
\textit{grad{}} will return the spatial gradients of an object. 
\textit{grad{}} will return the spatial gradients of an object. 
51 
For a rank 0 solution, this is of the form; 
For a rank 0 solution, this is of the form; 
52 
\begin{equation} 
\begin{equation} 
73 
$X$. Notice however that the general form contains the term $X 
$X$. Notice however that the general form contains the term $X 
74 
_{i,j}$\footnote{This is the first derivative in the $j^{th}$ 
_{i,j}$\footnote{This is the first derivative in the $j^{th}$ 
75 
direction for the $i^{th}$ component of the solution.}, 
direction for the $i^{th}$ component of the solution.}, 
76 
hence for a rank 0 object there is no need to do more than calculate the 
hence for a rank 0 object there is no need to do more then calculate the 
77 
gradient and submit it to the solver. In the case of the rank 1 or greater 
gradient and submit it to the solver. In the case of the rank 1 or greater 
78 
object, it is nesscary to calculate the trace also. This is the sum of the 
object, it is necessary to calculate the trace also. This is the sum of the 
79 
diagonal in \autoref{eqn:gradrank1}. 
diagonal in \autoref{eqn:gradrank1}. 
80 


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


103 
\section{Numerical Solution Stability} \label{sec:nsstab} 
\section{Numerical Solution Stability} \label{sec:nsstab} 
104 
Unfortunately, the wave equation belongs to a class of equations called 
Unfortunately, the wave equation belongs to a class of equations called 
105 
\textbf{stiff} PDEs. These types of equations can be difficult to solve 
\textbf{stiff} PDEs. These types of equations can be difficult to solve 
106 
numerically as they tend to oscilate about the exact solution which can 
numerically as they tend to oscillate about the exact solution which can 
107 
eventually lead to a catastrophic failure in the solution. To counter this 
eventually lead to a catastrophic failure. To counter this problem, explicitly 
108 
problem, explicitly stable schemes like 
stable schemes like the backwards Euler method and correct parameterisation of 
109 
the backwards Euler method are required. There are two variables which must be 
the problem are required. 
110 
considered for stability when numerically trying to solve the wave equation. 

111 
For linear media, the two variables are related via; 
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 
\begin{equation} \label{eqn:freqvel} 
\begin{equation} \label{eqn:freqvel} 
115 
f=\frac{v}{\lambda} 
f=\frac{v}{\lambda} 
116 
\end{equation} 
\end{equation} 
117 
The velocity $v$ that a wave travels in a medium is an important variable. For 
The velocity $v$ that a wave travels in a medium is an important variable. For 
118 
stability the analytical wave must not propagate faster than the numerical wave 
stability the analytical wave must not propagate faster then the numerical wave 
119 
is able to, and in general, needs to be much slower than the numerical wave. 
is able to, and in general, needs to be much slower then the numerical wave. 
120 
For example, a line 100m long is discretised into 1m intervals or 101 nodes. If 
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 
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 
the wave between each node will be 0.01 seconds. The time step, must therefore 
123 
be significantly less than this. Of the order $10E4$ would be appropriate. 
be significantly less then this. Of the order $10E4$ would be appropriate. 
124 


125 
The wave frequency content also plays a part in numerical stability. The 
The wave frequency content also plays a part in numerical stability. The 
126 
nyquistsampling theorem states that a signals bandwidth content will be 
nyquistsampling theorem states that a signals bandwidth content will be 
127 
accurately represented when an equispaced sampling rate $f _{n}$ is 
accurately represented when an equispaced sampling rate $f _{n}$ is 
128 
equal to or greater than twice the maximum frequency of the signal 
equal to or greater then twice the maximum frequency of the signal 
129 
$f_{s}$, or; 
$f_{s}$, or; 
130 
\begin{equation} \label{eqn:samptheorem} 
\begin{equation} \label{eqn:samptheorem} 
131 
f_{n} \geqslant f_{s} 
f_{n} \geqslant f_{s} 
132 
\end{equation} 
\end{equation} 
133 
For example a 50Hz signal will require a sampling rate greater than 100Hz or 
For example a 50Hz signal will require a sampling rate greater then 100Hz or 
134 
one sample every 0.01 seconds. The wave equation relies on a spatial frequency, 
one sample every 0.01 seconds. The wave equation relies on a spatial frequency, 
135 
thus the sampling theorem in this case applies to the solution mesh spacing. In 
thus the sampling theorem in this case applies to the solution mesh spacing. In 
136 
this way, the frequency content of the input signal directly affects the time 
this way, the frequency content of the input signal directly affects the time 
178 
\item The time iteration and control parameters need to be defined. 
\item The time iteration and control parameters need to be defined. 
179 
\item The PDE is initialised with source and boundary conditions. 
\item The PDE is initialised with source and boundary conditions. 
180 
\item The time loop is started and the PDE is solved at consecutive time steps. 
\item The time loop is started and the PDE is solved at consecutive time steps. 
181 
\item All or select solutions are saved to file for visualisation lated on. 
\item All or select solutions are saved to file for visualisation later on. 
182 
\end{enumerate} 
\end{enumerate} 
183 


184 
Parts of the script which warrant more attention are the definition of the 
Parts of the script which warrant more attention are the definition of the 
191 
respect to time. This is often a good opportunity to introduce a source to the 
respect to time. This is often a good opportunity to introduce a source to the 
192 
solution. This model has the source located at it's centre. The source should 
solution. This model has the source located at it's centre. The source should 
193 
be smooth and cover a number of samples to satisfy the frequency stability 
be smooth and cover a number of samples to satisfy the frequency stability 
194 
criterion. Small sources will generate high frequency signals. Here, the source 
criterion. Small sources will generate high frequency signals. Here, when using 
195 
is defined by a cosine function. 
a rectangular domain, the source is defined by a cosine function. 
196 
\begin{python} 
\begin{python} 
197 
U0=0.01 # amplitude of point source 
U0=0.01 # amplitude of point source 
198 
xc=[500,500] #location of point source 
xc=[500,500] #location of point source 
203 
whereNegative(length(xxc)src_radius) 
whereNegative(length(xxc)src_radius) 
204 
u_m1=u 
u_m1=u 
205 
\end{python} 
\end{python} 

When using a rectangular domain 

206 


207 
\subsection{Visualising the Source} 
\subsection{Visualising the Source} 
208 
There are two options for visualising the source. The first is to export the 
There are two options for visualising the source. The first is to export the 
209 
initial conditions of the model to VTK, which can be interpreted as a scalar 
initial conditions of the model to VTK, which can be interpreted as a scalar 
210 
suface in mayavi. The second is to take a cross section of the model. 
surface in Mayavi2. The second is to take a cross section of the model which 
211 

will require the \textit{Locator} function. 

For the later, we will require the \textit{Locator} function. 

212 
First \verb!Locator! must be imported; 
First \verb!Locator! must be imported; 
213 
\begin{python} 
\begin{python} 
214 
from esys.escript.pdetools import Locator 
from esys.escript.pdetools import Locator 
217 
to the point or points of interest. 
to the point or points of interest. 
218 


219 
It is now necessary to build a list of $(x,y)$ locations that specify where are 
It is now necessary to build a list of $(x,y)$ locations that specify where are 
220 
model slice will go. This is easily implemeted with a loop; 
model slice will go. This is easily implemented with a loop; 
221 
\begin{python} 
\begin{python} 
222 
cut_loc=[] 
cut_loc=[] 
223 
src_cut=[] 
src_cut=[] 
231 
src=Locator(mydomain,src_cut) 
src=Locator(mydomain,src_cut) 
232 
src_cut=src.getValue(u) 
src_cut=src.getValue(u) 
233 
\end{python} 
\end{python} 
234 
It is then a trivial task to plot and save the output using \mpl. 
It is then a trivial task to plot and save the output using \mpl 
235 

(\autoref{fig:cxsource}). 
236 
\begin{python} 
\begin{python} 
237 
pl.plot(cut_loc,src_cut) 
pl.plot(cut_loc,src_cut) 
238 
pl.axis([xc[0]src_radius*3,xc[0]+src_radius*3,0.,2*U0]) 
pl.axis([xc[0]src_radius*3,xc[0]+src_radius*3,0.,2*U0]) 
307 
\section{Stability Investigation} 
\section{Stability Investigation} 
308 
It is now prudent to investigate the stability limitations of this problem. 
It is now prudent to investigate the stability limitations of this problem. 
309 
First, we let the frequency content of the source be very small. If we define 
First, we let the frequency content of the source be very small. If we define 
310 
the source as a cosine input, than the wavlength of the input is equal to the 
the source as a cosine input, then the wavlength of the input is equal to the 
311 
radius of the source. Let this value be 5 meters. Now, if the maximum velocity 
radius of the source. Let this value be 5 meters. Now, if the maximum velocity 
312 
of the model is $c=380.0ms^{1}$ then the source 
of the model is $c=380.0ms^{1}$ then the source 
313 
frequency is $f_{r} = \frac{380.0}{5} = 76.0 Hz$. This is a worst case 
frequency is $f_{r} = \frac{380.0}{5} = 76.0 Hz$. This is a worst case 