1 

2 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
3 
% 
4 
% Copyright (c) 20032009 by University of Queensland 
5 
% Earth Systems Science Computational Center (ESSCC) 
6 
% http://www.uq.edu.au/esscc 
7 
% 
8 
% Primary Business: Queensland, Australia 
9 
% Licensed under the Open Software License version 3.0 
10 
% http://www.opensource.org/licenses/osl3.0.php 
11 
% 
12 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
13 

14 
\section{Seismic Wave Propagation in Two Dimensions} 
15 
\editor{This chapter aims to expand the readers understanding of escript by modelling the wave equations. 
16 
Challenges will include a second order differential (multiple initial conditions). A new PDE to fit to the general form. Movement to a 3D problem (maybe)??? } 
17 

18 
\esc can be used to model the propgation of waves through a medium. In this example we will see how the wave equation can be implemented using \esc and solved for in two dimension. Our domain is defined by a thin sheet that has dimensions $x$ and $y$ and to model waves we will introduce a point source displacement at our time zero. The affects of this displacement should propagate radially from the source and eventually be reflected from the boundaries of the model. 
19 

20 
The code described in this section can be found in \fileex{wavesolver2d001.py} 
21 
In a similar manner to the previous chapter the first step to creating our script is to import the necessary modules and functions. Following this the PDE and control variables must be defined. This includes the domain dimensions and type, the time scale and the time step. To ensure stability the time step can be calcuated such that it satisfies the Courant stability criteria \editor{MORE HERE ONCE METHOD FINALISED}. Considering the complexity of the computational solution to the wave equation it is proudant to consider how many steps will need to be solved. This example includes an acknowledgement clause 
22 
\begin{verbatim} 
23 
#Check to make sure number of time steps is not too large. 
24 
print "Time step size= ",h, "Expected number of outputs= ",tend/h 
25 
proceeder = raw_input("Is this ok?(y/n)") 
26 
#Exit if user thinks too many outputs. 
27 
if proceeder == "n": 
28 
sys.exit() 
29 
\end{verbatim} 
30 
This requires that the user knows the number of itterations that will be required to solve the model for the time period \verb 0 to \verb tend . The command \verb sys.exit() is used here to halt the script if the input to preceeder is \verb n and thus prevent a forced crash of the script should its projected solve time be too large. 
31 

32 
To solve this PDE we are going to introduce the concept of a python library. A library is useful as it allows a user to store defined functions that can be called to solve generic problems. The 2D wave equation satisfies this criteria. The first step is to create a new python file which we have called \verb cblib.py within this file we can set all of the necessary includes to make things easier in the future. Other advantages of libraries include a reduction in the duplication of code and the ability to modularise functions and variables. 
33 

34 
\verbwavesolver2d.py 
35 

36 
Wave propagation in the earth can be described by the wave equation: 
37 
\begin{equation} \label{eqn:wav} \index{wave equation} 
38 
\rho \frac{\partial^{2}u\hackscore{i}}{\partial t^2}  \frac{\partial \sigma _{ij}}{\partial_{j}} = 0 
39 
\end{equation} 
40 
where $\sigma$ represents stress and is given by: 
41 
\begin{equation} \label{eqn:sigw} 
42 
\sigma \hackscore{ij} = \lambda u\hackscore{k,k} \delta\hackscore{ij} + \mu ( u\hackscore{i,j} + u\hackscore{j,i}) 
43 
\end{equation} 
44 
$\lambda$ and $\mu$ are the Lame Coefficients. Specifically $\mu$ is the bulk modulus. The \ref{eqn:wav} can be written with the Einstein summation convention as: 
45 
\begin{equation} 
46 
\rho u\hackscore{j,tt} = \sigma\hackscore{ij,j} 
47 
\end{equation} 
48 
For this model we will specify the boundary conditions such that the normal of the stress from the boundary is zero. 
49 
\begin{eqnarray} \label{eqn:bdw} 
50 
\sigma \hackscore{ij}n\hackscore{j}=0 
51 
\end{eqnarray} 
52 
As the wave equation has a double time derivative, it is not sufficient to only stipulate the initial conditions for one time step. Two time steps must be specified so that the equation can be solved. For this simple example $u$ (\verb u ) and $u(t1)$ (\verb u_m1 ) will be the same but if both of these condititions are known, they can specified. It should be noted here that if multiple time steps are known or understood in the begining of a modle, they can be added to the simulation manually. The solver is then able to continue the model from where the data ends. Alternatively, if the source motion is understood, its position can be corrected for each itteration to create a more accurate recreation of an event. 
53 

54 
