35 
In a similar process to the previous chapter, we will use the acceleration 
In a similar process to the previous chapter, we will use the acceleration 
36 
solution to solve this PDE. By substituting $a$ directly for 
solution to solve this PDE. By substituting $a$ directly for 
37 
$\frac{\partial^{2}u_{i}}{\partial t^2}$ we can derive the 
$\frac{\partial^{2}u_{i}}{\partial t^2}$ we can derive the 
38 
displacement solution. Using $a$ we can see that \autoref{eqn:wav} becomes 
acceleration solution. Using $a$ we can see that \autoref{eqn:wav} becomes 
39 
\begin{equation} \label{eqn:wava} 
\begin{equation} \label{eqn:wava} 
40 
\rho a_{i}  \frac{\partial 
\rho a_{i}  \frac{\partial 
41 
\sigma_{ij}}{\partial x_{j}} = 0 
\sigma_{ij}}{\partial x_{j}} = 0 
42 
\end{equation} 
\end{equation} 
43 
Thus the problem will be solved for acceleration and then converted to 
Thus the problem will be solved for acceleration and then converted to 
44 
displacement using the backwards difference approximation. 
displacement using the backwards difference approximation as for the acoustic 
45 

example in the previous chapter. 
46 


47 
Consider now the stress $\sigma$. One can see that the stress consists of two 
Consider now the stress $\sigma$. One can see that the stress consists of two 
48 
distinct terms: 
distinct terms: 
59 
kronecker delta function with dimensions equivalent to $u$. The second term 
kronecker delta function with dimensions equivalent to $u$. The second term 
60 
\autoref{eqn:sigtrans} is the sum of $u$ with its own transpose. Putting these 
\autoref{eqn:sigtrans} is the sum of $u$ with its own transpose. Putting these 
61 
facts together we see that the spatial differential of the stress is given by the 
facts together we see that the spatial differential of the stress is given by the 
62 
gradient of $u$ and the aforementioned opperations. This value is then submitted 
gradient of $u$ and the aforementioned operations. This value is then submitted 
63 
to the \esc PDE as $X$. 
to the \esc PDE as $X$. 
64 
\begin{python} 
\begin{python} 
65 
g=grad(u); stress=lam*trace(g)*kmat+mu*(g+transpose(g)) 
g=grad(u); stress=lam*trace(g)*kmat+mu*(g+transpose(g)) 
69 
calculated so that the memory variables can be updated for the next time 
calculated so that the memory variables can be updated for the next time 
70 
iteration. 
iteration. 
71 
\begin{python} 
\begin{python} 
72 
accel = mypde.getSolution() #get PDE solution for accelleration 
accel = mypde.getSolution() #get PDE solution for acceleration 
73 
u_p1=(2.*uu_m1)+h*h*accel #calculate displacement 
u_p1=(2.*uu_m1)+h*h*accel #calculate displacement 
74 
u_m1=u; u=u_p1 # shift values by 1 
u_m1=u; u=u_p1 # shift values by 1 
75 
\end{python} 
\end{python} 
99 
limits the number of outputs and increases the speed of the solver. 
limits the number of outputs and increases the speed of the solver. 
100 


101 
\section{Multithreading} 
\section{Multithreading} 
102 
The wave equation solution can be quite demanding on cpu time. Enhancements can 
The wave equation solution can be quite demanding on CPU time. Enhancements can 
103 
be made by accessing multiple threads or cores on your computer. This does not 
be made by accessing multiple threads or cores on your computer. This does not 
104 
require any modification to the solution script and only comes into play when 
require any modification to the solution script and only comes into play when 
105 
escript is called from the shell. To use multiple threads \esc is called using 
\esc is called from the shell. To use multiple threads \esc is called using 
106 
the \verb!t! option with an interger argument for the number of threads 
the \verb!t! option with an integer argument for the number of threads 
107 
required. For example 
required. For example 
108 
\begin{verbatim} 
\begin{verbatim} 
109 
$escript t 4 example08a.py 
$escript t 4 example08a.py 
110 
\end{verbatim} 
\end{verbatim} 
111 
would call the script in this section and solve it using 4 threads. 
would call the script in this section and solve it using 4 threads. 
112 


113 
The computation times on an increasing number of cores is outlines in 
The computation times on an increasing number of cores is outlined in 
114 
\autoref{tab:wpcores}. 
\autoref{tab:wpcores}. 
115 


116 
\begin{table}[ht] 
\begin{table}[ht] 
136 
\section{Vector source on the boundary} 
\section{Vector source on the boundary} 
137 
\sslist{example08b.py} 
\sslist{example08b.py} 
138 
For this particular example, we will introduce the source by applying a 
For this particular example, we will introduce the source by applying a 
139 
displacment to the boundary during the initial time steps. The source will again 
displacement to the boundary during the initial time steps. The source will 
140 
be 
again be 
141 
a radially propagating wave but due to the vector nature of the PDE used, a 
a radially propagating wave but due to the vector nature of the PDE used, a 
142 
direction will need to be applied to the source. 
direction will need to be applied to the source. 
143 


160 
\end{python} 
\end{python} 
161 
where \verb xc is the source point on the boundary of the model. Note that 
where \verb xc is the source point on the boundary of the model. Note that 
162 
because the source is specifically located on the boundary, we have used the 
because the source is specifically located on the boundary, we have used the 
163 
\verb!FunctionOnBoundary! call to ensure the nodes are located upon the 
\verb!FunctionOnBoundary! call to ensure the nodes used to define the source 
164 
boundary only. These boundary nodes are passed to source as \verb!xb!. The 
are also located upon the boundary. These boundary nodes are passed to 
165 
source direction is then defined as an $(x,y)$ array and multiplied by the 
source as \verb!xb!. The source direction is then defined as an $(x,y)$ array and multiplied by the 
166 
source function. The directional array must have a magnitude of $\left 1 
source function. The directional array must have a magnitude of $\left 1 
167 
\right $ otherwise the amplitude of the source will become modified. For this 
\right $ otherwise the amplitude of the source will become modified. For this 
168 
example, the source is directed in the $y$ direction. 
example, the source is directed in the $y$ direction. 
222 
It is quite 
It is quite 
223 
simple to implement a source which is smooth in time. In addition to the 
simple to implement a source which is smooth in time. In addition to the 
224 
original source function the only extra requirement is a time function. For 
original source function the only extra requirement is a time function. For 
225 
this example the time variant source will be the derivative of a gausian curve 
this example the time variant source will be the derivative of a Gaussian curve 
226 
defined by the required dominant frequency (\autoref{fig:tvsource}). 
defined by the required dominant frequency (\autoref{fig:tvsource}). 
227 
\begin{python} 
\begin{python} 
228 
#Creating the time function of the source. 
#Creating the time function of the source. 
265 
\end{python} 
\end{python} 
266 


267 
Finally, for the length of the source, we are required to update each new 
Finally, for the length of the source, we are required to update each new 
268 
solution in the itterative section of the solver. This is done via; 
solution in the iterative section of the solver. This is done via; 
269 
\begin{python} 
\begin{python} 
270 
# increment loop values 
# increment loop values 
271 
t=t+h; n=n+1 
t=t+h; n=n+1 
278 


279 
\section{Absorbing Boundary Conditions} 
\section{Absorbing Boundary Conditions} 
280 
To mitigate the effect of the boundary on the model, absorbing boundary 
To mitigate the effect of the boundary on the model, absorbing boundary 
281 
conditions can be introduced. These conditions effectively dampen the wave 
conditions can be introduced. These conditions effectively dampen the wave energy 
282 
energy as they approach the bounday and thus prevent that energy from being 
as they approach the boundary and thus prevent that energy from being reflected. 
283 
reflected. This type of approach is used typically when a model only represents 
This type of approach is typically used when a model is shrunk to decrease 
284 
a small portion of the entire model, which in reality may have infinite bounds. 
computational requirements. In practise this applies to almost all models, 
285 
It is inpractical to calculate the solution for an infinite model and thus ABCs 
especially in earth sciences where the entire planet or a large enough 
286 
allow us the create an approximate solution with small to zero boundary effects 
portional of it cannot be modelled efficiently when considering small scale 
287 
on a model with a solvable size. 
problems. It is impractical to calculate the solution for an infinite model and thus ABCs allow 
288 

us the create an approximate solution with small to zero boundary effects on a 
289 

model with a solvable size. 
290 


291 
To dampen the waves, the method of \citet{Cerjan1985} 
To dampen the waves, the method of \citet{Cerjan1985} 
292 
where the solution and the stress are multiplied by a damping function defined 
where the solution and the stress are multiplied by a damping function defined 
348 
For stiff problems like the wave equation it is often prudent to implement 
For stiff problems like the wave equation it is often prudent to implement 
349 
second order meshing. This creates a more accurate mesh approximation with some 
second order meshing. This creates a more accurate mesh approximation with some 
350 
increased processing cost. To turn second order meshing on, the \verb!rectangle! 
increased processing cost. To turn second order meshing on, the \verb!rectangle! 
351 
function accpets an \verb!order! keyword argument. 
function accepts an \verb!order! keyword argument. 
352 
\begin{python} 
\begin{python} 
353 
domain=Rectangle(l0=mx,l1=my,n0=ndx, n1=ndy,order=2) # create the domain 
domain=Rectangle(l0=mx,l1=my,n0=ndx, n1=ndy,order=2) # create the domain 
354 
\end{python} 
\end{python} 
392 
\section{Pycad example} 
\section{Pycad example} 
393 
\sslist{example08c.py} 
\sslist{example08c.py} 
394 
To make the problem more interesting we will now introduce an interface to the 
To make the problem more interesting we will now introduce an interface to the 
395 
middle of the domain. Infact we will use the same domain as we did for heat flux 
middle of the domain. In fact we will use the same domain as we did fora 
396 
in \autoref{CHAP HEAT 2}. The domain contains a syncline with two set of 
different set of material properties on either side of the interface. 

material properties on either side of the interface. 

397 


398 
\begin{figure}[ht] 
\begin{figure}[ht] 
399 
\begin{center} 
\begin{center} 
417 
rho.setTaggedValue("top",rho1) 
rho.setTaggedValue("top",rho1) 
418 
rho.setTaggedValue("bottom",rho2) 
rho.setTaggedValue("bottom",rho2) 
419 
\end{python} 
\end{python} 
420 
Don't forget that teh source boudnary must also be tagged and added so it can be 
Don't forget that the source boundary must also be tagged and added so it can 
421 
referenced 
be referenced 
422 
\begin{python} 
\begin{python} 
423 
# Add the subdomains and flux boundaries. 
# Add the subdomains and flux boundaries. 
424 
d.addItems(PropertySet("top",tblock),PropertySet("bottom",bblock),\ 
d.addItems(PropertySet("top",tblock),PropertySet("bottom",bblock),\ 
425 
PropertySet("linetop",l30)) 
PropertySet("linetop",l30)) 
426 
\end{python} 
\end{python} 
427 
It is now possible to solve the script as in the previous examples. 
It is now possible to solve the script as in the previous examples 
428 

(\autoref{fig:ex08cres}). 
429 


430 
\begin{figure}[ht] 
\begin{figure}[ht] 
431 
\centering 
\centering 