# Diff of /trunk/doc/user/wave.tex

revision 580 by lkettle, Fri Mar 3 03:33:07 2006 UTC revision 581 by lkettle, Wed Mar 8 05:48:51 2006 UTC
# Line 20  On the boundary the normal stress is giv Line 20  On the boundary the normal stress is giv
20  \begin{eqnarray} \label{WAVE natural}  \begin{eqnarray} \label{WAVE natural}
21  \sigma\hackscore{ij}n\hackscore{j}=0  \sigma\hackscore{ij}n\hackscore{j}=0
22  \end{eqnarray}  \end{eqnarray}
23  for all time $t>0$. At initial time $t=0$ the displacement  for all time $t>0$.
24  $u\hackscore{i}$ and the velocity $u\hackscore{i,t}$ are give:
25    At initial time $t=0$ the displacement
26    $u\hackscore{i}$ and the velocity $u\hackscore{i,t}$ are given:
27  \begin{eqnarray} \label{WAVE initial}  \begin{eqnarray} \label{WAVE initial}
28  u\hackscore{i}(0,x)=0 & \mbox{ and } & u\hackscore{i,t}(0,x)=0 \;  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 \;
29  \end{eqnarray}  \end{eqnarray}
30  for all $x$ in the domain.  for all $x$ in the domain.
31
32  The points on the left face of the domain, ie. points $x=(x\hackscore{i})$ with $x\hackscore{0}=0$,  Here we are modelling a point source at the point $x\hackscore C$, in the numerical solution we
33  are moved into the $x\hackscore{2}$ direction by $u^{max}$. This is expressed by the constraint   set the initial displacement to be $U\hackscore 0$ in a sphere of small radius around the point
34  \begin{eqnarray} \label{WAVE impact}   $x\hackscore C$.
u\hackscore{2}(x,t)=s(t) \mbox{ for } x=(x\hackscore{i}) with  x\hackscore{0}=0
\end{eqnarray}
where
\begin{eqnarray} \label{WAVE impact s}
s(t)=u^{max} (1-e^{-(\tau^{-1}t)^3})
\end{eqnarray}
$\tau$ measures the time needed to reach the maximum displacement $u^{max}$ (
actually $\tau (ln 2)^{\frac{1}{3}}) \approx 0.9 \tau$ is the time until $\frac{u^{max}}{2}$ is reached).
The smaller $\tau$ the faster $u^{max}$ is reached.
35
36  We use an explicit time-integration scheme to calculate the displacement field $u$ at  We use an explicit time-integration scheme to calculate the displacement field $u$ at
37  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$  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$
# Line 46  which is defined by Line 39  which is defined by
39  \begin{eqnarray} \label{WAVE dyn 2}  \begin{eqnarray} \label{WAVE dyn 2}
40  u^{(n)}=2u^{(n-1)}-u^{(n-2)} + h^2 a^{(n)} \\  u^{(n)}=2u^{(n-1)}-u^{(n-2)} + h^2 a^{(n)} \\
41  \end{eqnarray}  \end{eqnarray}
42  for all $n=1,2,\ldots$. It is designed to solve a system of equations of the form  for all $n=2,3,\ldots$. It is designed to solve a system of equations of the form
43  \begin{eqnarray} \label{WAVE dyn 2b}  \begin{eqnarray} \label{WAVE dyn 2b}
44  u\hackscore{,tt}=G(u)  u\hackscore{,tt}=G(u)
45  \end{eqnarray}  \end{eqnarray}
# Line 62  and boundary conditions Line 55  and boundary conditions
55  \end{eqnarray}  \end{eqnarray}
56  derived from \eqn{WAVE natural} where  derived from \eqn{WAVE natural} where
57  \begin{eqnarray} \label{WAVE dyn 3a}  \begin{eqnarray} \label{WAVE dyn 3a}
58  \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})  \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}.
Additionally $a^{(n)}$ has to fullfill the constraint
\begin{eqnarray}\label{WAVE dyn 4}
a^{(n)}\hackscore{2}(x)= s\hackscore{,tt}(t^{(n)})\mbox{ where } s\hackscore{,tt}(t)=\frac{u^{max}}{\tau^2}
(6\tau^{-1}t- 9(\tau^{-1}t)^4) e^{-(\tau^{-1}t)^3}
59  \end{eqnarray}  \end{eqnarray}
60  derived from \eqn{WAVE impact}.  Now we have converted our problem for displacement, $u^{(n)}$, into a problem for
61    acceleration, $a^(n)$, which now depends
62    on the solution at the previous two time steps, $u^{(n-1)}$  and $u^{(n-2)}$.
63
64  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
65  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
66  the \LinearPDE class is discussed in \Sec{SEC LinearPDE}. The form which are relevant here are    the \LinearPDE class is discussed in \Sec{SEC LinearPDE}. The forms which are relevant here are
67  \label{WAVE dyn 100}  \label{WAVE dyn 100}
68  D\hackscore{ij} u\hackscore{j} = - X\hackscore{ij,j}\; .  D\hackscore{ij} a^{(n)}\hackscore{j} = - X\hackscore{ij,j}\; .
69
70  Implicitly, the natural boundary condition  Implicitly, the natural boundary condition
71  \label{WAVE dyn 101}  \label{WAVE dyn 101}
# Line 86  is assumed. The \LinearPDE object allows Line 75  is assumed. The \LinearPDE object allows
75  \label{WAVE dyn 101}  \label{WAVE dyn 101}
76  u\hackscore{i}(x)=r\hackscore{i}(x) \mbox{ for } x \mbox{ with } q\hackscore{i}>0  u\hackscore{i}(x)=r\hackscore{i}(x) \mbox{ for } x \mbox{ with } q\hackscore{i}>0
77
78  where $r$ and $q$ are given functions.  where $r$ and $q$ are given functions. However in this problem we don't need to define
79    any constraints for our problem.
80
81  With $u=a^{(n)}$ we can identify the values to be assigned to $D$, $X$, $r$ and $q$:  With $u=a^{(n)}$ we can identify the values to be assigned to $D$, $X$:
82  \label{WAVE  dyn 6}  \label{WAVE  dyn 6}
83  \begin{array}{ll}  \begin{array}{ll}
84  D\hackscore{ij}=\rho \delta\hackscore{ij}&  D\hackscore{ij}=\rho \delta\hackscore{ij}&
85  X\hackscore{ij}=-\sigma^{(n-1)}\hackscore{ij} \\  X\hackscore{ij}=-\sigma^{(n-1)}\hackscore{ij} \\
r\hackscore{i}=\delta\hackscore{i2} \cdot s\hackscore{,tt}(t^{(n)}) & q\hackscore{i}=\delta\hackscore{i2}  \cdot x\hackscore{0}.\texttt{whereZero()}
86  \end{array}  \end{array}
87
88  Remember that \var{y.whereZero()} returns a function which is $1$ for those points $x$ where $y$ is zero
89  (up to a certain tolerance, see \Sec{SEC ESCRIPT DATA}) and $0$ for those points where $y$ is not equal to zero.
90    %Remember that \var{y.whereZero()} returns a function which is $1$ for those points $x$ where $y$ is zero
91    %(up to a certain tolerance, see \Sec{SEC ESCRIPT DATA}) and $0$ for those points where $y$ is not equal to zero.
92
93  In the following script defines a the function \function{wavePropagation} which  The following script defines a the function \function{wavePropagation} which
94  implements the Verlet scheme to solve our wave propagation problem.  implements the Verlet scheme to solve our wave propagation problem.
95  The argument \var{domain} which is a \Domain class object  The argument \var{domain} which is a \Domain class object
96  defines the domain of the problem. \var{h} and \var{tend} is the step size  defines the domain of the problem. \var{h} and \var{tend} are the step size
97  and the time mark after which the simulation will be terminated. \var{lam}, \var{mu} and  and the time mark after which the simulation will be terminated respectively. \var{lam}, \var{mu} and
98  \var{rho} are material properties. \var{s_tt} is a function which returns $s\hackscore{,tt}$ at a given time:  \var{rho} are material properties.
99  \begin{python}  \begin{python}
100  from esys.linearPDEs import LinearPDE  from esys.linearPDEs import LinearPDE
101  from numarray import identity  from numarray import identity,zeros,ones
102  from esys.escript import *  from esys.escript import *
103  def wavePropagation(domain,h,tend,lam,mu,rho,s_tt):  from esys.escript.pdetools import Locator
104    def wavePropagation(domain,h,tend,lam,mu,rho,U0):
105     x=domain.getX()     x=domain.getX()
106     # ... open new PDE ...     # ... open new PDE ...
107     mypde=LinearPDE(domain)     mypde=LinearPDE(domain)
108     mypde.setLumpingOn()     mypde.setSolverMethod(LinearPDE.LUMPING)
109     kronecker=identity(mypde.getDim())     kronecker=identity(mypde.getDim())
110     mypde.setValue(D=kronecker*rho, \     #  spherical source at middle of bottom face
111                    q=x[0].whereZero()*kronecker[1,:])     xc=[width/2.,width/2.,0.]
112       # define small radius around point xc
113       # Lsup(x) returns the maximum value of the argument x
116       dunit=numarray.array([1.,0.,0.]) # defines direction of point source
117       mypde.setValue(D=kronecker*rho)
118     # ... set initial values ....     # ... set initial values ....
119     n=0     n=0
120     u=Vector(0,ContinuousFunction(domain))     # initial value of displacement at point source is constant (U0=0.01)
121     u_last=Vector(0,ContinuousFunction(domain))     # for first two time steps
124     t=0     t=0
125       # define the location of the point source
126       L=Locator(domain,xc)
127       # find potential at point source
128       u_pc=L.getValue(u)
129       print "u at point charge=",u_pc
130       u_pc_x = u_pc[0]
131       u_pc_y = u_pc[1]
132       u_pc_z = u_pc[2]
133       # open file to save displacement at point source
134       u_pc_data=open('./data/U_pc.out','w')
135       u_pc_data.write("%f %f %f %f\n"%(t,u_pc_x,u_pc_y,u_pc_z))
136     while t<tend:     while t<tend:
137       # ... get current stress ....       # ... get current stress ....
139       stress=lam*trace(g)*kronecker+mu*(g+transpose(g))       stress=lam*trace(g)*kronecker+mu*(g+transpose(g))
140       # ... get new acceleration ....       # ... get new acceleration ....
141       mypde.setValue(X=-stress,r=s_tt(t+h)*kronecker[1,:])       mypde.setValue(X=-stress)
142       a=mypde.getSolution()       a=mypde.getSolution()
143       # ... get new displacement ...       # ... get new displacement ...
144       u_new=2*u-u_last+h**2*a       u_new=2*u-u_last+h**2*a
# Line 137  def wavePropagation(domain,h,tend,lam,mu Line 148  def wavePropagation(domain,h,tend,lam,mu
148       t+=h       t+=h
149       n+=1       n+=1
150       print n,"-th time step t ",t       print n,"-th time step t ",t
151       print "a=",inf(a),sup(a)       L=Locator(domain,xc)
152       print "u=",inf(u),sup(u)       u_pc=L.getValue(u)
153       # ... save current acceleration in units of gravity       print "u at point charge=",u_pc
154       if n%10==0 : (length(a)/9.81).saveDX("/tmp/res/u.%i.dx"%(n/10))       u_pc_x=u_pc[0]
155         u_pc_y=u_pc[1]
156         u_pc_z=u_pc[2]
157         # save displacements at point source to file for t > 0
158         u_pc_data.write("%f %f %f %f\n"%(t,u_pc_x,u_pc_y,u_pc_z))
159         # ... save current acceleration in units of gravity and displacements
160         if n==1 or n%10==0: saveVTK("./data/usoln.%i.vtu"%(n/10),acceleration=length(a)/9.81,
161         displacement = length(u), Ux = u[0] )
162       u_pc_data.close()
163  \end{python}  \end{python}
164  Notice that  Notice that
165  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 from time $t$ are set outside the \code{while}
166  loop. This allows the \LinearPDE class to ruse information as much as possible  loop. This allows the \LinearPDE class to reuse information as much as possible
167  when iterating over time.    when iterating over time.
168
169  We have seen most of the functions in previous examples but there are some new functions here:  We have seen most of the functions in previous examples but there are some new functions here:
# Line 155  for \finley the statement \code{grad(gra Line 174  for \finley the statement \code{grad(gra
174  and \code{transpose(g)} returns the transposed of the matrix \var{g} (ie.  and \code{transpose(g)} returns the transposed of the matrix \var{g} (ie.
175  $\var{transpose(g)[i,j]}=\var{g[j,i]}$).  $\var{transpose(g)[i,j]}=\var{g[j,i]}$).
176
177  The initialization of \var{u} and \var{u_last} needs some explanation. The statement  We initialize the values of \code{u} and \code{u_last} to be $U\hackscore 0$ for a small
178  \begin{python}  sphere of radius, \code{src_radius} around the point source located at $x\hackscore C$.
179  u=Vector(0.,ContinuousFunction(domain))  These satisfy the initial conditions defined in \eqn{WAVE initial}.
180  \end{python}
181  declares \var{u} as a vector-valued function of its coordinates  The output \code{U_pc.out} and \code{vtu} files are saved in a directory called \code{data}.
182  in the domain (i.e. with two or three components in the case of  The \code{U_pc.out} stores 4 columns of data: $t,u\hackscore x,u\hackscore y,u\hackscore z$
183  a two or three dimensional domain, repectively). The first argument defines the value of the function which is equal  respectively, where $u\hackscore x,u\hackscore y,u\hackscore z$ are the $x\hackscore 0,x\hackscore 1,x\hackscore 2$ components of
184  to $0$ for all vectot components. The second argument \var{ContinuousFunction(domain)}  the displacement vector $u$ at the point source. These can be
185  specifies the smoothness of the function. In our case we want the displacement field to be  plotted easily using any plotting package. In gnuplot the command:
186  a continuous function over the \Domain \var{domain}. On other option would be  \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,
187  \var{Function(domain)} in which case continuity is not garanteed. In fact the  'U_pc.out' u 1:4 title 'U_z' w l lw 2} will reproduce Figure \ref{WAVE FIG 1}.
188  argument of the \function{grad} function has to be of the typ \var{ContinuousFunction(domain)} while  %\begin{python}
189  returned function is of the var{Function(domain)} type. For more details see \Sec{SEC ESCRIPT DATA}.  %u=Vector(0.,ContinuousFunction(domain))
190    %\end{python}
191    %declares \var{u} as a vector-valued function of its coordinates
192    %in the domain (i.e. with two or three components in the case of
193    %a two or three dimensional domain, repectively). The first argument defines the value of the function which is equal
194    %to $0$ for all vectot components. The second argument \var{ContinuousFunction(domain)}
195    %specifies the smoothness of the function. In our case we want the displacement field to be
196    %a continuous function over the \Domain \var{domain}. On other option would be
197    %\var{Function(domain)} in which case continuity is not garanteed. In fact the
198    %argument of the \function{grad} function has to be of the typ \var{ContinuousFunction(domain)} while
199    %returned function is of the var{Function(domain)} type. For more details see \Sec{SEC ESCRIPT DATA}.
200
201  The statement \code{myPDE.setLumpingOn()}  The statement \code{myPDE.setLumpingOn()}
202  switches on the usage of an aggressive approximation of the PDE operator, in this case  switches on the usage of an aggressive approximation of the PDE operator, in this case
203  formed by the coefficient \var{D} (actually the discrete  formed by the coefficient \var{D} (actually the discrete
204  version of the PDE operator is approximated by a diagonal matrix). The approximation allows, at the costs of  version of the PDE operator is approximated by a diagonal matrix). The approximation allows, at the costs of
206  solving of the PDE when the solution is requested. In the general case of the presents of \var{A}, \var{B} or \var  solving of the PDE when the solution is requested. In the general case, the presence of \var{A}, \var{B} or \var
207  {C} \code{setLumpingOn} will produce wrong results.  {C} \code{setLumpingOn} will produce wrong results.
208
It is pointed out that the value of the constraint \var{r} for the acceleration \var{a} at time $t=0$ is consistent
with the initial value of the acceleration which is identical zero. Otherwise, a discontinuity of \var{a} in time
direction would lead to heavy oscillations in the displacement. It is also pointed out, that at time $t=0$ the
constraint defined by \eqn{WAVE impact} fulfills
the initial conditions in \eqn{WAVE initial}.

209
210  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
211  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 actually allows us to use lumping).
212  The problem becomes so simple because we use the stress from the last time step rather then the stress which  The problem becomes so simple because we use the stress from the last time step rather then the stress which is
213  actually is 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
214  scheme \index{explicit scheme} \index{time integration!explicit}. The  scheme \index{explicit scheme} \index{time integration!explicit}. The
215  backward Euler scheme we have used in \Chap{DIFFUSION CHAP} is  backward Euler scheme we have used in \Chap{DIFFUSION CHAP} is
216  an example of an implicit schemes  an example of an implicit schemes
217  \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
218  each variable at a particular time rather then values from previous time steps. This will lead to problem which is  each variable at a particular time rather then values from previous time steps. This will lead to a problem which is
219  more expensive to solve, in particular for non-linear problems.  more expensive to solve, in particular for non-linear problems.
220  Although  Although
221  explicit time integration schemes are cheep to finalize a single time step, they need a significant smaller time  explicit time integration schemes are cheap to finalize a single time step, they need significantly smaller time
222  step then implicit schemes and can suffer from stability problems. Therefor they need a  steps then implicit schemes and can suffer from stability problems. Therefore they need a
223  very careful selection of the time step size $h$.  very careful selection of the time step size $h$.
224
225  An easy, heuristic way of choosing an appropriate  An easy, heuristic way of choosing an appropriate
# Line 207  $\sqrt{\frac{\lambda+2\mu}{\rho}}$ so on Line 230  $\sqrt{\frac{\lambda+2\mu}{\rho}}$ so on
230  \begin{eqnarray}\label{WAVE dyn 66}  \begin{eqnarray}\label{WAVE dyn 66}
231  h= \frac{1}{5} \sqrt{\frac{\rho}{\lambda+2\mu}} \Delta x  h= \frac{1}{5} \sqrt{\frac{\rho}{\lambda+2\mu}} \Delta x
232  \end{eqnarray}  \end{eqnarray}
233  where $\Delta x$ is the cell diameter. The factor $\frac{1}{5}$ is a safety factor considering the heuristic of  where $\Delta x$ is the cell diameter. The factor $\frac{1}{5}$ is a safety factor considering the heuristics of
234  the formula.  the formula.
235
236  The following script uses the \code{wavePropagtion} function to solve the  The following script uses the \code{wavePropagation} function to solve the
237  wave equation over a block in crust of the earth. The depth of the block is  wave equation for a point source located at the bottom face of a block. The width of the block in
238  $10km$ and the width in each direction is $100km$. The depth of the block is alligned  each direction on the bottom face is $10\mbox{km}$ ($x\hackscore 0$ and $x\hackscore 1$ directions ie. \code{l0} and \code{l1}).
239  with the $x\hackscore{0}$-direction. \var{ne} gives the number of elements in $x\hackscore{0}$-direction. The number  \code{ne} gives the number of elements in $x\hackscore{0}$ and $x\hackscore 1$ directions.
240  of elements in $x\hackscore{1}$- and $x\hackscore{2}$-direction is chosen such that the elements  The depth of the block is aligned with the $x\hackscore{2}$-direction.
241  become cubic. The impact is $2m$ acting within a time of about $10sec$:  The depth (\code{l2}) of the block in
242    the $x\hackscore{2}$-direction is chosen so that there are 10 elements and the magnitude of the
243    of the depth is chosen such that the elements
244    become cubic. We chose 10 for the number of elements in $x\hackscore{2}$-direction so that the
245    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
246    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]$.
247  \begin{python}  \begin{python}
248  ne=10           # number of cells in x_0-direction  from esys.finley import Brick
249  depth=10000.   # length in x_0-direction  ne=32          # number of cells in x_0 and x_1 directions
250  width=100000.  # length in x_1 and x_2 direction  width=10000.  # length in x_0 and x_1 directions
251  lam=3.462e9  lam=3.462e9
252  mu=3.462e9  mu=3.462e9
253  rho=1154.  rho=1154.
tau=10.
umax=2.
254  tend=60  tend=60
255  h=1./5.*sqrt(rho/(lam+2*mu))*(depth/ne)  h=(1./5.)*sqrt(rho/(lam+2*mu))*(width/ne)
def s_tt(t): return umax/tau**2*(6*t/tau-9*(t/tau)**4)*exp(-(t/tau)**3)

256  print "time step size = ",h  print "time step size = ",h
257  mydomain=Brick(ne,int(width/depth)*ne,int(width/depth)*ne,l0=depth,l1=width,l2=width)  U0=0.01 # amplitude of point source
wavePropagation(mydomain,h,tend,lam,mu,rho,s_tt)
\end{python}
The script is available as \file{wave.py} in the \ExampleDirectory \index{scripts!\file{wave.py}}.
258
259    mydomain=Brick(ne,ne,10,l0=width,l1=width,l2=10.*width/32.)
260  % \begin{figure}  wavePropagation(mydomain,h,tend,lam,mu,rho,U0)
261  % \centerline{\includegraphics[width=\figwidth]{DiffusionDomain}}  \end{python}
262  % \caption{Temperature Diffusion Problem with Circular Heat Source}  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
263  % \label{WAVE FIG 1}  working directory.
264  % \end{figure}  To visualize the results from the data directory:
265    \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.
266  \begin{figure}  Again use \code{Configure Data} to move backwards and forward in time, and
267  % \centerline{\includegraphics[width=\figwidth]{DiffusionRes1}}  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.
268  % \centerline{\includegraphics[width=\figwidth]{DiffusionRes16}}
269  % \centerline{\includegraphics[width=\figwidth]{DiffusionRes32}}   \begin{figure}{t!}
270  %\centerline{\includegraphics[width=\figwidth]{DiffusionRes48}}   \centerline{\includegraphics[width=3.in,angle=270]{WavePC.ps}}
271  \caption{Selected time steps of a wave propagation over a $10km \times 100km \times 100km$   \caption{Amplitude at Point source}
272  with time step size $h=0.0666$. Color represents the acceleration.}   \label{WAVE FIG 1}
273     \end{figure}
274
275    \begin{figure}{t}
276    \begin{center}
277    \includegraphics[width=2in,angle=270]{Wavet0.ps}
278    \includegraphics[width=2in,angle=270]{Wavet1.ps}
279    \includegraphics[width=2in,angle=270]{Wavet3.ps}
280    \includegraphics[width=2in,angle=270]{Wavet10.ps}
281    \includegraphics[width=2in,angle=270]{Wavet30.ps}
282    \includegraphics[width=2in,angle=270]{Wavet288.ps}
283    \end{center}
284    \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
285    from a point source initially at $(5\mbox{km},5\mbox{km},0)$
286    with time step size $h=0.02083$. Color represents the displacement.
287    Here the view is oriented onto the bottom face.}
288  \label{WAVE FIG 2}  \label{WAVE FIG 2}
289  \end{figure}  \end{figure}

Legend:
 Removed from v.580 changed lines Added in v.581