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. The \refEq{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 at the begining of a model, they can be added to the simulation manually. The solver can then continue the model from the end of the known data. Alternatively, if the source motion is understood, its position can be corrected for each iteration 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 stepping 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 of unit length such that $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 limitations 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 in our source directionality. This still introduced assumptions and removes realistic wave motions both P and S 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 two 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 by using a cosine taper. Our source terms then become; 
73 
\begin{verbatim} 
74 
u=U0*(cos(length(xxc)*3.1415/src_radius)+1)*\ 
75 
whereNegative(length(xxc)src_radius)*dunit 
76 
u_m1=u 
77 
\end{verbatim} 
78 

79 
It is often useful to know the values of PDE at certain locations in the model. To acheive this we are going to use a new generic function called \verb cbphones which allows us to specify receiver locations to record the PDE values at those points. The function \verb cbphones as the arguments; 
80 
\begin{verbatim} 
81 
# domain : domain of model 
82 
# U : Current time state displacement solution. 
83 
# phones : Geophone Locations 
84 
# dim : model dimesions 
85 
# savepath: where to output the data files local is default 
86 
\end{verbatim} 
87 
\editor{not generic as of yet but may move to make cbphones and the phones positioning a serious part of wavesolver 2d} 
88 
We have chosen to have three receivers and they are called using; 
89 
\begin{verbatim} 
90 
u_pot = cbphones(domain,u,[[0,500],[250,500],[400,500]],2) 
91 
\end{verbatim} 
92 
This places the receivers on the surface at the source location and two locations further along the top of the model. The output \verb u_pot can then be split and saved to file using the following command; 
93 
\begin{verbatim} 
94 
u_pc_data=open(os.path.join(savepath,'U_pc.out'),'w') 
95 
u_pc_data.write("%f %f %f %f %f %f %f\n"%(t,u_pc_x1,u_pc_y1,u_pc_x2,u_pc_y2,u_pc_x3,u_pc_y3)) 
96 
\end{verbatim} 
97 
Convieniently this saves the time, x direction displacement and y direction displacement values for these locations. Now that the initial conditions have been defined we can tackle the task of solving the wave equation for the number of required time steps. To do this we require a while loop and form of the wave equation which fits our general linear PDE form. We start with the form of the equation for stress \refEq{eqn:sigw}. We can define the kronecker matrix using the domain and take the derivative of \verb u via the function \verbgrad(u) . As $\lambda$ and $\mu$ are constants we can now define $\sigma$; 
98 
\begin{verbatim} 
99 
g=grad(u) 
100 
stress=lam*trace(g)*kmat+mu*(g+transpose(g)) 
101 
\end{verbatim} 
102 
Solving for the double time derivative of u on the LHS of \refEq{eqn:wav} required us to use the centred difference forumlua which returns; 
103 
\begin{equation} 
104 
u^n = 2u^{n1}u^{n2}+h^2 \biggl(\frac{\partial ^2 u}{\partial t^2}\biggr)^n 
105 
\end{equation} 
106 
Substituting for the double time derivative we see; 
107 
\begin{equation} 
108 
\rho u^n = 2\rho u^{n1} \rho u^{n2}+h^2 \sigma \hackscore{ij,j} ^n 
109 
\end{equation} 
110 
This fits the general form $Du=X \hackscore{j,j} + Y$ where $D=rho$; $Y=2\rho u^{n1} \rho u^{n2}$ and $X=h^2 \sigma \hackscore{ij,j} ^n$. \verb D does not vary between time steps can be defined before our iteration loop via; 
111 
\begin{verbatim} 
112 
mypde.setValue(D=kmat*rho) 
113 
\end{verbatim} 
114 
The values for \verb u must be refreshed after each iteration and are thus defined within our while loop via; 
115 
\begin{verbatim} 
116 
mypde.setValue(X=stress*(h*h),Y=(rho*2*urho*u_m1)) 
117 
\end{verbatim} 
118 
With each iteration we update \verb u and migrate our answers into the correct variables. The iterative values must also be updated as well as the response from our receiver locations. All this is acheived via; 
119 
\begin{verbatim} 
120 
u_p1 = mypde.getSolution() 
121 
u_m1=u 
122 
u=u_p1 
123 
t+=h 
124 
n+=1 
125 
u_pot = cbphones(domain,u,[[125.,250.],[250.,250.],[250.,375.]],2) 
126 
# save displacements at point source to file for t > 0 
127 
u_pc_data.write("%f %f %f %f %f %f %f\n"%(t,u_pc_x1,u_pc_y1,u_pc_x2,u_pc_y2,u_pc_x3,u_pc_y3)) 
128 
\end{verbatim} 
129 
With an appropriate file saving output we now have a working generic solver for our 2D wave equation problem. We have also included our two new generic programs in the \verb cblib library so they can be more simply imported along with their own dependancies to our test script. 
130 

131 
Writing a test program will allow us to more easily pass the variables required to the solver to generate an output solution. Our testing code described in this section can be found in \fileex{wavesolver2d001.py}. In a similar manner to the previous chapter the first step to creating our script is to import the necessary modules and functions including our new library file. 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. Our test script thus includes an acknowledgement clause; 
132 
\begin{verbatim} 
133 
#Check to make sure number of time steps is not too large. 
134 
print "Time step size= ",h, "Expected number of outputs= ",tend/h 
135 
proceeder = raw_input("Is this ok?(y/n)") 
136 
#Exit if user thinks too many outputs. 
137 
if proceeder == "n": 
138 
sys.exit() 
139 
\end{verbatim} 
140 
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. 