/[escript]/trunk/doc/user/wave.tex
ViewVC logotype

Contents of /trunk/doc/user/wave.tex

Parent Directory Parent Directory | Revision Log Revision Log


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

Properties

Name Value
svn:eol-style native
svn:keywords Author Date Id Revision

  ViewVC Help
Powered by ViewVC 1.1.26