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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 2503 - (show annotations)
Wed Jul 1 05:05:08 2009 UTC (9 years, 9 months ago) by jfenwick
File MIME type: application/x-tex
File size: 15859 byte(s)
Fixing doc build so examples tar and zip are built properly.
Fixed image filenames in user guide.

Please do not put extensions on image files.


1
2 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
3 %
4 % Copyright (c) 2003-2008 by University of Queensland
5 % Earth Systems Science Computational Center (ESSCC)
6 % http://www.uq.edu.au/esscc
7 %
8 % 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 %
12 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
13
14
15 \section{3-D Wave Propagation}
16 \label{WAVE CHAP}
17
18 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 \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 for all time $t>0$.
37
38 \begin{figure}[t!]
39 \centerline{\includegraphics[angle=-90,width=4.in]{figures/waveu}}
40 \caption{Displacement at Source Point}
41 \label{WAVE FIG 1.2}
42 \end{figure}
43
44 \begin{figure}[t!]
45 \centerline{\includegraphics[angle=-90,width=4.in]{figures/wavea}}
46 \caption{Acceleration at Source Point}
47 \label{WAVE FIG 1.1}
48 \end{figure}
49
50 Here we are modelling a point source at the point $x\hackscore C$ in the $x\hackscore{0}$-direction
51 which raise from zero to a maximum displacement $U\hackscore 0$ with in $\alpha$ seconds and then falls back to zero over the same period. In mathematical terms we use
52 \begin{eqnarray} \label{WAVE source}
53 u\hackscore 0(x\hackscore C,t)= U\hackscore 0 \frac{t^2}{\alpha^2} e^{1-\frac{t^2}{\alpha^2}} \
54 \end{eqnarray}
55 for all $t\ge0$. In the simulations we will choose $\alpha=0.3$, see Figure~\ref{WAVE FIG 1.2}
56 and will apply the source as a constraint in a in a sphere of small radius around the point
57 $x\hackscore C$.
58
59 We use an explicit time integration scheme to calculate the displacement field $u$ at
60 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$
61 which is defined by
62 \begin{eqnarray} \label{WAVE dyn 2}
63 u^{(n)}=2u^{(n-1)}-u^{(n-2)} + h^2 a^{(n)} \\
64 \end{eqnarray}
65 for all $n=2,3,\ldots$. It is designed to solve a system of equations of the form
66 \begin{eqnarray} \label{WAVE dyn 2b}
67 u\hackscore{,tt}=G(u)
68 \end{eqnarray}
69 where one sets $a^{(n)}=G(u^{(n-1)})$.
70
71 In our case $a^{(n)}$ is given by
72 \begin{eqnarray}\label{WAVE dyn 3}
73 \rho a^{(n)}\hackscore{i}=\sigma^{(n-1)}\hackscore{ij,j}
74 \end{eqnarray}
75 and boundary conditions
76 \begin{eqnarray} \label{WAVE natural at n}
77 \sigma^{(n-1)}\hackscore{ij}n\hackscore{j}=0
78 \end{eqnarray}
79 derived from \eqn{WAVE natural} where
80 \begin{eqnarray} \label{WAVE dyn 3a}
81 \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}).
82 \end{eqnarray}
83 We also need to apply the constraint
84 \begin{eqnarray} \label{WAVE dyn 3aa}
85 a^{(n)}\hackscore{0}(x\hackscore C,t)= U\hackscore{0}
86 \frac{2}{\alpha^2} \left( 1- 5 \frac{t^2}{\alpha^2} +2 \frac{t^4}{\alpha^4}
87 \right) e^{1-\frac{t^2}{\alpha^2}}
88 \end{eqnarray}
89 which is derived from equation~\ref{WAVE source} by calculating the second order time derivative,
90 see~\ref{WAVE FIG 1.2}. Now we have converted our problem for displacement, $u^{(n)}$, into a problem for
91 acceleration, $a^{(n)}$, which now depends
92 on the solution at the previous two time steps, $u^{(n-1)}$ and $u^{(n-2)}$.
93
94 In each time step we have to solve this problem to get the acceleration $a^{(n)}$, and we will
95 use the \LinearPDE class of the \linearPDEs to do so. The general form of the PDE defined through
96 the \LinearPDE class is discussed in \Sec{SEC LinearPDE}. The form which is relevant here is
97 \begin{equation}\label{WAVE dyn 100}
98 D\hackscore{ij} a^{(n)}\hackscore{j} = - X\hackscore{ij,j}\; .
99 \end{equation}
100 The natural boundary condition
101 \begin{equation}\label{WAVE dyn 101}
102 n\hackscore{j}X\hackscore{ij} =0
103 \end{equation}
104 is used.
105 With $u=a^{(n)}$ we can identify the values to be assigned to $D$ and $X$:
106 \begin{equation}\label{WAVE dyn 6}
107 \begin{array}{ll}
108 D\hackscore{ij}=\rho \delta\hackscore{ij}&
109 X\hackscore{ij}=-\sigma^{(n-1)}\hackscore{ij} \\
110 \end{array}
111 \end{equation}
112 Moreover we need to define the location $r$ where the constraint~\ref{WAVE dyn 3aa} is applied. We will apply
113 the constraint on a small sphere of radius $R$ around $x\hackscore C$ (we will use 3p.c. of the width of the domain):
114 \begin{equation}\label{WAVE dyn 6 1}
115 r\hackscore{i}(x) =
116 \left\{
117 \begin{array}{rc}
118 1 & \|x-x_c\|\le R \\
119 0 & \mbox{otherwise}
120 \end{array}
121 \right.
122 \end{equation}
123 The following script defines a the function \function{wavePropagation} which
124 implements the Verlet scheme to solve our wave propagation problem.
125 The argument \var{domain} which is a \Domain class object
126 defines the domain of the problem. \var{h} and \var{tend} are the time step size
127 and the end time of the simulation. \var{lam}, \var{mu} and
128 \var{rho} are material properties.
129 \begin{python}
130 def wavePropagation(domain,h,tend,lam,mu,rho,U0):
131 x=domain.getX()
132 # ... open new PDE ...
133 mypde=LinearPDE(domain)
134 mypde.getSolverOptions().setSolverMethod(mypde.getSolverOptions().LUMPING)
135 kronecker=identity(mypde.getDim())
136
137 # spherical source at middle of bottom face
138
139 xc=[width/2.,width/2.,0.]
140 # define small radius around point xc
141 # Lsup(x) returns the maximum value of the argument x
142 src_radius = 0.03*width
143 print "src_radius = ",src_radius
144
145 dunit=numpy.array([1.,0.,0.]) # defines direction of point source
146
147 mypde.setValue(D=kronecker*rho, q=whereNegative(length(x-xc)-src_radius)*dunit)
148 # ... set initial values ....
149 n=0
150 # initial value of displacement at point source is constant (U0=0.01)
151 # for first two time steps
152 u=Vector(0.,Solution(domain))
153 u_last=Vector(0.,Solution(domain))
154 t=0
155
156 # define the location of the point source
157 L=Locator(domain,xc)
158 # find potential at point source
159 u_pc=L.getValue(u)
160 print "u at point charge=",u_pc
161 # open file to save displacement at point source
162 u_pc_data=FileWriter('./data/U_pc.out')
163 u_pc_data.write("%f %f %f %f\n"%(t,u_pc[0],u_pc[1],u_pc[2]))
164
165 while t<tend:
166 t+=h
167 # ... get current stress ....
168 g=grad(u)
169 stress=lam*trace(g)*kronecker+mu*(g+transpose(g))
170 # ... get new acceleration ....
171 amplitude=U0*2*exp(1)/alpha**2*(1-5*(t/alpha)**2+2*(t/alpha)**4)*exp(-(t/alpha)**2)
172 mypde.setValue(X=-stress, r=dunit*amplitude)
173 a=mypde.getSolution()
174 # ... get new displacement ...
175 u_new=2*u-u_last+h**2*a
176 # ... shift displacements ....
177 u_last=u
178 u=u_new
179 n+=1
180 print n,"-th time step t ",t
181 u_pc=L.getValue(u)
182 print "u at point charge=",u_pc
183 # save displacements at point source to file for t > 0
184 u_pc_data.write("%f %f %f %f\n"%(t,u_pc[0],u_pc[1],u_pc[2]))
185
186 # ... save current acceleration in units of gravity and displacements
187 if n==1 or n%10==0: saveVTK("./data/usoln.%i.vtu"%(n/10),acceleration=length(a)/9.81,
188 displacement = length(u), tensor = stress, Ux = u[0] )
189
190 u_pc_data.close()
191 \end{python}
192 Notice that
193 all coefficients of the PDE which are independent of time $t$ are set outside the \code{while}
194 loop. This is very efficient since it allows the \LinearPDE class to reuse information as much as possible
195 when iterating over time.
196
197 The statement
198 \begin{python}
199 mypde.getSolverOptions().setSolverMethod(mypde.getSolverOptions().LUMPING)
200 \end{python}
201 switches on the use of an aggressive approximation of the PDE operator as a diagonal matrix
202 formed from the coefficient \var{D}.
203 The approximation allows, at the cost of
204 additional error, very fast
205 solution of the PDE. When using lumping the presence of \var{A}, \var{B} or \var{C} will produce wrong results.
206
207 There are a few new \escript functions in this example:
208 \code{grad(u)} returns the gradient $u\hackscore{i,j}$ of $u$ (in fact \var{grad(g)[i,j]}=$=u\hackscore{i,j}$).
209 There are restrictions on the argument of the \function{grad} function, for instance
210 the statement \code{grad(grad(u))} will raise an exception.
211 \code{trace(g)} returns the sum of the main diagonal elements \var{g[k,k]} of \var{g}
212 and \code{transpose(g)} returns the matrix transpose of \var{g} (ie. $\var{transpose(g)[i,j]}=\var{g[j,i]}$).
213
214 We initialize the values of \code{u} and \code{u_last} to be zero. It is important
215 to initialize both with the \SolutionFS \FunctionSpace as they have to be seen as solutions of PDEs from previous time steps. In fact, the \function{grad} does not accept arguments with a certain \FunctionSpace, for more details see \Sec{SEC ESCRIPT DATA}.
216
217 The \class{Locator} is designed to extract values at a given location (in this case $x\hackscore C$) from functions such as the displacement vector \code{u}. Typically the \class{Locator} is used in the following form:
218 \begin{python}
219 L=Locator(domain,xc)
220 u=...
221 u_pc=L.getValue(u)
222 \end{python}
223 The return value \code{u_pc} is the value of \code{u} at the location \code{xc}\footnote{In fact the finite element node which is closest to the given position. The usage of \class{Locator} is MPI save.}.
224
225
226 The output \code{U_pc.out} and \code{vtu} files are saved in a directory called \code{data}.
227 The \code{U_pc.out} stores 4 columns of data: $t,u\hackscore x,u\hackscore y,u\hackscore z$
228 respectively, where $u\hackscore x,u\hackscore y,u\hackscore z$ are the $x\hackscore 0,x\hackscore 1,x\hackscore 2$ components of
229 the displacement vector $u$ at the point source. These can be
230 plotted easily using any plotting package. In gnuplot the command:
231 \begin{verbatim}
232 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,
233 'U_pc.out' u 1:4 title 'U_z' w l lw 2
234 \end{verbatim}
235 will reproduce Figure~\ref{WAVE FIG 1}. It is pointed out that we are not using the
236 standart \PYTHON \function{open} to create the file \code{U_pc.out} as it is not
237 when running \escript under MPI, see chapter~\ref{EXECUTION} for more details.
238
239 One of the big advantages of the Verlet scheme is the fact that the problem to be solved
240 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).
241 The problem becomes so simple because we use the stress from the last time step rather then the stress which is
242 actually present at the current time step. Schemes using this approach are called an explicit time integration
243 schemes \index{explicit scheme} \index{time integration!explicit}. The
244 backward Euler scheme we have used in \Chap{DIFFUSION CHAP} is
245 an example of an implicit scheme
246 \index{implicit scheme} \index{time integration!implicit}. In this case one uses the actual status of
247 each variable at a particular time rather then values from previous time steps. This will lead to a problem which is
248 more expensive to solve, in particular for non-linear problems.
249 Although
250 explicit time integration schemes are cheap to finalize a single time step, they need significantly smaller time
251 steps then implicit schemes and can suffer from stability problems. Therefore they need a
252 very careful selection of the time step size $h$.
253
254 An easy, heuristic way of choosing an appropriate
255 time step size is the Courant condition \index{Courant condition} \index{explicit scheme!Courant condition}
256 which says that within a time step a information should not travel further than a cell used in the
257 discretization scheme. In the case of the wave equation the velocity of a (p-) wave is given as
258 $\sqrt{\frac{\lambda+2\mu}{\rho}}$ so one should choose $h$ from
259 \begin{eqnarray}\label{WAVE dyn 66}
260 h= \frac{1}{5} \sqrt{\frac{\rho}{\lambda+2\mu}} \Delta x
261 \end{eqnarray}
262 where $\Delta x$ is the cell diameter. The factor $\frac{1}{5}$ is a safety factor considering the heuristics of
263 the formula.
264
265 The following script uses the \code{wavePropagation} function to solve the
266 wave equation for a point source located at the bottom face of a block. The width of the block in
267 each direction on the bottom face is $10\mbox{km}$ ($x\hackscore 0$ and $x\hackscore 1$ directions ie. \code{l0} and \code{l1}).
268 The \code{ne} gives the number of elements in $x\hackscore{0}$ and $x\hackscore 1$ directions.
269 The depth of the block is aligned with the $x\hackscore{2}$-direction.
270 The depth (\code{l2}) of the block in
271 the $x\hackscore{2}$-direction is chosen so that there are 10 elements and the magnitude of the
272 of the depth is chosen such that the elements
273 become cubic. We chose 10 for the number of elements in $x\hackscore{2}$-direction so that the
274 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
275 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]$.
276 \begin{python}
277 from esys.finley import Brick
278 ne=32 # number of cells in x_0 and x_1 directions
279 width=10000. # length in x_0 and x_1 directions
280 lam=3.462e9
281 mu=3.462e9
282 rho=1154.
283 tend=60
284 U0=0.01 # amplitude of point source
285
286 mydomain=Brick(ne,ne,10,l0=width,l1=width,l2=10.*width/32.)
287 h=(1./5.)*inf(sqrt(rho/(lam+2*mu))*inf(domain.getSize())
288 print "time step size = ",h
289 wavePropagation(mydomain,h,tend,lam,mu,rho,U0)
290 \end{python}
291 The \function{domain.getSize()} return the local element size $\Delta x$. The
292 \function{inf} makes sure that the Courant condition~\ref{WAVE dyn 66} olds everywhere in the domain.
293
294 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
295 working directory.
296 To visualize the results from the data directory:
297 \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.
298 Again use \code{Configure Data} to move backwards and forward in time, and
299 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.
300
301 \begin{figure}[t!]
302 \centerline{\includegraphics[width=4.in]{figures/WavePC}}
303 \caption{Amplitude at Point source}
304 \label{WAVE FIG 1}
305 \end{figure}
306
307 \begin{figure}[t]
308 \begin{center}
309 \includegraphics[width=2in]{figures/Wavet0}
310 \includegraphics[width=2in]{figures/Wavet1}
311 \includegraphics[width=2in]{figures/Wavet3}
312 \includegraphics[width=2in]{figures/Wavet10}
313 \includegraphics[width=2in]{figures/Wavet30}
314 \includegraphics[width=2in]{figures/Wavet288}
315 \end{center}
316 \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
317 from a point source initially at $(5\mbox{km},5\mbox{km},0)$
318 with time step size $h=0.02083$. Color represents the displacement.
319 Here the view is oriented onto the bottom face.
320 \label{WAVE FIG 2}}
321 \end{figure}

Properties

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

  ViewVC Help
Powered by ViewVC 1.1.26