doc/user/wave.tex


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% Copyright (c) 2003-2008 by University of Queensland
% Earth Systems Science Computational Center (ESSCC)
% http://www.uq.edu.au/esscc
%
% Primary Business: Queensland, Australia
% Licensed under the Open Software License version 3.0
% http://www.opensource.org/licenses/osl-3.0.php
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


\section{3-D Wave Propagation}
\label{WAVE CHAP}

In this next example we want to calculate the displacement field $u\hackscore{i}$ for any time $t>0$ by solving the wave equation:
\index{wave equation}
\begin{eqnarray}\label{WAVE general problem}
\rho u\hackscore{i,tt} - \sigma\hackscore{ij,j}=0
\end{eqnarray}
in a three dimensional block of length $L$ in $x\hackscore{0}$
and $x\hackscore{1}$ direction and height $H$
in $x\hackscore{2}$ direction. $\rho$ is the known density which may be a function of its location.
$\sigma\hackscore{ij}$ is the stress field \index{stress} which in case of an isotropic, linear elastic material is given by
\begin{eqnarray} \label{WAVE stress}
\sigma\hackscore{ij} & = & \lambda u\hackscore{k,k} \delta\hackscore{ij} + \mu ( u\hackscore{i,j} + u\hackscore{j,i})
\end{eqnarray}
where $\lambda$ and $\mu$ are the Lame coefficients 
\index{Lame coefficients} and $\delta\hackscore{ij}$ denotes the Kronecker symbol\index{Kronecker symbol}.
On the boundary the normal stress is given by
\begin{eqnarray} \label{WAVE natural}
\sigma\hackscore{ij}n\hackscore{j}=0
\end{eqnarray}
for all time $t>0$. 

At initial time $t=0$ the displacement 
$u\hackscore{i}$ and the velocity $u\hackscore{i,t}$ are given:
\begin{eqnarray} \label{WAVE initial}
u\hackscore{i}(0,x)= \left\{\begin{array}{cc} U\hackscore{0} & \mbox{for $x$ at point charge, $x\hackscore{C}$} \\ 0 & \mbox{elsewhere}\end{array}  \right. & \mbox{ and } & u\hackscore{i,t}(0,x)=0 \; 
\end{eqnarray}
for all $x$ in the domain. 

Here we are modelling a point source at the point $x\hackscore C$, in the numerical solution we
 set the initial displacement to be $U\hackscore 0$ in a sphere of small radius around the point
 $x\hackscore C$.  

We use an explicit time integration scheme to calculate the displacement field $u$ at 
certain time marks $t^{(n)}$ where $t^{(n)}=t^{(n-1)}+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$
which is defined by
\begin{eqnarray} \label{WAVE dyn 2}
u^{(n)}=2u^{(n-1)}-u^{(n-2)} + h^2 a^{(n)} \\
\end{eqnarray}
for all $n=2,3,\ldots$. It is designed to solve a system of equations of the form
\begin{eqnarray} \label{WAVE dyn 2b} 
u\hackscore{,tt}=G(u)
\end{eqnarray}
where one sets $a^{(n)}=G(u^{(n-1)})$.

In our case $a^{(n)}$ is given by
\begin{eqnarray}\label{WAVE dyn 3}
\rho a^{(n)}\hackscore{i}=\sigma^{(n-1)}\hackscore{ij,j}
\end{eqnarray}
and boundary conditions
\begin{eqnarray} \label{WAVE natural at n}
\sigma^{(n-1)}\hackscore{ij}n\hackscore{j}=0
\end{eqnarray}
derived from \eqn{WAVE natural} where 
\begin{eqnarray} \label{WAVE dyn 3a}
\sigma\hackscore{ij}^{(n-1)} & = & \lambda u^{(n-1)}\hackscore{k,k} \delta\hackscore{ij} + \mu ( u^{(n-1)}\hackscore{i,j} + u^{(n-1)}\hackscore{j,i}).
\end{eqnarray}
Now we have converted our problem for displacement, $u^{(n)}$, into a problem for 
acceleration, $a^(n)$, which now depends 
on the solution at the previous two time steps, $u^{(n-1)}$  and $u^{(n-2)}$.

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

With $u=a^{(n)}$ we can identify the values to be assigned to $D$ and $X$:
\begin{equation}\label{WAVE  dyn 6}
\begin{array}{ll}
D\hackscore{ij}=\rho \delta\hackscore{ij}&
X\hackscore{ij}=-\sigma^{(n-1)}\hackscore{ij} \\
\end{array}
\end{equation}


%Remember that \var{y.whereZero()} returns a function which is $1$ for those points $x$ where $y$ is zero
%(up to a certain tolerance, see \Sec{SEC ESCRIPT DATA}) and $0$ for those points where $y$ is not equal to zero.
 
The following script defines a the function \function{wavePropagation} which
implements the Verlet scheme to solve our wave propagation problem. 
The argument \var{domain} which is a \Domain class object
defines the domain of the problem. \var{h} and \var{tend} are the time step size
and the end time of the simulation. \var{lam}, \var{mu} and 
\var{rho} are material properties. 
\begin{python}
from esys.linearPDEs import LinearPDE
from numarray import identity,zeros,ones
from esys.escript import *
from esys.escript.pdetools import Locator
def wavePropagation(domain,h,tend,lam,mu,rho,U0):
   x=domain.getX()
   # ... open new PDE ...
   mypde=LinearPDE(domain)
   mypde.setSolverMethod(LinearPDE.LUMPING)
   kronecker=identity(mypde.getDim())
   #  spherical source at middle of bottom face
   xc=[width/2.,width/2.,0.] 
   # define small radius around point xc
   # Lsup(x) returns the maximum value of the argument x
   src_radius = 0.1*Lsup(domain.getSize())
   print "src_radius = ",src_radius
   dunit=numarray.array([1.,0.,0.]) # defines direction of point source
   mypde.setValue(D=kronecker*rho)
   # ... set initial values ....
   n=0
   # initial value of displacement at point source is constant (U0=0.01)
   # for first two time steps
   u=U0*whereNegative(length(x-xc)-src_radius)*dunit
   u_last=U0*whereNegative(length(x-xc)-src_radius)*dunit
   t=0
   # define the location of the point source
   L=Locator(domain,xc)
   # find potential at point source
   u_pc=L.getValue(u)
   print "u at point charge=",u_pc
   u_pc_x = u_pc[0]
   u_pc_y = u_pc[1]
   u_pc_z = u_pc[2]
   # open file to save displacement at point source
   u_pc_data=open('./data/U_pc.out','w')
   u_pc_data.write("%f %f %f %f\n"%(t,u_pc_x,u_pc_y,u_pc_z))
   while t<tend:
     # ... get current stress ....
     g=grad(u)
     stress=lam*trace(g)*kronecker+mu*(g+transpose(g))
     # ... get new acceleration ....
     mypde.setValue(X=-stress)
     a=mypde.getSolution()
     # ... get new displacement ...
     u_new=2*u-u_last+h**2*a
     # ... shift displacements ....
     u_last=u
     u=u_new
     t+=h
     n+=1
     print n,"-th time step t ",t
     L=Locator(domain,xc)
     u_pc=L.getValue(u)
     print "u at point charge=",u_pc
     u_pc_x=u_pc[0]
     u_pc_y=u_pc[1]
     u_pc_z=u_pc[2]
     # save displacements at point source to file for t > 0
     u_pc_data.write("%f %f %f %f\n"%(t,u_pc_x,u_pc_y,u_pc_z))
     # ... save current acceleration in units of gravity and displacements
     if n==1 or n%10==0: saveVTK("./data/usoln.%i.vtu"%(n/10),acceleration=length(a)/9.81,
     displacement = length(u), Ux = u[0] )
   u_pc_data.close()
\end{python}
Notice that 
all coefficients of the PDE which are independent of time $t$ are set outside the \code{while} 
loop. This is very efficient since it allows the \LinearPDE class to reuse information as much as possible 
when iterating over time.  
 
there are a few new \escript functions in this example: 
\code{grad(u)} returns the gradient $u\hackscore{i,j}$ of $u$ (in fact \var{grad(g)[i,j]}=$=u\hackscore{i,j}$).
There are restrictions on the argument of the \function{grad} function, for instance
the statement \code{grad(grad(u))} will raise an exception.
\code{trace(g)} returns the sum of the main diagonal elements \var{g[k,k]} of \var{g} 
and \code{transpose(g)} returns the matrix transpose of \var{g} (ie. $\var{transpose(g)[i,j]}=\var{g[j,i]}$). 

We initialize the values of \code{u} and \code{u_last} to be $U\hackscore 0$ for a small 
sphere of radius, \code{src_radius} around the point source located at $x\hackscore C$.
These satisfy the initial conditions defined in \eqn{WAVE initial}.

The output \code{U_pc.out} and \code{vtu} files are saved in a directory called \code{data}.
The \code{U_pc.out} stores 4 columns of data: $t,u\hackscore x,u\hackscore y,u\hackscore z$ 
respectively, where $u\hackscore x,u\hackscore y,u\hackscore z$ are the $x\hackscore 0,x\hackscore 1,x\hackscore 2$ components of 
the displacement vector $u$ at the point source. These can be
plotted easily using any plotting package. In gnuplot the command:
\code{plot 'U_pc.out' u 1:2 title 'U_x' w l lw 2, 'U_pc.out' u 1:3 title 'U_y' w l lw 2, 
'U_pc.out' u 1:4 title 'U_z' w l lw 2} will reproduce Figure \ref{WAVE FIG 1}.
%\begin{python}
%u=Vector(0.,ContinuousFunction(domain))
%\end{python}
%declares \var{u} as a vector-valued function of its coordinates 
%in the domain (i.e. with two or three components in the case of 
%a two or three dimensional domain, respectively). The first argument defines the value of the function which is equal
%to $0$ for all vectot components. The second argument \var{ContinuousFunction(domain)} 
%specifies the `smoothness` of the function. In our case we want the displacement field to be 
%a continuous function over the \Domain \var{domain}. On other option would be 
%\var{Function(domain)} in which case continuity is not guaranteed. In fact the 
%argument of the \function{grad} function has to be of the typ \var{ContinuousFunction(domain)} while
%returned function is of the var{Function(domain)} type. For more details see \Sec{SEC ESCRIPT DATA}. 

The statement \code{myPDE.setLumpingOn()}
switches on the use of an aggressive approximation of the PDE operator as a diagonal matrix
formed from the coefficient \var{D}.
The approximation allows, at the cost of 
additional error, very fast 
solution of the PDE. When using \code{setLumpingOn} the presence of \var{A}, \var{B} or \var
{C} will produce wrong results.


One of the big advantages of the Verlet scheme is the fact that the problem to be solved 
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).
The problem becomes so simple because we use the stress from the last time step rather then the stress which is
actually present at the current time step. Schemes using this approach are called an explicit time integration 
schemes \index{explicit scheme} \index{time integration!explicit}. The 
backward Euler scheme we have used in \Chap{DIFFUSION CHAP} is 
an example of an implicit scheme
\index{implicit scheme} \index{time integration!implicit}. In this case one uses the actual status of 
each variable at a particular time rather then values from previous time steps. This will lead to a problem which is 
more expensive to solve, in particular for non-linear problems. 
Although 
explicit time integration schemes are cheap to finalize a single time step, they need significantly smaller time
steps then implicit schemes and can suffer from stability problems. Therefore they need a 
very careful selection of the time step size $h$.

An easy, heuristic way of choosing an appropriate
time step size is the Courant condition \index{Courant condition} \index{explicit scheme!Courant condition}
which says that within a time step a information should not travel further than a cell used in the 
discretization scheme. In the case of the wave equation the velocity of a (p-) wave is given as
$\sqrt{\frac{\lambda+2\mu}{\rho}}$ so one should choose $h$ from
\begin{eqnarray}\label{WAVE dyn 66}
h= \frac{1}{5} \sqrt{\frac{\rho}{\lambda+2\mu}} \Delta x
\end{eqnarray}
where $\Delta x$ is the cell diameter. The factor $\frac{1}{5}$ is a safety factor considering the heuristics of 
the formula.

The following script uses the \code{wavePropagation} function to solve the
wave equation for a point source located at the bottom face of a block. The width of the block in 
each direction on the bottom face is $10\mbox{km}$ ($x\hackscore 0$ and $x\hackscore 1$ directions ie. \code{l0} and \code{l1}).
The \code{ne} gives the number of elements in $x\hackscore{0}$ and $x\hackscore 1$ directions.  
The depth of the block is aligned with the $x\hackscore{2}$-direction. 
The depth (\code{l2}) of the block in
the $x\hackscore{2}$-direction is chosen so that there are 10 elements and the magnitude of the
of the depth is chosen such that the elements 
become cubic. We chose 10 for the number of elements in $x\hackscore{2}$-direction so that the 
computation would be faster. \code{Brick($n\hackscore 0,  n\hackscore 1, n\hackscore 2,l\hackscore 0,l\hackscore 1,l\hackscore 2$)} is an \finley function which creates a rectangular mesh 
with $n\hackscore 0 \times n\hackscore 1 \times n\hackscore2$ elements over the brick $[0,l\hackscore 0] \times [0,l\hackscore 1] \times [0,l\hackscore 2]$.
\begin{python}
from esys.finley import Brick
ne=32          # number of cells in x_0 and x_1 directions
width=10000.  # length in x_0 and x_1 directions
lam=3.462e9
mu=3.462e9
rho=1154.
tend=60
h=(1./5.)*sqrt(rho/(lam+2*mu))*(width/ne)
print "time step size = ",h
U0=0.01 # amplitude of point source

mydomain=Brick(ne,ne,10,l0=width,l1=width,l2=10.*width/32.)
wavePropagation(mydomain,h,tend,lam,mu,rho,U0)
\end{python}
The script is available as \file{wave.py} in the \ExampleDirectory \index{scripts!\file{wave.py}}. Before running the code make sure you have created a directory called \code{data} in the current
working directory.
To visualize the results from the data directory: 
\begin{verbatim} mayavi -d usoln.1.vtu -m SurfaceMap &\end{verbatim} You can rotate this figure by clicking on it with the mouse and moving it around.
Again use \code{Configure Data} to move backwards and forward in time, and 
also to choose the results (acceleration, displacement or $u\hackscore x$) by using \code{Select Scalar}. Figure \ref{WAVE FIG 2} shows the results for the displacement at various time steps.

\begin{figure}[t!]
\centerline{\includegraphics[width=3.in,angle=270]{figures/WavePC.eps}}
\caption{Amplitude at Point source}
\label{WAVE FIG 1}
\end{figure}

\begin{figure}[t]
\begin{center}
\includegraphics[width=2in,angle=270]{figures/Wavet0.eps}
\includegraphics[width=2in,angle=270]{figures/Wavet1.eps}
\includegraphics[width=2in,angle=270]{figures/Wavet3.eps}
\includegraphics[width=2in,angle=270]{figures/Wavet10.eps}
\includegraphics[width=2in,angle=270]{figures/Wavet30.eps}
\includegraphics[width=2in,angle=270]{figures/Wavet288.eps}
\end{center}
\caption{Selected time steps ($n = 0,1,30,100,300$ and 2880) of a wave propagation over a $10\mbox{km} \times 10\mbox{km} \times 3.125\mbox{km}$ block 
from a point source initially at $(5\mbox{km},5\mbox{km},0)$
with time step size $h=0.02083$. Color represents the displacement.
Here the view is oriented onto the bottom face.
\label{WAVE FIG 2}}
\end{figure}
1	ksteube	1811
2			%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
3	ksteube	1316	%
4	ksteube	1811	% Copyright (c) 2003-2008 by University of Queensland
5			% Earth Systems Science Computational Center (ESSCC)
6			% http://www.uq.edu.au/esscc
7	gross	625	%
8	ksteube	1811	% Primary Business: Queensland, Australia
9			% Licensed under the Open Software License version 3.0
10			% http://www.opensource.org/licenses/osl-3.0.php
11	gross	625	%
12	ksteube	1811	%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
13	ksteube	1316
14	ksteube	1811
15	ksteube	1318	\section{3-D Wave Propagation}
16	jgs	107	\label{WAVE CHAP}
17
18	artak	1971	In this next example we want to calculate the displacement field $u\hackscore{i}$ for any time $t>0$ by solving the wave equation:
19			\index{wave equation}
20	jgs	107	\begin{eqnarray}\label{WAVE general problem}
21			\rho u\hackscore{i,tt} - \sigma\hackscore{ij,j}=0
22			\end{eqnarray}
23			in a three dimensional block of length $L$ in $x\hackscore{0}$
24			and $x\hackscore{1}$ direction and height $H$
25			in $x\hackscore{2}$ direction. $\rho$ is the known density which may be a function of its location.
26			$\sigma\hackscore{ij}$ is the stress field \index{stress} which in case of an isotropic, linear elastic material is given by
27			\begin{eqnarray} \label{WAVE stress}
28			\sigma\hackscore{ij} & = & \lambda u\hackscore{k,k} \delta\hackscore{ij} + \mu ( u\hackscore{i,j} + u\hackscore{j,i})
29			\end{eqnarray}
30			where $\lambda$ and $\mu$ are the Lame coefficients
31			\index{Lame coefficients} and $\delta\hackscore{ij}$ denotes the Kronecker symbol\index{Kronecker symbol}.
32			On the boundary the normal stress is given by
33			\begin{eqnarray} \label{WAVE natural}
34			\sigma\hackscore{ij}n\hackscore{j}=0
35			\end{eqnarray}
36	lkettle	581	for all time $t>0$.
37
38			At initial time $t=0$ the displacement
39			$u\hackscore{i}$ and the velocity $u\hackscore{i,t}$ are given:
40	jgs	107	\begin{eqnarray} \label{WAVE initial}
41	lkettle	581	u\hackscore{i}(0,x)= \left\{\begin{array}{cc} U\hackscore{0} & \mbox{for $x$ at point charge, $x\hackscore{C}$} \\ 0 & \mbox{elsewhere}\end{array} \right. & \mbox{ and } & u\hackscore{i,t}(0,x)=0 \;
42	jgs	107	\end{eqnarray}
43			for all $x$ in the domain.
44
45	lkettle	581	Here we are modelling a point source at the point $x\hackscore C$, in the numerical solution we
46			set the initial displacement to be $U\hackscore 0$ in a sphere of small radius around the point
47			$x\hackscore C$.
48	jgs	107
49	ksteube	1316	We use an explicit time integration scheme to calculate the displacement field $u$ at
50	jgs	107	certain time marks $t^{(n)}$ where $t^{(n)}=t^{(n-1)}+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$
51			which is defined by
52			\begin{eqnarray} \label{WAVE dyn 2}
53	lkettle	575	u^{(n)}=2u^{(n-1)}-u^{(n-2)} + h^2 a^{(n)} \\
54	jgs	107	\end{eqnarray}
55	lkettle	581	for all $n=2,3,\ldots$. It is designed to solve a system of equations of the form
56	jgs	107	\begin{eqnarray} \label{WAVE dyn 2b}
57			u\hackscore{,tt}=G(u)
58			\end{eqnarray}
59	gross	660	where one sets $a^{(n)}=G(u^{(n-1)})$.
60	jgs	107
61			In our case $a^{(n)}$ is given by
62			\begin{eqnarray}\label{WAVE dyn 3}
63			\rho a^{(n)}\hackscore{i}=\sigma^{(n-1)}\hackscore{ij,j}
64			\end{eqnarray}
65			and boundary conditions
66	gross	593	\begin{eqnarray} \label{WAVE natural at n}
67	jgs	107	\sigma^{(n-1)}\hackscore{ij}n\hackscore{j}=0
68			\end{eqnarray}
69			derived from \eqn{WAVE natural} where
70			\begin{eqnarray} \label{WAVE dyn 3a}
71	lkettle	581	\sigma\hackscore{ij}^{(n-1)} & = & \lambda u^{(n-1)}\hackscore{k,k} \delta\hackscore{ij} + \mu ( u^{(n-1)}\hackscore{i,j} + u^{(n-1)}\hackscore{j,i}).
72	jgs	107	\end{eqnarray}
73	lkettle	581	Now we have converted our problem for displacement, $u^{(n)}$, into a problem for
74			acceleration, $a^(n)$, which now depends
75			on the solution at the previous two time steps, $u^{(n-1)}$ and $u^{(n-2)}$.
76	jgs	107
77	ksteube	1316	In each time step we have to solve this problem to get the acceleration $a^{(n)}$, and we will
78	gross	565	use the \LinearPDE class of the \linearPDEs to do so. The general form of the PDE defined through
79	ksteube	1316	the \LinearPDE class is discussed in \Sec{SEC LinearPDE}. The form which is relevant here is
80	jgs	107	\begin{equation}\label{WAVE dyn 100}
81	lkettle	581	D\hackscore{ij} a^{(n)}\hackscore{j} = - X\hackscore{ij,j}\; .
82	jgs	107	\end{equation}
83	ksteube	1316	The natural boundary condition
84	jgs	107	\begin{equation}\label{WAVE dyn 101}
85			n\hackscore{j}X\hackscore{ij} =0
86			\end{equation}
87	ksteube	1316	is used.
88	jgs	107
89	ksteube	1316	With $u=a^{(n)}$ we can identify the values to be assigned to $D$ and $X$:
90	jgs	107	\begin{equation}\label{WAVE dyn 6}
91			\begin{array}{ll}
92			D\hackscore{ij}=\rho \delta\hackscore{ij}&
93			X\hackscore{ij}=-\sigma^{(n-1)}\hackscore{ij} \\
94			\end{array}
95			\end{equation}
96	lkettle	581
97
98			%Remember that \var{y.whereZero()} returns a function which is $1$ for those points $x$ where $y$ is zero
99			%(up to a certain tolerance, see \Sec{SEC ESCRIPT DATA}) and $0$ for those points where $y$ is not equal to zero.
100	jgs	107
101	lkettle	581	The following script defines a the function \function{wavePropagation} which
102	jgs	107	implements the Verlet scheme to solve our wave propagation problem.
103			The argument \var{domain} which is a \Domain class object
104	ksteube	1316	defines the domain of the problem. \var{h} and \var{tend} are the time step size
105			and the end time of the simulation. \var{lam}, \var{mu} and
106	lkettle	581	\var{rho} are material properties.
107	jgs	107	\begin{python}
108			from esys.linearPDEs import LinearPDE
109	lkettle	581	from numarray import identity,zeros,ones
110	jgs	107	from esys.escript import *
111	lkettle	581	from esys.escript.pdetools import Locator
112			def wavePropagation(domain,h,tend,lam,mu,rho,U0):
113	jgs	107	x=domain.getX()
114			# ... open new PDE ...
115			mypde=LinearPDE(domain)
116	lkettle	581	mypde.setSolverMethod(LinearPDE.LUMPING)
117	jgs	110	kronecker=identity(mypde.getDim())
118	lkettle	581	# spherical source at middle of bottom face
119			xc=[width/2.,width/2.,0.]
120			# define small radius around point xc
121			# Lsup(x) returns the maximum value of the argument x
122			src_radius = 0.1*Lsup(domain.getSize())
123			print "src_radius = ",src_radius
124			dunit=numarray.array([1.,0.,0.]) # defines direction of point source
125			mypde.setValue(D=kronecker*rho)
126	jgs	107	# ... set initial values ....
127	jgs	110	n=0
128	lkettle	581	# initial value of displacement at point source is constant (U0=0.01)
129			# for first two time steps
130			u=U0whereNegative(length(x-xc)-src_radius)dunit
131			u_last=U0whereNegative(length(x-xc)-src_radius)dunit
132	jgs	107	t=0
133	lkettle	581	# define the location of the point source
134			L=Locator(domain,xc)
135			# find potential at point source
136			u_pc=L.getValue(u)
137			print "u at point charge=",u_pc
138			u_pc_x = u_pc[0]
139			u_pc_y = u_pc[1]
140			u_pc_z = u_pc[2]
141			# open file to save displacement at point source
142			u_pc_data=open('./data/U_pc.out','w')
143			u_pc_data.write("%f %f %f %f\n"%(t,u_pc_x,u_pc_y,u_pc_z))
144	jgs	107	while t<tend:
145			# ... get current stress ....
146			g=grad(u)
147	jgs	110	stress=lamtrace(g)kronecker+mu*(g+transpose(g))
148			# ... get new acceleration ....
149	lkettle	581	mypde.setValue(X=-stress)
150	jgs	107	a=mypde.getSolution()
151			# ... get new displacement ...
152			u_new=2u-u_last+h2a
153			# ... shift displacements ....
154			u_last=u
155			u=u_new
156			t+=h
157	jgs	110	n+=1
158			print n,"-th time step t ",t
159	lkettle	581	L=Locator(domain,xc)
160			u_pc=L.getValue(u)
161			print "u at point charge=",u_pc
162			u_pc_x=u_pc[0]
163			u_pc_y=u_pc[1]
164			u_pc_z=u_pc[2]
165			# save displacements at point source to file for t > 0
166			u_pc_data.write("%f %f %f %f\n"%(t,u_pc_x,u_pc_y,u_pc_z))
167			# ... save current acceleration in units of gravity and displacements
168			if n==1 or n%10==0: saveVTK("./data/usoln.%i.vtu"%(n/10),acceleration=length(a)/9.81,
169			displacement = length(u), Ux = u[0] )
170			u_pc_data.close()
171	jgs	107	\end{python}
172	jgs	110	Notice that
173	ksteube	1316	all coefficients of the PDE which are independent of time $t$ are set outside the \code{while}
174			loop. This is very efficient since it allows the \LinearPDE class to reuse information as much as possible
175	jgs	110	when iterating over time.
176
177	ksteube	1316	there are a few new \escript functions in this example:
178	jgs	107	\code{grad(u)} returns the gradient $u\hackscore{i,j}$ of $u$ (in fact \var{grad(g)[i,j]}=$=u\hackscore{i,j}$).
179	ksteube	1316	There are restrictions on the argument of the \function{grad} function, for instance
180			the statement \code{grad(grad(u))} will raise an exception.
181	jgs	107	\code{trace(g)} returns the sum of the main diagonal elements \var{g[k,k]} of \var{g}
182	ksteube	1316	and \code{transpose(g)} returns the matrix transpose of \var{g} (ie. $\var{transpose(g)[i,j]}=\var{g[j,i]}$).
183	jgs	107
184	lkettle	581	We initialize the values of \code{u} and \code{u_last} to be $U\hackscore 0$ for a small
185			sphere of radius, \code{src_radius} around the point source located at $x\hackscore C$.
186			These satisfy the initial conditions defined in \eqn{WAVE initial}.
187	jgs	107
188	lkettle	581	The output \code{U_pc.out} and \code{vtu} files are saved in a directory called \code{data}.
189			The \code{U_pc.out} stores 4 columns of data: $t,u\hackscore x,u\hackscore y,u\hackscore z$
190			respectively, where $u\hackscore x,u\hackscore y,u\hackscore z$ are the $x\hackscore 0,x\hackscore 1,x\hackscore 2$ components of
191			the displacement vector $u$ at the point source. These can be
192			plotted easily using any plotting package. In gnuplot the command:
193			\code{plot 'U_pc.out' u 1:2 title 'U_x' w l lw 2, 'U_pc.out' u 1:3 title 'U_y' w l lw 2,
194			'U_pc.out' u 1:4 title 'U_z' w l lw 2} will reproduce Figure \ref{WAVE FIG 1}.
195			%\begin{python}
196			%u=Vector(0.,ContinuousFunction(domain))
197			%\end{python}
198			%declares \var{u} as a vector-valued function of its coordinates
199			%in the domain (i.e. with two or three components in the case of
200	ksteube	1316	%a two or three dimensional domain, respectively). The first argument defines the value of the function which is equal
201	lkettle	581	%to $0$ for all vectot components. The second argument \var{ContinuousFunction(domain)}
202			%specifies the `smoothness` of the function. In our case we want the displacement field to be
203			%a continuous function over the \Domain \var{domain}. On other option would be
204	ksteube	1316	%\var{Function(domain)} in which case continuity is not guaranteed. In fact the
205	lkettle	581	%argument of the \function{grad} function has to be of the typ \var{ContinuousFunction(domain)} while
206			%returned function is of the var{Function(domain)} type. For more details see \Sec{SEC ESCRIPT DATA}.
207
208	jgs	110	The statement \code{myPDE.setLumpingOn()}
209	ksteube	1316	switches on the use of an aggressive approximation of the PDE operator as a diagonal matrix
210			formed from the coefficient \var{D}.
211			The approximation allows, at the cost of
212			additional error, very fast
213			solution of the PDE. When using \code{setLumpingOn} the presence of \var{A}, \var{B} or \var
214			{C} will produce wrong results.
215	jgs	110
216	jgs	107
217	jgs	110	One of the big advantages of the Verlet scheme is the fact that the problem to be solved
218	ksteube	1316	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).
219	lkettle	581	The problem becomes so simple because we use the stress from the last time step rather then the stress which is
220			actually present at the current time step. Schemes using this approach are called an explicit time integration
221	ksteube	1316	schemes \index{explicit scheme} \index{time integration!explicit}. The
222	jgs	110	backward Euler scheme we have used in \Chap{DIFFUSION CHAP} is
223	ksteube	1316	an example of an implicit scheme
224	jgs	110	\index{implicit scheme} \index{time integration!implicit}. In this case one uses the actual status of
225	lkettle	581	each variable at a particular time rather then values from previous time steps. This will lead to a problem which is
226	jgs	110	more expensive to solve, in particular for non-linear problems.
227			Although
228	lkettle	581	explicit time integration schemes are cheap to finalize a single time step, they need significantly smaller time
229			steps then implicit schemes and can suffer from stability problems. Therefore they need a
230	jgs	110	very careful selection of the time step size $h$.
231	jgs	107
232	jgs	110	An easy, heuristic way of choosing an appropriate
233	jgs	107	time step size is the Courant condition \index{Courant condition} \index{explicit scheme!Courant condition}
234			which says that within a time step a information should not travel further than a cell used in the
235			discretization scheme. In the case of the wave equation the velocity of a (p-) wave is given as
236			$\sqrt{\frac{\lambda+2\mu}{\rho}}$ so one should choose $h$ from
237			\begin{eqnarray}\label{WAVE dyn 66}
238			h= \frac{1}{5} \sqrt{\frac{\rho}{\lambda+2\mu}} \Delta x
239			\end{eqnarray}
240	lkettle	581	where $\Delta x$ is the cell diameter. The factor $\frac{1}{5}$ is a safety factor considering the heuristics of
241	jgs	107	the formula.
242
243	lkettle	581	The following script uses the \code{wavePropagation} function to solve the
244			wave equation for a point source located at the bottom face of a block. The width of the block in
245			each direction on the bottom face is $10\mbox{km}$ ($x\hackscore 0$ and $x\hackscore 1$ directions ie. \code{l0} and \code{l1}).
246	artak	1971	The \code{ne} gives the number of elements in $x\hackscore{0}$ and $x\hackscore 1$ directions.
247	lkettle	581	The depth of the block is aligned with the $x\hackscore{2}$-direction.
248			The depth (\code{l2}) of the block in
249			the $x\hackscore{2}$-direction is chosen so that there are 10 elements and the magnitude of the
250			of the depth is chosen such that the elements
251			become cubic. We chose 10 for the number of elements in $x\hackscore{2}$-direction so that the
252			computation would be faster. \code{Brick($n\hackscore 0, n\hackscore 1, n\hackscore 2,l\hackscore 0,l\hackscore 1,l\hackscore 2$)} is an \finley function which creates a rectangular mesh
253			with $n\hackscore 0 \times n\hackscore 1 \times n\hackscore2$ elements over the brick $[0,l\hackscore 0] \times [0,l\hackscore 1] \times [0,l\hackscore 2]$.
254	jgs	107	\begin{python}
255	lkettle	581	from esys.finley import Brick
256			ne=32 # number of cells in x_0 and x_1 directions
257			width=10000. # length in x_0 and x_1 directions
258	jgs	110	lam=3.462e9
259			mu=3.462e9
260	jgs	107	rho=1154.
261	jgs	110	tend=60
262	lkettle	581	h=(1./5.)sqrt(rho/(lam+2mu))*(width/ne)
263			print "time step size = ",h
264			U0=0.01 # amplitude of point source
265	jgs	110
266	lkettle	581	mydomain=Brick(ne,ne,10,l0=width,l1=width,l2=10.*width/32.)
267			wavePropagation(mydomain,h,tend,lam,mu,rho,U0)
268	jgs	107	\end{python}
269	lkettle	581	The script is available as \file{wave.py} in the \ExampleDirectory \index{scripts!\file{wave.py}}. Before running the code make sure you have created a directory called \code{data} in the current
270			working directory.
271			To visualize the results from the data directory:
272			\begin{verbatim} mayavi -d usoln.1.vtu -m SurfaceMap &\end{verbatim} You can rotate this figure by clicking on it with the mouse and moving it around.
273			Again use \code{Configure Data} to move backwards and forward in time, and
274			also to choose the results (acceleration, displacement or $u\hackscore x$) by using \code{Select Scalar}. Figure \ref{WAVE FIG 2} shows the results for the displacement at various time steps.
275	jgs	107
276	jfenwick	2104	\begin{figure}[t!]
277	gross	599	\centerline{\includegraphics[width=3.in,angle=270]{figures/WavePC.eps}}
278	gross	593	\caption{Amplitude at Point source}
279			\label{WAVE FIG 1}
280			\end{figure}
281	jgs	107
282	jfenwick	2104	\begin{figure}[t]
283	lkettle	581	\begin{center}
284	gross	599	\includegraphics[width=2in,angle=270]{figures/Wavet0.eps}
285			\includegraphics[width=2in,angle=270]{figures/Wavet1.eps}
286			\includegraphics[width=2in,angle=270]{figures/Wavet3.eps}
287			\includegraphics[width=2in,angle=270]{figures/Wavet10.eps}
288			\includegraphics[width=2in,angle=270]{figures/Wavet30.eps}
289			\includegraphics[width=2in,angle=270]{figures/Wavet288.eps}
290	lkettle	581	\end{center}
291			\caption{Selected time steps ($n = 0,1,30,100,300$ and 2880) of a wave propagation over a $10\mbox{km} \times 10\mbox{km} \times 3.125\mbox{km}$ block
292			from a point source initially at $(5\mbox{km},5\mbox{km},0)$
293			with time step size $h=0.02083$. Color represents the displacement.
294	gross	707	Here the view is oriented onto the bottom face.
295			\label{WAVE FIG 2}}
296	jgs	110	\end{figure}
Name	Value
svn:eol-style	native
svn:keywords	Author Date Id Revision