12 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
13 


14 
\section{Seismic Wave Propagation in Two Dimensions} 
\section{Seismic Wave Propagation in Two Dimensions} 
15 
\editor{This chapter aims to expand the readers understanding of escript by modelling the wave equations. 
\editor{This chapter aims to expand the readers understanding of escript by 
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)??? } 
modelling the wave equations. 
17 

Challenges will include a second order differential (multiple initial 
18 

conditions). A new PDE to fit to the general form. Movement to a 3D problem 
19 

(maybe)??? } 
20 


21 

\sslist{example08a.py} 
22 


23 

We will now expand upon the previous chapter by introducing a vector form of 
24 

the wave equation. This means that the waves will have both a scalar magnitude, 
25 

but also a direction. This type of scenario is apparent in wave forms that 
26 

exhibit compressional and transverse particle motion. A common type of wave 
27 

that obeys this principle are seismic waves. 
28 


29 
\verbwavesolver2d.py 
Wave propagation in the earth can be described by the elastic wave equation: 




Wave propagation in the earth can be described by the wave equation: 

30 
\begin{equation} \label{eqn:wav} \index{wave equation} 
\begin{equation} \label{eqn:wav} \index{wave equation} 
31 
\rho \frac{\partial^{2}u\hackscore{i}}{\partial t^2}  \frac{\partial \sigma\hackscore{ij}}{\partial\hackscore{j}} = 0 
\rho \frac{\partial^{2}u\hackscore{i}}{\partial t^2}  \frac{\partial 
32 

\sigma\hackscore{ij}}{\partial x\hackscore{j}} = 0 
33 
\end{equation} 
\end{equation} 
34 
where $\sigma$ represents stress and is given by: 
where $\sigma$ is the stress given by: 
35 
\begin{equation} \label{eqn:sigw} 
\begin{equation} \label{eqn:sigw} 
36 
\sigma \hackscore{ij} = \lambda u\hackscore{k,k} \delta\hackscore{ij} + \mu ( u\hackscore{i,j} + u\hackscore{j,i}) 
\sigma \hackscore{ij} = \lambda u\hackscore{k,k} \delta\hackscore{ij} + \mu ( 
37 

u\hackscore{i,j} + u\hackscore{j,i}) 
38 
\end{equation} 
\end{equation} 
39 
$\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: 
where $\lambda$ and $\mu$ are the Lame Coefficients. Specifically $\mu$ is the 
40 

bulk modulus. The \refEq{eqn:wav} can be written with the Einstein summation 
41 

convention as: 
42 
\begin{equation} 
\begin{equation} 
43 
\rho u\hackscore{i,tt} = \sigma\hackscore{ij,j} 
\rho u\hackscore{i,tt} = \sigma\hackscore{ij,j} 
44 
\end{equation} 
\end{equation} 

For this model we will specify the boundary conditions such that the normal of the stress from the boundary is zero. 


\begin{eqnarray} \label{eqn:bdw} 


\sigma \hackscore{ij}n\hackscore{j}=0 


\end{eqnarray} 


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. 





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; 


\begin{python} 


domain : domain to solve over 


h : time step 


tend : end time 


lam, mu : lames constants for domain 


rho : density of domain 


U0 : magnitude of source 


xc : source location in domain (Vector) 


savepath: where to output the data files 


\end{python} 


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; 


\begin{python} 


mypde.setSolverMethod(LinearPDE.LUMPING) 


\end{python} 


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. 





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. 





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; 


\begin{python} 


src_radius = 50 


\end{python} 


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; 


\begin{enumerate} 


\item Along the xaxis: $U=[dx=1,dy=0]$ 


\item Along the yaxis: $U=[dx=0,dy=1]$ 


\item At 45deg: $U=[dx=\frac{1}{\sqrt2},dy=\frac{1}{\sqrt2}]$ 


\end{enumerate} 


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. 


For our example we will use; 


\begin{python} 


dunit=numarray.array([1.,0.]) 


\end{python} 


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; 


\begin{python} 


u=U0*(cos(length(xxc)*3.1415/src_radius)+1)*\ 


whereNegative(length(xxc)src_radius)*dunit 


u_m1=u 


\end{python} 

45 


46 
It is often useful to know the values of PDE at certain locations in the model. 
In a similar process to the previous chapter, we will use the acceleration 
47 
To acheive this we are going to use a new generic function 
solution to solve this PDE. By substituting $a$ directly for 
48 
called \verb cbphones which allows us to specify receiver locations to record 
$\frac{\partial^{2}u\hackscore{i}}{\partial t^2}$ we can derive the 
49 
the PDE values at those points. The function \verb cbphones as the arguments; 
displacement solution. Using $a$ \refEq{eqn:wav} becomes; 
50 
\begin{python} 
\begin{equation} \label{eqn:wava} 
51 
# domain : domain of model 
\rho a\hackscore{i}  \frac{\partial 
52 
# U : Current time state displacement solution. 
\sigma\hackscore{ij}}{\partial x\hackscore{j}} = 0 

# phones : Geophone Locations 


# dim : model dimesions 


# savepath: where to output the data files local is default 


\end{python} 


\editor{not generic as of yet but may move to make cbphones and the phones 


positioning a serious part of wavesolver 2d} 


We have chosen to have three receivers and they are called using; 


\begin{python} 


u_pot = cbphones(domain,u,[[0,500],[250,500],[400,500]],2) 


\end{python} 


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; 


\begin{python} 


u_pc_data=open(os.path.join(savepath,'U_pc.out'),'w') 


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)) 


\end{python} 


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$; 


\begin{python} 


g=grad(u) 


stress=lam*trace(g)*kmat+mu*(g+transpose(g)) 


\end{python} 


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; 


\begin{equation} 


u^n = 2u^{n1}u^{n2}+h^2 \biggl(\frac{\partial ^2 u}{\partial t^2}\biggr)^n 


\end{equation} 


Substituting for the double time derivative we see; 


\begin{equation} 


\rho u^n = 2\rho u^{n1} \rho u^{n2}+h^2 \sigma \hackscore{ij,j} ^n 

53 
\end{equation} 
\end{equation} 
54 
This fits the general form $Du=X \hackscore{j,j} + Y$ where $D=rho$; $Y=2\rho 

55 
u^{n1} \rho u^{n2}$ and $X=h^2 \sigma \hackscore{ij,j} ^n$. \verb D does not 
\section{Vector source on the boundary} 
56 
vary between time steps can be defined before our iteration loop via; 
For this particular example, we will introduce the source by applying a 
57 
\begin{python} 
displacment to the boundary during the initial time steps. The source will again 
58 
mypde.setValue(D=kmat*rho) 
be 
59 
\end{python} 
a radially propagating wave but due to the vector nature of the PDE used, a 
60 
The values for \verb u must be refreshed after each iteration and are thus 
direction will need to be applied to the source. 
61 
defined within our while loop via; 

62 
\begin{python} 
The first step is to choose an amplitude and create the source as in the 
63 
mypde.setValue(X=stress*(h*h),Y=(rho*2*urho*u_m1)) 
previous chapter. 
64 
\end{python} 
\begin{python} 
65 
With each iteration we update \verb u and migrate our answers into the correct 
src_length = 20; print "src_length = ",src_length 
66 
variables. The iterative values must also be updated as well as the response 
# set initial values for first two time steps with source terms 
67 
from our receiver locations. All this is acheived via; 
y=U0*(cos(length(xxc)*3.1415/src_length)+1)*whereNegative(length(xxc)src_leng 
68 
\begin{python} 
th) 
69 
u_p1 = mypde.getSolution() 
\end{python} 
70 

where \verb xc is the source point on the boundary of the model. The source 
71 

direction is then defined as an $(x,y)$ array and multiplied by the source 
72 

function. The directional array must have a magnitude of $1$ otherwise the 
73 

amplitude of the source will become modified. For this example, the source is 
74 

directed in the $y$ direction. 
75 

\begin{python} 
76 

src_dir=numpy.array([0.,1.]) # defines direction of point source as down 
77 

y=y*src_dir 
78 

\end{python} 
79 

The function can then be applied as a boundary condition by setting it equal to 
80 

$y$ in the general form. 
81 

\begin{python} 
82 

mypde.setValue(y=y) #set the source as a function on the boundary 
83 

\end{python} 
84 

Because we are no longer using the source to define our initial condition to 
85 

the model, we must set the model state to zero for the first two time steps. 
86 

\begin{python} 
87 

# initial value of displacement at point source is constant (U0=0.01) 
88 

# for first two time steps 
89 

u=[0.0,0.0]*whereNegative(x) 
90 
u_m1=u 
u_m1=u 

u=u_p1 


t+=h 


n+=1 


u_pot = cbphones(domain,u,[[125.,250.],[250.,250.],[250.,375.]],2) 


# save displacements at point source to file for t > 0 


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)) 


\end{python} 


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. 





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; 


\begin{python} 


#Check to make sure number of time steps is not too large. 


print "Time step size= ",h, "Expected number of outputs= ",tend/h 


proceeder = raw_input("Is this ok?(y/n)") 


#Exit if user thinks too many outputs. 


if proceeder == "n": 


sys.exit() 

91 
\end{python} 
\end{python} 
92 
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. 

93 

If the source will introduce energy to the system over a period longer than one 
94 

or two time steps (ie the initial conditions), $y$ can be updated during the 
95 

iteration stage. 