ViewVC logotype

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

Parent Directory Parent Directory | Revision Log Revision Log

Revision 1318 - (show annotations)
Wed Sep 26 04:39:14 2007 UTC (13 years, 3 months ago) by ksteube
File MIME type: application/x-tex
File size: 14991 byte(s)
A very quick edit of chapter 3 of the User Guide. More editing needed.

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


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

  ViewVC Help
Powered by ViewVC 1.1.26