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 
\verbwavesolver2d.py 
19 

20 
Wave propagation in the earth can be described by the wave equation: 
21 
\begin{equation} \label{eqn:wav} \index{wave equation} 
22 
\rho \frac{\partial^{2}u\hackscore{i}}{\partial t^2}  \frac{\partial \sigma\hackscore{ij}}{\partial\hackscore{j}} = 0 
23 
\end{equation} 
24 
where $\sigma$ represents stress and is given by: 
25 
\begin{equation} \label{eqn:sigw} 
26 
\sigma \hackscore{ij} = \lambda u\hackscore{k,k} \delta\hackscore{ij} + \mu ( u\hackscore{i,j} + u\hackscore{j,i}) 
27 
\end{equation} 
28 
$\lambda$ and $\mu$ are the Lame Coefficients. Specifically $\mu$ is the bulk modulus. Equation \ref{eqn:wav} can be written with the Einstein summation convention as: 
29 
\begin{equation} 
30 
\rho u\hackscore{i,tt} = \sigma\hackscore{ij,j} 
31 
\end{equation} 
32 
For this model we will specify the boundary conditions such that the normal of the stress from the boundary is zero. 
33 
\begin{eqnarray} \label{eqn:bdw} 
34 
\sigma \hackscore{ij}n\hackscore{j}=0 
35 
\end{eqnarray} 
36 
To solve this PDE we are going to write a more generic solution routine than we have in the previous chapters. There are a number of advantages to this approach. Firstly, by writing a subroutine that will solve a 2D wave propagation problem it reduces the amount of code duplication that may occur. When errors arrise one need only ammend the subroutine rather than all iterations of it that may have been created. This saves time and effort in the long run. 
37 

38 
To create the our subroutine we will need to import all our necessary libraries again. This is as for previous examples. Then we can define our script and the variables it will take. Our subroutine is located in \verb/examples/cblib/wavesolver2d.py . The arguments of the subroutine are; 
39 
\begin{verbatim} 
40 
domain : domain to solve over 
41 
h : time step 
42 
tend : end time 
43 
lam, mu : lames constants for domain 
44 
rho : density of domain 
45 
U0 : magnitude of source 
46 
xc : source location in domain (Vector) 
47 
savepath: where to output the data files 
48 
\end{verbatim} 
49 
There are a few differences between the wave equation and the heat diffusion problem of the previous chapter. While the nodes are defined the same way with the function \verb getX and the PDE is still linear; one must consider the solution method. Without the line; 
50 
\begin{verbatim} 
51 
mypde.setSolverMethod(LinearPDE.LUMPING) 
52 
\end{verbatim} 
53 
the PDE would take a significant amount of time to solve. The \verb LUMPING functionality implements an aggressive approximation for the $D$ coefficient matrix of the \esc linear PDE general form. While \verb LUMPING introduces additional error to the solution it can significantly reduce the solution time. Care should be taken however, as this function can only be used when the $A$, $B$ and $C$ coefficients of the general form are zero. 
54 

55 
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 example $u$ (\verb u ) and $u(t1)$ (\verb u_m1 ) will be the same but if both of these condititions are known, they can be specified individually. It should be noted here that if multiple time steps are known for the begining of a model, they can be added to the simulation manually. The solver can then continue the model from the known data. Alternatively, if the source motion is understood, its position can be corrected for each itteration to create a more accurate recreation of an event. 
56 

57 
The source in this example will induce a radially propagating wave. A small displacement will be applied to the medium about a singularity which we have called \verb xc , this is the source location. We start by giving the source some spatial magnitude by defining a small radius about \verb xc which is affected. The \verb src_radius needs to cover a significant portion of grid nodes, otherwise the waves generated will suffer from dispersion due to an inadequate grid step size. If the source is small, the grid steping must reflect the size of the source for more accurate results. Our radius will be; 
58 
\begin{verbatim} 
59 
src_radius = 50 
60 
\end{verbatim} 
61 
Now that the extent of the source has been allocated it needs two more things; a direction and a magnitude. We can choose a direction based on the 360 degrees that exist in a full circle. If we take $\theta=0$ to be the xaxis and move counter clockwise then we can create a directional vector $U=[dx,dy]$ where $tan(\theta) = dy/dx$. It is also necessary to ensure that our directional vector is a unit $U=1$ which implies $\sqrt{dx^2+dy^2}=1$. By doing this we ensure that no accidental scaling is introduced to our source term. Here are three examples of different directions which satisfy the above conditions; 
62 
\begin{enumerate} 
63 
\item Along the xaxis: $U=[dx=1,dy=0]$ 
64 
\item Along the yaxis: $U=[dx=0,dy=1]$ 
65 
\item At 45deg: $U=[dx=\frac{1}{\sqrt2},dy=\frac{1}{\sqrt2}]$ 
66 
\end{enumerate} 
67 
There are limitation to specifying the source in this manner. Realistically we would not expect a 2D surface source to move form side to side as an isotropic source makes more sense. \editor{I am not sure here how to create an isotropic source function.}. In the 3D case things are not quite so bad. Normally we are interested in the pwaves that are directed dowwards and thus we need not have any x or y component to our source directionality. This still introduced assumptions and removes realistic wave motions both s and p from the model. 
68 
For our example we will use; 
69 
\begin{verbatim} 
70 
dunit=numarray.array([1.,0.]) 
71 
\end{verbatim} 
72 
Next we must define the values of our entire domain for the first and second time step. For the purposes of this example it is sufficient to have these to timesteps as equal. Setting the source is similar to earlier problems where we can use \esc functions to set specific areas of the domain to certain values. We must also smooth our sourse to its surrounds to prevent ?diffusion? errors. This is acheived using a cosine taper. Our source terms then become; 
73 
\begin{verbatim} 
74 
u=U0*(cos(length(xxc)*3.1415/src_radius)+1)*whereNegative(length(xxc)src_radius)*dunit 
75 
u_m1=u 
76 
\end{verbatim} 
77 

78 
Write now about phones.py, work on fitting form of wave equation to source terms. 
79 

80 

81 
In this example we will see how the wave equation can be implemented using \esc and solved for in two dimensions. 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 time zero. The affects of this displacement should propagate radially from the source and eventually be reflected from the boundaries of the model. 
82 

83 

84 

85 

86 

87 

88 
The code described in this section can be found in \fileex{wavesolver2d001.py} 
89 
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 
90 
\begin{verbatim} 
91 
#Check to make sure number of time steps is not too large. 
92 
print "Time step size= ",h, "Expected number of outputs= ",tend/h 
93 
proceeder = raw_input("Is this ok?(y/n)") 
94 
#Exit if user thinks too many outputs. 
95 
if proceeder == "n": 
96 
sys.exit() 
97 
\end{verbatim} 
98 
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. 
99 

100 
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. 
101 

102 

103 
