ViewVC logotype

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

Parent Directory Parent Directory | Revision Log Revision Log

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


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

  ViewVC Help
Powered by ViewVC 1.1.26