1 

2 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
3 
% 
4 
% Copyright (c) 20032010 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 
16 
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 
Wave propagation in the earth can be described by the elastic wave equation: 
30 
\begin{equation} \label{eqn:wav} \index{wave equation} 
31 
\rho \frac{\partial^{2}u\hackscore{i}}{\partial t^2}  \frac{\partial 
32 
\sigma\hackscore{ij}}{\partial x\hackscore{j}} = 0 
33 
\end{equation} 
34 
where $\sigma$ is the stress given by: 
35 
\begin{equation} \label{eqn:sigw} 
36 
\sigma \hackscore{ij} = \lambda u\hackscore{k,k} \delta\hackscore{ij} + \mu ( 
37 
u\hackscore{i,j} + u\hackscore{j,i}) 
38 
\end{equation} 
39 
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} 
43 
\rho u\hackscore{i,tt} = \sigma\hackscore{ij,j} 
44 
\end{equation} 
45 

46 
In a similar process to the previous chapter, we will use the acceleration 
47 
solution to solve this PDE. By substituting $a$ directly for 
48 
$\frac{\partial^{2}u\hackscore{i}}{\partial t^2}$ we can derive the 
49 
displacement solution. Using $a$ \refEq{eqn:wav} becomes; 
50 
\begin{equation} \label{eqn:wava} 
51 
\rho a\hackscore{i}  \frac{\partial 
52 
\sigma\hackscore{ij}}{\partial x\hackscore{j}} = 0 
53 
\end{equation} 
54 

55 
\section{Vector source on the boundary} 
56 
For this particular example, we will introduce the source by applying a 
57 
displacment to the boundary during the initial time steps. The source will again 
58 
be 
59 
a radially propagating wave but due to the vector nature of the PDE used, a 
60 
direction will need to be applied to the source. 
61 

62 
The first step is to choose an amplitude and create the source as in the 
63 
previous chapter. 
64 
\begin{python} 
65 
src_length = 20; print "src_length = ",src_length 
66 
# set initial values for first two time steps with source terms 
67 
y=U0*(cos(length(xxc)*3.1415/src_length)+1)*whereNegative(length(xxc)src_leng 
68 
th) 
69 
\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 
91 
\end{python} 
92 

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. 