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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 6686 - (show annotations)
Thu Jun 14 04:01:43 2018 UTC (2 years, 4 months ago) by aellery
File MIME type: application/x-tex
File size: 17768 byte(s)
Fixed some minor typos
1
2 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
3 % Copyright (c) 2003-2018 by The University of Queensland
4 % http://www.uq.edu.au
5 %
6 % Primary Business: Queensland, Australia
7 % Licensed under the Apache License, version 2.0
8 % http://www.apache.org/licenses/LICENSE-2.0
9 %
10 % Development until 2012 by Earth Systems Science Computational Center (ESSCC)
11 % Development 2012-2013 by School of Earth Sciences
12 % Development from 2014 by Centre for Geoscience Computing (GeoComp)
13 %
14 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
15
16 \section{Wave Propagation}
17 \label{WAVE CHAP}
18
19 In this next example we want to calculate the displacement field $u_{i}$ for any time $t>0$ by solving the wave equation:\index{wave equation}
20 \begin{eqnarray}\label{WAVE general problem}
21 \rho u_{i,tt} - \sigma_{ij,j}=0
22 \end{eqnarray}
23 in a three dimensional block of length $L$ in $x_{0}$
24 and $x_{1}$ direction and height $H$ in $x_{2}$ direction.
25 $\rho$ is the known density which may be a function of its location.
26 $\sigma_{ij}$ is the stress field\index{stress} which in case of an
27 isotropic, linear elastic material is given by
28 \begin{eqnarray}\label{WAVE stress}
29 \sigma_{ij} = \lambda u_{k,k} \delta_{ij} + \mu (u_{i,j} + u_{j,i})
30 \end{eqnarray}
31 where $\lambda$ and $\mu$ are the Lam\'e coefficients\index{Lam\'e coefficients}
32 and $\delta_{ij}$ denotes the Kronecker symbol\index{Kronecker symbol}.
33 On the boundary the normal stress is given by
34 \begin{eqnarray}\label{WAVE natural}
35 \sigma_{ij}n_{j}=0
36 \end{eqnarray}
37 for all time $t>0$.
38
39 \begin{figure}[t!]
40 \centerline{\includegraphics[width=14cm]{waveu}}
41 \caption{Input Displacement at Source Point ($\alpha=0.7$, $t_{0}=3$, $U_{0}=1$).}
42 \label{WAVE FIG 1.2}
43 \end{figure}
44
45 \begin{figure}
46 \centerline{\includegraphics[width=14cm]{wavea}}
47 \caption{Input Acceleration at Source Point ($\alpha=0.7$, $t_{0}=3$, $U_{0}=1$).}
48 \label{WAVE FIG 1.1}
49 \end{figure}
50
51 Here we are modelling a point source at the point $x_C$ in the
52 $x_{0}$-direction which is a negative pulse of amplitude
53 $U_{0}$ followed by the same positive pulse.
54 In mathematical terms we use
55 \begin{eqnarray} \label{WAVE source}
56 u_{0}(x_C,t)= U_{0} \; \sqrt{2} \frac{(t-t_{0})}{\alpha} e^{\frac{1}{2}-\frac{(t-t_{0})^2}{\alpha^2}} \
57 \end{eqnarray}
58 for all $t\ge0$ where $\alpha$ is the width of the pulse and $t_{0}$
59 is the time when the pulse changes from negative to positive.
60 In the simulations we will choose $\alpha=0.3$ and $t_{0}=2$ (see
61 \fig{WAVE FIG 1.2}) and apply the source as a constraint in a sphere of small
62 radius around the point $x_C$.
63
64 We use an explicit time integration scheme to calculate the displacement field
65 $u$ at certain time marks $t^{(n)}$, where $t^{(n)}=t^{(n-1)}+h$ with time
66 step size $h>0$.
67 In the following the upper index ${(n)}$ refers to values at time $t^{(n)}$.
68 We use the Verlet scheme\index{Verlet scheme} with constant time step size $h$
69 which is defined by
70 \begin{eqnarray} \label{WAVE dyn 2}
71 u^{(n)}=2u^{(n-1)}-u^{(n-2)} + h^2 a^{(n)} \\
72 \end{eqnarray}
73 for all $n=2,3,\ldots$. It is designed to solve a system of equations of the form
74 \begin{eqnarray} \label{WAVE dyn 2b}
75 u_{,tt}=G(u)
76 \end{eqnarray}
77 where one sets $a^{(n)}=G(u^{(n-1)})$.
78
79 In our case $a^{(n)}$ is given by
80 \begin{eqnarray}\label{WAVE dyn 3}
81 \rho a^{(n)}_{i}=\sigma^{(n-1)}_{ij,j}
82 \end{eqnarray}
83 and boundary conditions
84 \begin{eqnarray} \label{WAVE natural at n}
85 \sigma^{(n-1)}_{ij}n_{j}=0
86 \end{eqnarray}
87 derived from \eqn{WAVE natural} where
88 \begin{eqnarray} \label{WAVE dyn 3a}
89 \sigma_{ij}^{(n-1)} & = & \lambda u^{(n-1)}_{k,k} \delta_{ij} + \mu ( u^{(n-1)}_{i,j} + u^{(n-1)}_{j,i}).
90 \end{eqnarray}
91 We also need to apply the constraint
92 \begin{eqnarray} \label{WAVE dyn 3aa}
93 a^{(n)}_{0}(x_C,t)= U_{0}
94 \frac{\sqrt(2.)}{\alpha^2} (4\frac{(t-t_{0})^3}{\alpha^3}-6\frac{t-t_{0}}{\alpha})e^{\frac{1}{2}-\frac{(t-t_{0})^2}{\alpha^2}}
95 \end{eqnarray}
96 which is derived from equation~\ref{WAVE source} by calculating the second
97 order time derivative (see \fig{WAVE FIG 1.1}).
98 Now we have converted our problem for displacement, $u^{(n)}$, into a problem
99 for acceleration, $a^{(n)}$, which depends on the solution at the previous two
100 time steps $u^{(n-1)}$ and $u^{(n-2)}$.
101
102 In each time step we have to solve this problem to get the acceleration $a^{(n)}$,
103 and we will use the \LinearPDE class of the \linearPDEs package to do so.
104 The general form of the PDE defined through the \LinearPDE class is discussed
105 in \Sec{SEC LinearPDE}. The form which is relevant here is
106 \begin{equation}\label{WAVE dyn 100}
107 D_{ij} a^{(n)}_{j} = - X_{ij,j}\; .
108 \end{equation}
109 The natural boundary condition
110 \begin{equation}\label{WAVE dyn 101}
111 n_{j}X_{ij} =0
112 \end{equation}
113 is used.
114 With $u=a^{(n)}$ we can identify the values to be assigned to $D$ and $X$:
115 \begin{equation}\label{WAVE dyn 6}
116 \begin{array}{ll}
117 D_{ij}=\rho \delta_{ij}&
118 X_{ij}=-\sigma^{(n-1)}_{ij}\\
119 \end{array}
120 \end{equation}
121 Moreover we need to define the location $r$ where the constraint~\ref{WAVE dyn 3aa} is applied.
122 We will apply the constraint on a small sphere of radius $R$ around
123 $x_C$ (we will use 3\% of the width of the domain):
124 \begin{equation}\label{WAVE dyn 6 1}
125 q_{i}(x) =
126 \left\{
127 \begin{array}{l l}
128 1 & \quad \text{where $\|x-x_c\|\le R$}\\
129 0 & \quad \text{otherwise.}
130 \end{array}
131 \right.
132 \end{equation}
133 The following script defines the function \function{wavePropagation} which
134 implements the Verlet scheme to solve our wave propagation problem.
135 The argument \var{domain} which is a \Domain class object defines the domain
136 of the problem.
137 \var{h} and \var{tend} are the time step size and the end time of the simulation.
138 \var{lam}, \var{mu} and \var{rho} are material properties.
139 \begin{python}
140 def wavePropagation(domain,h,tend,lam,mu,rho, x_c, src_radius, U0):
141 # lists to collect displacement at point source which is returned
142 # to the caller
143 ts, u_pc0, u_pc1, u_pc2 = [], [], [], []
144
145 x=domain.getX()
146 # ... open new PDE ...
147 mypde=LinearPDE(domain)
148 mypde.getSolverOptions().setSolverMethod(SolverOptions.HRZ_LUMPING)
149 kronecker=identity(mypde.getDim())
150 dunit=numpy.array([1., 0., 0.]) # defines direction of point source
151 mypde.setValue(D=kronecker*rho, q=whereNegative(length(x-xc)-src_radius)*dunit)
152 # ... set initial values ....
153 n=0
154 # for first two time steps
155 u=Vector(0., Solution(domain))
156 u_last=Vector(0., Solution(domain))
157 t=0
158 # define the location of the point source
159 L=Locator(domain, xc)
160 # find potential at point source
161 u_pc=L.getValue(u)
162 print("u at point charge=",u_pc)
163 ts.append(t)
164 u_pc0.append(u_pc[0])
165 u_pc1.append(u_pc[1])
166 u_pc2.append(u_pc[2])
167
168 while t<tend:
169 t+=h
170 # ... get current stress ....
171 g=grad(u)
172 stress=lam*trace(g)*kronecker+mu*(g+transpose(g))
173 # ... get new acceleration ....
174 amplitude=U0*(4*(t-t0)**3/alpha**3-6*(t-t0)/alpha)*sqrt(2.)/alpha**2 \
175 *exp(1./2.-(t-t0)**2/alpha**2)
176 mypde.setValue(X=-stress, r=dunit*amplitude)
177 a=mypde.getSolution()
178 # ... get new displacement ...
179 u_new=2*u-u_last+h**2*a
180 # ... shift displacements ....
181 u_last=u
182 u=u_new
183 n+=1
184 print(n,"-th time step, t=",t)
185 u_pc=L.getValue(u)
186 print("u at point charge=",u_pc)
187 # save displacements at point source to file for t > 0
188 ts.append(t)
189 u_pc0.append(u_pc[0])
190 u_pc1.append(u_pc[1])
191 u_pc2.append(u_pc[2])
192
193 # ... save current acceleration in units of gravity and displacements
194 if n==1 or n%10==0:
195 saveVTK("./data/usoln.%i.vtu"%(n/10), \
196 acceleration = length(a)/9.81, \
197 displacement = length(u), \
198 tensor = stress, Ux = u[0])
199
200 return ts, u_pc0, u_pc1, u_pc2
201 \end{python}
202 Notice that all coefficients of the PDE which are independent of time $t$ are
203 set outside the \code{while} loop.
204 This is for efficiency reasons since it allows the \LinearPDE class to reuse
205 information as much as possible when iterating over time.
206
207 The statement
208 \begin{python}
209 mypde.getSolverOptions().setSolverMethod(SolverOptions.HRZ_LUMPING)
210 \end{python}
211 enables the use of an aggressive approximation of the PDE operator as a
212 diagonal matrix formed from the coefficient \var{D}.
213 The approximation allows, at the cost of additional error, very fast solution
214 of the PDE, see also Section~\ref{LUMPING}.
215
216 There are a few new \escript functions in this example:
217 \code{grad(u)} returns the gradient $u_{i,j}$ of $u$ (in fact \var{grad(g)[i,j]}$==u_{i,j}$).
218 There are restrictions on the argument of the \function{grad} function, for
219 instance the statement \code{grad(grad(u))} will raise an exception.
220 \code{trace(g)} returns the sum of the main diagonal elements \var{g[k,k]} of
221 \var{g} and \code{transpose(g)} returns the matrix transpose of \var{g} (i.e.
222 $\var{transpose(g)[i,j]}==\var{g[j,i]}$).
223
224 We initialize the values of \code{u} and \code{u_last} to be zero.
225 It is important to initialize both with the \SolutionFS as they have to be
226 seen as solutions of PDEs from previous time steps.
227 In fact, the \function{grad} does not accept arguments with a certain
228 \FunctionSpace, for more details see \Sec{SEC ESCRIPT DATA}.
229
230 The \class{Locator} class is designed to extract values at a given location
231 (in this case $x_C$) from functions such as the displacement vector \code{u}.
232 Typically \class{Locator} is used in the following way:
233 \begin{python}
234 L=Locator(domain, xc)
235 u=...
236 u_pc=L.getValue(u)
237 \end{python}
238 The return value \code{u_pc} is the value of \code{u} at the location
239 \code{xc}\footnote{In fact, it is the finite element node which is closest to
240 the given position. The usage of \class{Locator} is MPI safe.}.
241 The values are collected in the lists \var{u_pc0}, \var{u_pc1} and \var{u_pc2}
242 together with the corresponding time marker in \var{ts}.
243 These values are handed back to the caller. Later we will show ways to access these data.
244
245 One of the big advantages of the Verlet scheme is the fact that the problem to
246 be solved in each time step is very simple and does not involve any spatial
247 derivatives (which is what allows us to use lumping in this simulation).
248 The problem becomes so simple because we use the stress from the last time
249 step rather than the stress which is actually present at the current time step.
250 Schemes using this approach are called explicit time integration schemes\index{explicit scheme}\index{time integration!explicit}.
251 The backward Euler scheme we have used in \Chap{DIFFUSION CHAP} is an example
252 of an implicit scheme\index{implicit scheme}\index{time integration!implicit}.
253 In this case one uses the actual status of each variable at a particular time
254 rather than values from previous time steps.
255 This will lead to a problem which is more expensive to solve, in particular
256 for non-linear cases.
257 Although explicit time integration schemes are cheap to finalize a single time
258 step, they need significantly smaller time steps than implicit schemes and can
259 suffer from stability problems.
260 Therefore they require a very careful selection of the time step size $h$.
261
262 An easy, heuristic way of choosing an appropriate time step size is the
263 Courant-Friedrichs-Lewy condition (CFL condition)\index{Courant condition}\index{explicit scheme!Courant condition}
264 which says that within a time step information should not travel further than
265 a cell used in the discretization scheme.
266 In the case of the wave equation the velocity of a (p-) wave is given as
267 $\sqrt{\frac{\lambda+2\mu}{\rho}}$ so one should choose $h$ from
268 \begin{eqnarray}\label{WAVE dyn 66}
269 h= \frac{1}{5} \sqrt{\frac{\rho}{\lambda+2\mu}} \Delta x
270 \end{eqnarray}
271 where $\Delta x$ is the cell diameter.
272 The factor $\frac{1}{5}$ is a safety factor considering the heuristics of the formula.
273
274 \begin{figure}[t]
275 \begin{center}
276 \includegraphics[width=2in]{Wave11}
277 \includegraphics[width=2in]{Wave22}
278 \includegraphics[width=2in]{Wave28}
279 \includegraphics[width=2in]{Wave32}
280 \includegraphics[width=2in]{Wave36}
281 \end{center}
282 \caption{Selected time steps ($n = 11, 22, 32, 36$) of a wave propagation over
283 a 10km $\times$ 10km $\times$ 3.125km block from a point source initially at
284 (5km, 5km, 0) with time step size $h=0.02083$. Color represents the
285 displacement. Here the view is oriented onto the bottom face.
286 \label{WAVE FIG 2}}
287 \end{figure}
288
289 The following script uses the \code{wavePropagation} function to solve the
290 wave equation for a point source located at the bottom face of a block.
291 The width of the block in each direction on the bottom face is 10km
292 ($x_0$ and $x_1$ directions, i.e. \code{l0} and \code{l1}).
293 The variable \code{ne} gives the number of elements in $x_{0}$ and $x_1$ directions.
294 The depth of the block is aligned with the $x_{2}$-direction.
295 The depth (\code{l2}) of the block in the $x_{2}$-direction is chosen
296 so that there are 10 elements, and the magnitude of the depth is chosen such
297 that the elements become cubic.
298 We chose 10 for the number of elements in the $x_{2}$-direction so
299 that the computation is faster. \code{Brick($n_0, n_1, n_2,l_0,l_1,l_2$)}
300 is an \finley function which creates a rectangular mesh with $n_0 \times n_1 \times n_2$
301 elements over the brick $[0,l_0] \times [0,l_1] \times [0,l_2]$.
302 \begin{python}
303 from esys.finley import Brick
304 ne = 32 # number of cells in x_0 and x_1 directions
305 width = 10000. # length in x_0 and x_1 directions
306 lam = 3.462e9
307 mu = 3.462e9
308 rho = 1154.
309 tend = 60
310 U0 = 1. # amplitude of point source
311 # spherical source at middle of bottom face
312 xc=[width/2.,width/2.,0.]
313 # define small radius around point xc
314 src_radius = 0.03*width
315 print("src_radius =",src_radius)
316 mydomain=Brick(ne, ne, 10, l0=width, l1=width, l2=10.*width/32.)
317 h=(1./5.)*inf(sqrt(rho/(lam+2*mu)))*inf(domain.getSize())
318 print("time step size =",h)
319 ts, u_pc0, u_pc1, u_pc2 = \
320 wavePropagation(mydomain, h, tend, lam, mu, rho, xc, src_radius, U0)
321 \end{python}
322 The \function{domain.getSize()} function returns the local element size $\Delta x$.
323 Using \function{inf} ensures that the CFL condition~\ref{WAVE dyn 66} \index{Courant number}\index{CFL condition} holds
324 everywhere in the domain.
325
326 The script is available as \file{wave.py} in the \ExampleDirectory\index{scripts!\file{wave.py}}.
327 To visualize the results from the data directory:
328 \begin{verbatim}
329 mayavi2 -d usoln.1.vtu -m Surface
330 \end{verbatim}
331 You can rotate this figure by clicking on it with the mouse and moving it around.
332 Again use \emph{Configure Data} to move backward and forward in time, and
333 also to choose the results (acceleration, displacement or $u_x$) by
334 using \emph{Select Scalar}.
335 \fig{WAVE FIG 2} shows the results for the displacement at various time steps.
336
337 \begin{figure}[t!]
338 \centerline{\includegraphics[width=14cm]{WavePC}}
339 \caption{Amplitude at Point source from the Simulation}
340 \label{WAVE FIG 1}
341 \end{figure}
342
343 It remains to show some possibilities to inspect the collected data
344 \var{u_pc0}, \var{u_pc1} and \var{u_pc2}.
345 One way is to write the data to a file and then use an external package such
346 as \gnuplot, LibreOffice Calc or Excel to read the data for further analysis.
347 The following code shows one possible way to write the data to the file \file{./data/U_pc.csv}:
348 \begin{python}
349 u_pc_data=FileWriter('./data/U_pc.csv')
350 for i in range(len(ts)):
351 u_pc_data.write("%f %f %f %f\n"%(ts[i],u_pc0[i],u_pc1[i],u_pc2[i]))
352 u_pc_data.close()
353 \end{python}
354 The file \code{U_pc.csv} stores 4 columns of data: $t,u_x,u_y,u_z$
355 respectively, where $u_x,u_y,u_z$ are the
356 $x_0,x_1,x_2$ components of the displacement
357 vector $u$ at the point source.
358 These can be plotted easily using any plotting package.
359 In \gnuplot the command:
360 \begin{verbatim}
361 plot 'U_pc.csv' u 1:2 title 'U_x' w l lw 2, 'U_pc.csv' \
362 u 1:3 title 'U_y' w l lw 2, \
363 'U_pc.csv' u 1:4 title 'U_z' w l lw 2
364 \end{verbatim}
365 will reproduce \fig{WAVE FIG 1} (As expected this is identical to the input
366 signal shown in \fig{WAVE FIG 1.2}).
367 It is pointed out that we are not using the standard \PYTHON \function{open}
368 to write to the file \code{U_pc.csv} as it is not safe when running \escript
369 under MPI, see Chapter~\ref{EXECUTION} for more details.
370
371 Alternatively, one can implement plotting the results at run time rather than
372 in a post-processing step.
373 This avoids the generation of an intermediate data file.
374 In {\it escript} the preferred way of creating 2D plots of time dependent data
375 is \MATPLOTLIB. The following script creates the plot and writes it into the
376 file \file{u_pc.png} in the PNG image format:
377 \begin{python}
378 import matplotlib.pyplot as plt
379 if getMPIRankWorld() == 0:
380 plt.title("Displacement at Point Source")
381 plt.plot(ts, u_pc0, '-', label="x_0", linewidth=1)
382 plt.plot(ts, u_pc1, '-', label="x_1", linewidth=1)
383 plt.plot(ts, u_pc2, '-', label="x_2", linewidth=1)
384 plt.xlabel('time')
385 plt.ylabel('displacement')
386 plt.legend()
387 plt.savefig('u_pc.png', format='png')
388 \end{python}
389 You can add \function{plt.show()} to create an interactive browser window.
390 Notice that by checking the condition \code{getMPIRankWorld()==0} the plot
391 is generated on one processor only (in this case the rank 0 processor) when
392 run under \MPI.
393
394 Both options for processing the point source data are included in the example
395 file \file{wave.py}. There are other options available to process these data
396 in particular through the \SCIPY package, e.g. Fourier transformations.
397 It is beyond the scope of this user's guide to document the usage of
398 \SCIPY for time series analysis but it is highly recommended to look in
399 relevant readily available documentation.
400

Properties

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

  ViewVC Help
Powered by ViewVC 1.1.26