ViewVC logotype

Contents of /branches/stage3.0/doc/user/wave.tex

Parent Directory Parent Directory | Revision Log Revision Log

Revision 2584 - (show annotations)
Wed Aug 5 02:44:51 2009 UTC (9 years, 6 months ago) by jfenwick
File MIME type: application/x-tex
File size: 18547 byte(s)
Merging Lutz' doco changes to release
2 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
3 %
4 % Copyright (c) 2003-2009 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 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
15 \section{3-D Wave Propagation}
16 \label{WAVE CHAP}
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$.
38 \begin{figure}[t!]
39 \centerline{\includegraphics[angle=-90,width=4.in]{figures/waveu}}
40 \caption{Input Displacement at Source Point ($\alpha=0.7$, $t\hackscore{0}=3$, $U\hackscore{0}=1$).}
41 \label{WAVE FIG 1.2}
42 \end{figure}
44 \begin{figure}[t!]
45 \centerline{\includegraphics[angle=-90,width=4.in]{figures/wavea}}
46 \caption{Input Acceleration at Source Point ($\alpha=0.7$, $t\hackscore{0}=3$, $U\hackscore{0}=1$).}
47 \label{WAVE FIG 1.1}
48 \end{figure}
50 Here we are modelling a point source at the point $x\hackscore C$ in the $x\hackscore{0}$-direction
51 which is a negative pulse of amplitude $U\hackscore{0}$ followed by the same
52 positive pulse. In mathematical terms we use
53 \begin{eqnarray} \label{WAVE source}
54 u\hackscore{0}(x\hackscore C,t)= U\hackscore{0} \; \sqrt{2} \frac{(t-t\hackscore{0})}{\alpha} e^{\frac{1}{2}-\frac{(t-t\hackscore{0})^2}{\alpha^2}} \
55 \end{eqnarray}
56 for all $t\ge0$ where $\alpha$ is the width of the pulse and $t\hackscore{0}$ is the time when
57 the pulse changes from negative to positive. In the simulations we will choose $\alpha=0.3$ and $t\hackscore{0}=2$, see Figure~\ref{WAVE FIG 1.2}
58 and will apply the source as a constraint in a in a sphere of small radius around the point
59 $x\hackscore C$.
61 We use an explicit time integration scheme to calculate the displacement field $u$ at
62 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$
63 which is defined by
64 \begin{eqnarray} \label{WAVE dyn 2}
65 u^{(n)}=2u^{(n-1)}-u^{(n-2)} + h^2 a^{(n)} \\
66 \end{eqnarray}
67 for all $n=2,3,\ldots$. It is designed to solve a system of equations of the form
68 \begin{eqnarray} \label{WAVE dyn 2b}
69 u\hackscore{,tt}=G(u)
70 \end{eqnarray}
71 where one sets $a^{(n)}=G(u^{(n-1)})$.
73 In our case $a^{(n)}$ is given by
74 \begin{eqnarray}\label{WAVE dyn 3}
75 \rho a^{(n)}\hackscore{i}=\sigma^{(n-1)}\hackscore{ij,j}
76 \end{eqnarray}
77 and boundary conditions
78 \begin{eqnarray} \label{WAVE natural at n}
79 \sigma^{(n-1)}\hackscore{ij}n\hackscore{j}=0
80 \end{eqnarray}
81 derived from \eqn{WAVE natural} where
82 \begin{eqnarray} \label{WAVE dyn 3a}
83 \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}).
84 \end{eqnarray}
85 We also need to apply the constraint
86 \begin{eqnarray} \label{WAVE dyn 3aa}
87 a^{(n)}\hackscore{0}(x\hackscore C,t)= U\hackscore{0}
88 \frac{\sqrt(2.)}{\alpha^2} (4\frac{(t-t\hackscore{0})^3}{\alpha^3}-6\frac{t-t\hackscore{0}}{\alpha})e^{\frac{1}{2}-\frac{(t-t\hackscore{0})^2}{\alpha^2}}
89 \end{eqnarray}
90 which is derived from equation~\ref{WAVE source} by calculating the second order time derivative,
91 see Figure~\ref{WAVE FIG 1.1}. Now we have converted our problem for displacement, $u^{(n)}$, into a problem for
92 acceleration, $a^{(n)}$, which now depends
93 on the solution at the previous two time steps, $u^{(n-1)}$ and $u^{(n-2)}$.
95 In each time step we have to solve this problem to get the acceleration $a^{(n)}$, and we will
96 use the \LinearPDE class of the \linearPDEs to do so. The general form of the PDE defined through
97 the \LinearPDE class is discussed in \Sec{SEC LinearPDE}. The form which is relevant here is
98 \begin{equation}\label{WAVE dyn 100}
99 D\hackscore{ij} a^{(n)}\hackscore{j} = - X\hackscore{ij,j}\; .
100 \end{equation}
101 The natural boundary condition
102 \begin{equation}\label{WAVE dyn 101}
103 n\hackscore{j}X\hackscore{ij} =0
104 \end{equation}
105 is used.
106 With $u=a^{(n)}$ we can identify the values to be assigned to $D$ and $X$:
107 \begin{equation}\label{WAVE dyn 6}
108 \begin{array}{ll}
109 D\hackscore{ij}=\rho \delta\hackscore{ij}&
110 X\hackscore{ij}=-\sigma^{(n-1)}\hackscore{ij} \\
111 \end{array}
112 \end{equation}
113 Moreover we need to define the location $r$ where the constraint~\ref{WAVE dyn 3aa} is applied. We will apply
114 the constraint on a small sphere of radius $R$ around $x\hackscore C$ (we will use 3p.c. of the width of the domain):
115 \begin{equation}\label{WAVE dyn 6 1}
116 q\hackscore{i}(x) =
117 \left\{
118 \begin{array}{rc}
119 1 & \|x-x_c\|\le R \\
120 0 & \mbox{otherwise}
121 \end{array}
122 \right.
123 \end{equation}
124 The following script defines a the function \function{wavePropagation} which
125 implements the Verlet scheme to solve our wave propagation problem.
126 The argument \var{domain} which is a \Domain class object
127 defines the domain of the problem. \var{h} and \var{tend} are the time step size
128 and the end time of the simulation. \var{lam}, \var{mu} and
129 \var{rho} are material properties.
130 \begin{python}
131 def wavePropagation(domain,h,tend,lam,mu,rho, x_c, src_radius, U0):
132 # lists to collect displacement at point source which is returned to the caller
133 ts, u_pc0,u_pc1,u_pc2=[], [], [], []
135 x=domain.getX()
136 # ... open new PDE ...
137 mypde=LinearPDE(domain)
138 mypde.getSolverOptions().setSolverMethod(mypde.getSolverOptions().LUMPING)
139 kronecker=identity(mypde.getDim())
140 dunit=numpy.array([1.,0.,0.]) # defines direction of point source
141 mypde.setValue(D=kronecker*rho, q=whereNegative(length(x-xc)-src_radius)*dunit)
142 # ... set initial values ....
143 n=0
144 # for first two time steps
145 u=Vector(0.,Solution(domain))
146 u_last=Vector(0.,Solution(domain))
147 t=0
148 # define the location of the point source
149 L=Locator(domain,xc)
150 # find potential at point source
151 u_pc=L.getValue(u)
152 print "u at point charge=",u_pc
153 # open file to save displacement at point source
154 u_pc_data=FileWriter('./data/U_pc.out')
155 ts.append(t); u_pc0.append(u_pc[0]), u_pc1.append(u_pc[1]), u_pc2.append(u_pc[2])
157 while t<tend:
158 t+=h
159 # ... get current stress ....
160 g=grad(u)
161 stress=lam*trace(g)*kronecker+mu*(g+transpose(g))
162 # ... get new acceleration ....
163 amplitude=U0*(4*(t-t0)**3/alpha**3-6*(t-t0)/alpha)*sqrt(2.)/alpha**2 \
164 *exp(1./2.-(t-t0)**2/alpha**2)
165 mypde.setValue(X=-stress, r=dunit*amplitude)
166 a=mypde.getSolution()
167 # ... get new displacement ...
168 u_new=2*u-u_last+h**2*a
169 # ... shift displacements ....
170 u_last=u
171 u=u_new
172 n+=1
173 print n,"-th time step t ",t
174 u_pc=L.getValue(u)
175 print "u at point charge=",u_pc
176 # save displacements at point source to file for t > 0
177 ts.append(t); u_pc0.append(u_pc[0]), u_pc1.append(u_pc[1]), \
178 u_pc2.append(u_pc[2])
180 # ... save current acceleration in units of gravity and displacements
181 if n==1 or n%10==0: saveVTK("./data/usoln.%i.vtu"%(n/10), \
182 acceleration=length(a)/9.81,
183 displacement = length(u), \
184 tensor = stress, Ux = u[0] )
186 return ts, u_pc0, u_pc1, u_pc2
187 \end{python}
188 Notice that
189 all coefficients of the PDE which are independent of time $t$ are set outside the \code{while}
190 loop. This is very efficient since it allows the \LinearPDE class to reuse information as much as possible
191 when iterating over time.
193 The statement
194 \begin{python}
195 mypde.getSolverOptions().setSolverMethod(mypde.getSolverOptions().LUMPING)
196 \end{python}
197 switches on the use of an aggressive approximation of the PDE operator as a diagonal matrix
198 formed from the coefficient \var{D}.
199 The approximation allows, at the cost of
200 additional error, very fast
201 solution of the PDE. When using lumping the presence of \var{A}, \var{B} or \var{C} will produce wrong results.
203 There are a few new \escript functions in this example:
204 \code{grad(u)} returns the gradient $u\hackscore{i,j}$ of $u$ (in fact \var{grad(g)[i,j]}=$=u\hackscore{i,j}$).
205 There are restrictions on the argument of the \function{grad} function, for instance
206 the statement \code{grad(grad(u))} will raise an exception.
207 \code{trace(g)} returns the sum of the main diagonal elements \var{g[k,k]} of \var{g}
208 and \code{transpose(g)} returns the matrix transpose of \var{g} (ie. $\var{transpose(g)[i,j]}=\var{g[j,i]}$).
210 We initialize the values of \code{u} and \code{u_last} to be zero. It is important
211 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}.
213 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:
214 \begin{python}
215 L=Locator(domain,xc)
216 u=...
217 u_pc=L.getValue(u)
218 \end{python}
219 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.}. The values
220 are collected in the lists \var{u_pc0}, \var{u_pc1} and \var{u_pc2} together with the
221 corresponding time marker in \var{ts}. The values are handed back to the caller. Later we will show to ways to
222 access these data.
224 One of the big advantages of the Verlet scheme is the fact that the problem to be solved
225 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).
226 The problem becomes so simple because we use the stress from the last time step rather then the stress which is
227 actually present at the current time step. Schemes using this approach are called an explicit time integration
228 schemes \index{explicit scheme} \index{time integration!explicit}. The
229 backward Euler scheme we have used in \Chap{DIFFUSION CHAP} is
230 an example of an implicit scheme
231 \index{implicit scheme} \index{time integration!implicit}. In this case one uses the actual status of
232 each variable at a particular time rather then values from previous time steps. This will lead to a problem which is
233 more expensive to solve, in particular for non-linear problems.
234 Although
235 explicit time integration schemes are cheap to finalize a single time step, they need significantly smaller time
236 steps then implicit schemes and can suffer from stability problems. Therefore they need a
237 very careful selection of the time step size $h$.
239 An easy, heuristic way of choosing an appropriate
240 time step size is the Courant condition \index{Courant condition} \index{explicit scheme!Courant condition}
241 which says that within a time step a information should not travel further than a cell used in the
242 discretization scheme. In the case of the wave equation the velocity of a (p-) wave is given as
243 $\sqrt{\frac{\lambda+2\mu}{\rho}}$ so one should choose $h$ from
244 \begin{eqnarray}\label{WAVE dyn 66}
245 h= \frac{1}{5} \sqrt{\frac{\rho}{\lambda+2\mu}} \Delta x
246 \end{eqnarray}
247 where $\Delta x$ is the cell diameter. The factor $\frac{1}{5}$ is a safety factor considering the heuristics of
248 the formula.
250 \begin{figure}[t]
251 \begin{center}
252 \includegraphics[width=2in]{figures/Wave11}
253 \includegraphics[width=2in]{figures/Wave22}
254 \includegraphics[width=2in]{figures/Wave28}
255 \includegraphics[width=2in]{figures/Wave32}
256 \includegraphics[width=2in]{figures/Wave36}
257 \end{center}
258 \caption{Selected time steps ($n = 11, 22, 32, 36$) of a wave propagation over a $10\mbox{km} \times 10\mbox{km} \times 3.125\mbox{km}$ block
259 from a point source initially at $(5\mbox{km},5\mbox{km},0)$
260 with time step size $h=0.02083$. Color represents the displacement.
261 Here the view is oriented onto the bottom face.
262 \label{WAVE FIG 2}}
263 \end{figure}
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=1. # amplitude of point source
285 # spherical source at middle of bottom face
286 xc=[width/2.,width/2.,0.]
287 # define small radius around point xc
288 src_radius = 0.03*width
289 print "src_radius = ",src_radius
290 mydomain=Brick(ne,ne,10,l0=width,l1=width,l2=10.*width/32.)
291 h=(1./5.)*inf(sqrt(rho/(lam+2*mu))*inf(domain.getSize())
292 print "time step size = ",h
293 ts, u_pc0, u_pc1, u_pc2 = \
294 wavePropagation(mydomain,h,tend,lam,mu,rho,xc, src_radius, U0)
295 \end{python}
296 The \function{domain.getSize()} return the local element size $\Delta x$. The
297 \function{inf} makes sure that the Courant condition~\ref{WAVE dyn 66} olds everywhere in the domain.
299 The script is available as \file{wave.py} in the \ExampleDirectory \index{scripts!\file{wave.py}}.
300 To visualize the results from the data directory:
301 \begin{verbatim}
302 mayavi2 -d usoln.1.vtu -m SurfaceMap &
303 \end{verbatim}
304 You can rotate this figure by clicking on it with the mouse and moving it around.
305 Again use \code{Configure Data} to move backwards and forward in time, and
306 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.
308 \begin{figure}[t!]
309 \centerline{\includegraphics[width=4.in, angle=-90]{figures/WavePC}}
310 \caption{Amplitude at Point source from the Simulation.}
311 \label{WAVE FIG 1}
312 \end{figure}
314 It remains to show some possibilities to inspect the collected data \var{u_pc0}, \var{u_pc1} and \var{u_pc2}.
315 One way is to write the data to a file and then use an external package such as \gnuplot, excel or OpenOffice to read the data for further analysis. The following code shows one possible form to write the data to the
316 file \file{./data/U_pc.out}:
317 \begin{python}
318 u_pc_data=FileWriter('./data/U_pc.out')
319 for i in xrange(len(ts)) :
320 u_pc_data.write("%f %f %f %f\n"%(ts[i],u_pc0[i],u_pc1[i],u_pc2[i]))
321 u_pc_data.close()
322 \end{python}
323 The \code{U_pc.out} stores 4 columns of data: $t,u\hackscore x,u\hackscore y,u\hackscore z$
324 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
325 plotted easily using any plotting package. In \gnuplot the command:
326 \begin{verbatim}
327 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,
328 'U_pc.out' u 1:4 title 'U_z' w l lw 2
329 \end{verbatim}
330 will reproduce Figure~\ref{WAVE FIG 1} (As expected this is identical to the input signal shown in Figure~\ref{WAVE FIG 1.2})2. It is pointed out that we are not using the
331 standard \PYTHON \function{open} to write to the file \code{U_pc.out} as it is not safe
332 when running \escript under MPI, see chapter~\ref{EXECUTION} for more details.
334 Alternatively, one can implement plotting the results at run time rather than in a post-processing step. This avoids
335 the generation of an intermediate data file. In {\it escript} the preferred way of creating 2D plots of
336 time dependent data is \MATPLOTLIB. The following script creates the plot and writes it into the
337 file \file{u_pc.png} in the PNG image format:
338 \begin{python}
339 import matplotlib.pyplot as plt
340 if getMPIRankWorld() == 0:
341 plt.title("Displacement at Point Source")
342 plt.plot(ts, u_pc0, '-', label="x_0", linewidth=1)
343 plt.plot(ts, u_pc1, '-', label="x_1", linewidth=1)
344 plt.plot(ts, u_pc2, '-', label="x_2", linewidth=1)
345 plt.xlabel('time')
346 plt.ylabel('displacement')
347 plt.legend()
348 plt.savefig('u_pc.png', format='png')
349 \end{python}
350 You can add the \function{plt.show()} to create a interactive browser window. Please not that
351 through the \code{getMPIRankWorld() == 0} statement the plot is generated on one processor only (in this case the rank 0 processor) when run under MPI.
353 Both options for processing the point source data are include in the example file \file{wave.py}. There other options available to process these data in particular through the \SCIPY package , eg Fourier transformations. It is beyond the scope of this users guide to document the usage of \SCIPY for time series analysis but is highly recommended that users use \SCIPY to any further data processing.


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

  ViewVC Help
Powered by ViewVC 1.1.26