1 

% 
2 
% $Id$ 
% $Id$ 
3 
% 
% 
4 
% Copyright © 2006 by ACcESS MNRF 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
5 
% \url{http://www.access.edu.au 
% 
6 
% Primary Business: Queensland, Australia. 
% Copyright 20032007 by ACceSS MNRF 
7 
% Licensed under the Open Software License version 3.0 
% Copyright 2007 by University of Queensland 
8 
% http://www.opensource.org/license/osl3.0.php 
% 
9 

% http://esscc.uq.edu.au 
10 

% Primary Business: Queensland, Australia 
11 

% Licensed under the Open Software License version 3.0 
12 

% http://www.opensource.org/licenses/osl3.0.php 
13 
% 
% 
14 

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
15 

% 
16 


17 
\section{3D Wave Propagation} 
\section{3D Wave Propagation} 
18 
\label{WAVE CHAP} 
\label{WAVE CHAP} 
19 


20 
We want to calculate the displacement field $u\hackscore{i}$ for any time $t>0$ by solving the wave equation 
In this next example we want to calculate the displacement field $u\hackscore{i}$ for any time $t>0$ by solving the wave equation 
21 
\index{wave equation}: 
\index{wave equation}: 
22 
\begin{eqnarray}\label{WAVE general problem} 
\begin{eqnarray}\label{WAVE general problem} 
23 
\rho u\hackscore{i,tt}  \sigma\hackscore{ij,j}=0 
\rho u\hackscore{i,tt}  \sigma\hackscore{ij,j}=0 
48 
set the initial displacement to be $U\hackscore 0$ in a sphere of small radius around the point 
set the initial displacement to be $U\hackscore 0$ in a sphere of small radius around the point 
49 
$x\hackscore C$. 
$x\hackscore C$. 
50 


51 
We use an explicit timeintegration scheme to calculate the displacement field $u$ at 
We use an explicit time integration scheme to calculate the displacement field $u$ at 
52 
certain time marks $t^{(n)}$ where $t^{(n)}=t^{(n1)}+h$ with time step size $h>0$. In the following the upper index ${(n)}$ refers to values at time $t^{(n)}$. We use the Verlet scheme \index{Verlet scheme} with constant time step size $h$ 
certain time marks $t^{(n)}$ where $t^{(n)}=t^{(n1)}+h$ with time step size $h>0$. In the following the upper index ${(n)}$ refers to values at time $t^{(n)}$. We use the Verlet scheme \index{Verlet scheme} with constant time step size $h$ 
53 
which is defined by 
which is defined by 
54 
\begin{eqnarray} \label{WAVE dyn 2} 
\begin{eqnarray} \label{WAVE dyn 2} 
76 
acceleration, $a^(n)$, which now depends 
acceleration, $a^(n)$, which now depends 
77 
on the solution at the previous two time steps, $u^{(n1)}$ and $u^{(n2)}$. 
on the solution at the previous two time steps, $u^{(n1)}$ and $u^{(n2)}$. 
78 


79 
In each time step we have to solve this problem to get the acceleration $a^{(n)}$ and we will 
In each time step we have to solve this problem to get the acceleration $a^{(n)}$, and we will 
80 
use the \LinearPDE class of the \linearPDEs to do so. The general form of the PDE defined through 
use the \LinearPDE class of the \linearPDEs to do so. The general form of the PDE defined through 
81 
the \LinearPDE class is discussed in \Sec{SEC LinearPDE}. The form which is relevant here are 
the \LinearPDE class is discussed in \Sec{SEC LinearPDE}. The form which is relevant here is 
82 
\begin{equation}\label{WAVE dyn 100} 
\begin{equation}\label{WAVE dyn 100} 
83 
D\hackscore{ij} a^{(n)}\hackscore{j} =  X\hackscore{ij,j}\; . 
D\hackscore{ij} a^{(n)}\hackscore{j} =  X\hackscore{ij,j}\; . 
84 
\end{equation} 
\end{equation} 
85 
Implicitly, the natural boundary condition 
The natural boundary condition 
86 
\begin{equation}\label{WAVE dyn 101} 
\begin{equation}\label{WAVE dyn 101} 
87 
n\hackscore{j}X\hackscore{ij} =0 
n\hackscore{j}X\hackscore{ij} =0 
88 
\end{equation} 
\end{equation} 
89 
is assumed. 
is used. 
90 


91 
With $u=a^{(n)}$ we can identify the values to be assigned to $D$, $X$: 
With $u=a^{(n)}$ we can identify the values to be assigned to $D$ and $X$: 
92 
\begin{equation}\label{WAVE dyn 6} 
\begin{equation}\label{WAVE dyn 6} 
93 
\begin{array}{ll} 
\begin{array}{ll} 
94 
D\hackscore{ij}=\rho \delta\hackscore{ij}& 
D\hackscore{ij}=\rho \delta\hackscore{ij}& 
103 
The following script defines a the function \function{wavePropagation} which 
The following script defines a the function \function{wavePropagation} which 
104 
implements the Verlet scheme to solve our wave propagation problem. 
implements the Verlet scheme to solve our wave propagation problem. 
105 
The argument \var{domain} which is a \Domain class object 
The argument \var{domain} which is a \Domain class object 
106 
defines the domain of the problem. \var{h} and \var{tend} are the step size 
defines the domain of the problem. \var{h} and \var{tend} are the time step size 
107 
and the time mark after which the simulation will be terminated respectively. \var{lam}, \var{mu} and 
and the end time of the simulation. \var{lam}, \var{mu} and 
108 
\var{rho} are material properties. 
\var{rho} are material properties. 
109 
\begin{python} 
\begin{python} 
110 
from esys.linearPDEs import LinearPDE 
from esys.linearPDEs import LinearPDE 
172 
u_pc_data.close() 
u_pc_data.close() 
173 
\end{python} 
\end{python} 
174 
Notice that 
Notice that 
175 
all coefficients of the PDE which are independent from time $t$ are set outside the \code{while} 
all coefficients of the PDE which are independent of time $t$ are set outside the \code{while} 
176 
loop. This allows the \LinearPDE class to reuse information as much as possible 
loop. This is very efficient since it allows the \LinearPDE class to reuse information as much as possible 
177 
when iterating over time. 
when iterating over time. 
178 


179 
We have seen most of the functions in previous examples but there are some new functions here: 
there are a few new \escript functions in this example: 
180 
\code{grad(u)} returns the gradient $u\hackscore{i,j}$ of $u$ (in fact \var{grad(g)[i,j]}=$=u\hackscore{i,j}$). 
\code{grad(u)} returns the gradient $u\hackscore{i,j}$ of $u$ (in fact \var{grad(g)[i,j]}=$=u\hackscore{i,j}$). 
181 
It is pointed out that in general there are restrictions on the argument of the \function{grad} function, for instance 
There are restrictions on the argument of the \function{grad} function, for instance 
182 
for \finley the statement \code{grad(grad(u))} will raise an exception. 
the statement \code{grad(grad(u))} will raise an exception. 
183 
\code{trace(g)} returns the sum of the main diagonal elements \var{g[k,k]} of \var{g} 
\code{trace(g)} returns the sum of the main diagonal elements \var{g[k,k]} of \var{g} 
184 
and \code{transpose(g)} returns the transposed of the matrix \var{g} (ie. 
and \code{transpose(g)} returns the matrix transpose of \var{g} (ie. $\var{transpose(g)[i,j]}=\var{g[j,i]}$). 

$\var{transpose(g)[i,j]}=\var{g[j,i]}$). 

185 


186 
We initialize the values of \code{u} and \code{u_last} to be $U\hackscore 0$ for a small 
We initialize the values of \code{u} and \code{u_last} to be $U\hackscore 0$ for a small 
187 
sphere of radius, \code{src_radius} around the point source located at $x\hackscore C$. 
sphere of radius, \code{src_radius} around the point source located at $x\hackscore C$. 
199 
%\end{python} 
%\end{python} 
200 
%declares \var{u} as a vectorvalued function of its coordinates 
%declares \var{u} as a vectorvalued function of its coordinates 
201 
%in the domain (i.e. with two or three components in the case of 
%in the domain (i.e. with two or three components in the case of 
202 
%a two or three dimensional domain, repectively). The first argument defines the value of the function which is equal 
%a two or three dimensional domain, respectively). The first argument defines the value of the function which is equal 
203 
%to $0$ for all vectot components. The second argument \var{ContinuousFunction(domain)} 
%to $0$ for all vectot components. The second argument \var{ContinuousFunction(domain)} 
204 
%specifies the `smoothness` of the function. In our case we want the displacement field to be 
%specifies the `smoothness` of the function. In our case we want the displacement field to be 
205 
%a continuous function over the \Domain \var{domain}. On other option would be 
%a continuous function over the \Domain \var{domain}. On other option would be 
206 
%\var{Function(domain)} in which case continuity is not garanteed. In fact the 
%\var{Function(domain)} in which case continuity is not guaranteed. In fact the 
207 
%argument of the \function{grad} function has to be of the typ \var{ContinuousFunction(domain)} while 
%argument of the \function{grad} function has to be of the typ \var{ContinuousFunction(domain)} while 
208 
%returned function is of the var{Function(domain)} type. For more details see \Sec{SEC ESCRIPT DATA}. 
%returned function is of the var{Function(domain)} type. For more details see \Sec{SEC ESCRIPT DATA}. 
209 


210 
The statement \code{myPDE.setLumpingOn()} 
The statement \code{myPDE.setLumpingOn()} 
211 
switches on the usage of an aggressive approximation of the PDE operator, in this case 
switches on the use of an aggressive approximation of the PDE operator as a diagonal matrix 
212 
formed by the coefficient \var{D} (actually the discrete 
formed from the coefficient \var{D}. 
213 
version of the PDE operator is approximated by a diagonal matrix). The approximation allows, at the costs of 
The approximation allows, at the cost of 
214 
an additional error, very fast 
additional error, very fast 
215 
solving of the PDE when the solution is requested. In the general case, the presence of \var{A}, \var{B} or \var 
solution of the PDE. When using \code{setLumpingOn} the presence of \var{A}, \var{B} or \var 
216 
{C} \code{setLumpingOn} will produce wrong results. 
{C} will produce wrong results. 
217 


218 


219 
One of the big advantages of the Verlet scheme is the fact that the problem to be solved 
One of the big advantages of the Verlet scheme is the fact that the problem to be solved 
220 
in each time step is very simple and does not involve any spatial derivatives (which actually allows us to use lumping). 
in each time step is very simple and does not involve any spatial derivatives (which is what allows us to use lumping in this simulation). 
221 
The problem becomes so simple because we use the stress from the last time step rather then the stress which is 
The problem becomes so simple because we use the stress from the last time step rather then the stress which is 
222 
actually present at the current time step. Schemes using this approach are called an explicit time integration 
actually present at the current time step. Schemes using this approach are called an explicit time integration 
223 
scheme \index{explicit scheme} \index{time integration!explicit}. The 
schemes \index{explicit scheme} \index{time integration!explicit}. The 
224 
backward Euler scheme we have used in \Chap{DIFFUSION CHAP} is 
backward Euler scheme we have used in \Chap{DIFFUSION CHAP} is 
225 
an example of an implicit schemes 
an example of an implicit scheme 
226 
\index{implicit scheme} \index{time integration!implicit}. In this case one uses the actual status of 
\index{implicit scheme} \index{time integration!implicit}. In this case one uses the actual status of 
227 
each variable at a particular time rather then values from previous time steps. This will lead to a problem which is 
each variable at a particular time rather then values from previous time steps. This will lead to a problem which is 
228 
more expensive to solve, in particular for nonlinear problems. 
more expensive to solve, in particular for nonlinear problems. 