ViewVC logotype

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

Parent Directory Parent Directory | Revision Log Revision Log

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