# Diff of /trunk/doc/cookbook/example07.tex

revision 3003 by ahallam, Wed Apr 7 02:29:57 2010 UTC revision 3054 by ahallam, Wed Jun 30 02:22:25 2010 UTC
# Line 15  Line 15
15
16  The acoustic wave equation governs the propagation of pressure waves. Wave  The acoustic wave equation governs the propagation of pressure waves. Wave
17  types that obey this law tend to travel in liquids or gases where shear waves  types that obey this law tend to travel in liquids or gases where shear waves
18  or longitudinal style wave motion is not possible. The obvious example is sound  or longitudinal style wave motion is not possible. An obvious example is sound
19  waves.  waves.
20
21  The acoustic wave equation is;  The acoustic wave equation is defined as;
22
23   \nabla ^2 p - \frac{1}{c^2} \frac{\partial ^2 p}{\partial t^2} = 0   \nabla ^2 p - \frac{1}{c^2} \frac{\partial ^2 p}{\partial t^2} = 0
24  \label{eqn:acswave}  \label{eqn:acswave}
25
26  where $p$ is the pressure, $t$ is the time and $c$ is the wave velocity.  where $p$ is the pressure, $t$ is the time and $c$ is the wave velocity.

27
28    \section{The Laplacian in \esc}
29    The Laplacian opperator which can be written as $\Delta$ or $\nabla^2$  is
30    calculated via the divergence of the gradient of the object, which is in this
31    example $p$. Thus we can write;
32
33     \nabla^2 p = \nabla \cdot \nabla p =
34        \sum\hackscore{i}^n
35        \frac{\partial^2 p}{\partial x^2\hackscore{i}}
36     \label{eqn:laplacian}
37
38    For the two dimensional case in Cartesian coordinates \refEq{eqn:laplacian}
39    becomes;
40
41     \nabla^2 p = \frac{\partial^2 p}{\partial x^2}
42               + \frac{\partial^2 p}{\partial y^2}
43
44
45  \section{Numerical Solution Stability}  In \esc the Laplacian is calculated using the divergence representation and the
46  Unfortunately, the wave equation is belongs to a class of equations called  intrinsic functions \textit{grad()} and \textit{trace()}. The fucntion
47  \textbf{stiff} PDEs. This types of equations can be difficult to solve  \textit{grad{}} will return the spatial gradients of an object.
48  numerically as they tend to oscilate about the exact solution and can  For a rank 0 solution, this is of the form;
49  eventually blow up. To counter this problem, explicitly stable schemes like the
50  backwards Euler method are required. In terms of the wave equation, the   \nabla p = \left[
51  analytical wave must not propagate faster than the numerical wave is able to,         \frac{\partial p}{\partial x \hackscore{0}},
52  and in general, needs to be much slower than the numerical wave.         \frac{\partial p}{\partial x \hackscore{1}}
53                      \right]
54    \label{eqn:grad}
55
56    Larger ranked solution objects will return gradient tensors. For example, a
57    pressure field which acts in the directions $p \hackscore{0}$ and $p 58 \hackscore{1}$ would return;
59
60      \nabla p = \begin{bmatrix}
61           \frac{\partial p \hackscore{0}}{\partial x \hackscore{0}} &
62            \frac{\partial p \hackscore{1}}{\partial x \hackscore{0}} \\
63          \frac{\partial p \hackscore{0}}{\partial x \hackscore{1}} &
64            \frac{\partial p \hackscore{1}}{\partial x \hackscore{1}}
65                      \end{bmatrix}
66    \label{eqn:gradrank1}
67
68
69    \refEq{eqn:grad} corresponds to the Linear PDE general form value
70    $X$. Notice however that the general form contains the term $X 71 \hackscore{i,j}$\footnote{This is the first derivative in the $j^{th}$
72    direction for the $i^{th}$ component of the solution.},
73    hence for a rank 0 object there is no need to do more than calculate the
74    gradient and submit it to the solver. In the case of the rank 1 or greater
75    object, it is nesscary to calculate the trace also. This is the sum of the
76    diagonal in \refeq{eqn:gradrank1}.
77
78    Thus when solving for equations containing the Laplacian one of two things must
79    be completed. If the object \verb p   is less than rank 1 the gradient is
80    calculated via;
81    \begin{python}
82     gradient=grad(p)
83    \end{python}
84    and if the object is greater thank or equal to a rank 1 tensor, the trace of
85    the gradient is calculated.
86    \begin{python}
87     gradient=trace(grad(p))
88    \end{python}
89    These valuse can then be submitted to the PDE solver via the general form term
90    $X$. The Laplacian is then computed in the solution process by taking the
91    divergence of $X$.
92
93    Note, if you are unsure about the rank of your tensor, the \textit{getRank}
94    command will return the rank of the PDE object.
95    \begin{python}
96     rank = p.getRank()
97    \end{python}
98
99
100    \section{Numerical Solution Stability} \label{sec:nsstab}
101    Unfortunately, the wave equation belongs to a class of equations called
102    \textbf{stiff} PDEs. These types of equations can be difficult to solve
103    numerically as they tend to oscilate about the exact solution which can
104    eventually lead to a catastrophic failure in the solution. To counter this
105    problem, explicitly stable schemes like
106    the backwards Euler method are required. There are two variables which must be
107    considered for stability when numerically trying to solve the wave equation.
108    For linear media, the two variables are related via;
109     \label{eqn:freqvel}
110    f=\frac{v}{\lambda}
111
112    The velocity $v$ that a wave travels in a medium is an important variable. For
113    stability the analytical wave must not propagate faster than the numerical wave
114    is able to, and in general, needs to be much slower than the numerical wave.
115  For example, a line 100m long is discretised into 1m intervals or 101 nodes. If  For example, a line 100m long is discretised into 1m intervals or 101 nodes. If
116  a wave enters with a propagation velocity of 100m/s then the travel time for  a wave enters with a propagation velocity of 100m/s then the travel time for
117  the wave between each node will be 0.01 seconds. The time step, must therefore  the wave between each node will be 0.01 seconds. The time step, must therefore
118  be significantly less than this. Of the order $10E-4$ would be appropriate.  be significantly less than this. Of the order $10E-4$ would be appropriate.
119
120  This requirement for very small step sizes makes stiff equations difficult to  The wave frequency content also plays a part in numerical stability. The
121    nyquist-sampling theorem states that a signals bandwidth content will be
122    accurately represented when an equispaced sampling rate $f \hackscore{n}$ is
123    equal to or greater than twice the maximum frequency of the signal
124    $f\hackscore{s}$, or;
125     \label{eqn:samptheorem}
126     f\hackscore{n} \geqslant f\hackscore{s}
127
128    For example a 50Hz signal will require a sampling rate greater than 100Hz or
129    one sample every 0.01 seconds. The wave equation relies on a spatial frequency,
130    thus the sampling theorem in this case applies to the solution mesh spacing. In
131    this way, the frequency content of the input signal directly affects the time
132    discretisation of the problem.
133
134    To accurately model the wave equation with high resolutions and velocities
135    means that very fine spatial and time discretisation is necessary for most
136    problems.
137    This requirement makes the wave equation arduous to
138  solve numerically due to the large number of time iterations required in each  solve numerically due to the large number of time iterations required in each
139  solution. Models with very high velocities and fine meshes will be the worst  solution. Models with very high velocities and frequencies will be the worst
140  affected by this problem.  affected by this problem.
141

142  \section{Displacement Solution}  \section{Displacement Solution}
143  \sslist{example07a.py}  \sslist{example07a.py}
144
# Line 59  second derivative; Line 151  second derivative;
151  substituting \refEq{eqn:centdiff} for $\frac{\partial ^2 p }{\partial t ^2}$  substituting \refEq{eqn:centdiff} for $\frac{\partial ^2 p }{\partial t ^2}$
152  in \refEq{eqn:acswave};  in \refEq{eqn:acswave};
153
154   \nabla ^2 p - \frac{1}{c^2h^2} \left[p_{(t+1)} - 2p_{(t)} + p_{(t-1)} \right]   \nabla ^2 p - \frac{1}{c^2h^2} \left[p\hackscore{(t+1)} - 2p\hackscore{(t)} +
155    p\hackscore{(t-1)} \right]
156  = 0  = 0
157  \label{eqn:waveu}  \label{eqn:waveu}
158
159  Rearranging for $p_{(t+1)}$;  Rearranging for $p_{(t+1)}$;
160
161   p_{(t+1)} = c^2 h^2 \nabla ^2 p_{(t)} +2p_{(t)} - p_{(t-1)}   p\hackscore{(t+1)} = c^2 h^2 \nabla ^2 p\hackscore{(t)} +2p\hackscore{(t)} -
162    p\hackscore{(t-1)}
163
164  this can be compared with the general form of the \modLPDE module and it  this can be compared with the general form of the \modLPDE module and it
165  becomes clear that $D=1$, $X=-c^2 h^2 \nabla ^2 p_{(t)}$ and $Y=2p_{(t)} - becomes clear that$D=1$,$X\hackscore{i,j}=-c^2 h^2 \nabla ^2 p_{(t)}$and 166 p_{(t-1)}$.  $Y=2p_{(t)} - p_{(t-1)}$.
167
168    The solution script is similar to others that we have created in previous
169    chapters. The general steps are;
170    \begin{enumerate}
171     \item The necessary libraries must be imported.
172     \item The domain needs to be defined.
173     \item The time iteration and control parameters need to be defined.
174     \item The PDE is initialised with source and boundary conditions.
175     \item The time loop is started and the PDE is solved at consecutive time steps.
176     \item All or select solutions are saved to file for visualisation lated on.
177    \end{enumerate}
178
179    Parts of the script which warrant more attention are the definition of the
180    source, visualising the source, the solution time loop and the VTK data export.
181
182    \subsection{Pressure Sources}
183    As the pressure is a scalar, one need only define the pressure for two
184    time steps prior to the start of the solution loop. Two known solutions are
185    required because the wave equation contains a double partial derivative with
186    respect to time. This is often a good opportunity to introduce a source to the
187    solution. This model has the source located at it's centre. The source should
188    be smooth and cover a number of samples to satisfy the frequency stability
189    criterion. Small sources will generate high frequency signals. Here, the source
190    is defined by a cosine function.
191    \begin{python}
192    U0=0.01 # amplitude of point source
193    xc=[500,500] #location of point source
194    # define small radius around point xc
195    src_radius = 30
196    # for first two time steps
197    u=U0*(cos(length(x-xc)*3.1415/src_radius)+1)*\
198        whereNegative(length(x-xc)-src_radius)
199    u_m1=u
200    \end{python}
201    When using a rectangular domain
202
203    \subsection{Visualising the Source}
204    There are two options for visualising the source. The first is to export the
205    initial conditions of the model to VTK, which can be interpreted as a scalar
206    suface in mayavi. The second is to take a cross section of the model.
207
208    For the later, we will require the \textit{Locator} function.
209    First \verb Locator  must be imported;
210    \begin{python}
211     from esys.escript.pdetools import Locator
212    \end{python}
213    The function can then be used on the domain to locate the nearest domain node
214    to the point or points of interest.
215
216    It is now necessary to build a list of $(x,y)$ locations that specify where are
217    model slice will go. This is easily implemeted with a loop;
218    \begin{python}
219    cut_loc=[]
220    src_cut=[]
221    for i in range(ndx/2-ndx/10,ndx/2+ndx/10):
222        cut_loc.append(xstep*i)
223        src_cut.append([xstep*i,xc[1]])
224    \end{python}
225    We then submit the output to \verb Locator  and finally return the appropriate
226    values using the \verb getValue  function.
227    \begin{python}
228     src=Locator(mydomain,src_cut)
229    src_cut=src.getValue(u)
230    \end{python}
231    It is then a trivial task to plot and save the output using \mpl.
232    \begin{python}
233     pl.plot(cut_loc,src_cut)
234    pl.axis([xc[0]-src_radius*3,xc[0]+src_radius*3,0.,2*U0])
235    pl.savefig(os.path.join(savepath,"source_line.png"))
236    \end{python}
237    \begin{figure}[h]
238     \centering
239    FIXME PLEASE!
240    % \includegraphics[width=6in]{figures/sourceline.png}
241     \caption{Cross section of the source function.}
242     \label{fig:cxsource}
243    \end{figure}
244
245
246    \subsection{Point Monitoring}
247    In the more general case where the solution mesh is irregular or specific
248    locations need to be monitored, it is simple enough to use the \textit{Locator}
249    function.
250    \begin{python}
251     rec=Locator(mydomain,[250.,250.])
252    \end{python}
253     When the solution \verb u  is updated we can extract the value at that point
254    via;
255    \begin{python}
256     u_rec=rec.getValue(u)
257    \end{python}
258    For consecutive time steps one can record the values from \verb u_rec  in an
259    array initialised as \verb u_rec0=[]  with;
260    \begin{python}
261      u_rec0.append(rec.getValue(u))
262    \end{python}
263
264    It can be useful to monitor the value at a single or multiple individual points
265    in the model during the modelling process. This is done using
266    the \verb Locator  function.
267
268
269  \section{Acceleration Solution}  \section{Acceleration Solution}
270  \sslist{example07b.py}  \sslist{example07b.py}
271
272  An alternative method is to solve for the acceleration $\frac{\partial ^2 An alternative method is to solve for the acceleration$\frac{\partial ^2
273  p}{\partial t^2}$directly, and derive the the displacement solution from the p}{\partial t^2}$ directly, and derive the displacement solution from the
274  PDE solution. \refEq{eqn:waveu} is thus modified;  PDE solution. \refEq{eqn:waveu} is thus modified;
275
276    \nabla ^2 p - \frac{1}{c^2} a = 0    \nabla ^2 p - \frac{1}{c^2} a = 0
277  \label{eqn:wavea}  \label{eqn:wavea}
278
279  and can be solved directly with $Y=0$ and $X=-c^2 \nabla ^2 p_{(t)}$.  and can be solved directly with $Y=0$ and $X=-c^2 \nabla ^2 p\hackscore{(t)}$.
280  After each iteration the displacement is re-evaluated via;  After each iteration the displacement is re-evaluated via;
281
282   p_{(t+1)}=2p_{(t)} - p_{(t-1)} + h^2a   p\hackscore{(t+1)}=2p\hackscore{(t)} - p\hackscore{(t-1)} + h^2a
283
284
285    \subsection{Lumping}
286    For \esc, the acceleration solution is prefered as it allows the use of matrix
287    lumping. Lumping or mass lumping as it is sometimes known, is the process of
288    aggressively approximating the density elements of a mass matrix into the main
289    diagonal. The use of Lumping is motivaed by the simplicity of diagonal matrix
290     inversion. As a result, Lumping can significantly reduce the computational
291    requirements of a problem. Care should be taken however, as this
292    function can only be used when the $A$, $B$ and $C$ coefficients of the
293    general form are zero.
294
295    To turn lumping on in \esc one can use the command;
296    \begin{python}
297     mypde.getSolverOptions().setSolverMethod(mypde.getSolverOptions().LUMPING)
298    \end{python}
299    It is also possible to check if lumping is set using;
300    \begin{python}
301      print mypde.isUsingLumping()
302    \end{python}
303
304    \section{Stability Investigation}
305    It is now prudent to investigate the stability limitations of this problem.
306    First, we let the frequency content of the source be very small. If we define
307    the source as a cosine input, than the wavlength of the input is equal to the
308    radius of the source. Let this value be 5 meters. Now, if the maximum velocity
309    of the model is $c=380.0ms^{-1}$ then the source
310    frequency is $f\hackscore{r} = \frac{380.0}{5} = 76.0 Hz$. This is a worst case
311    scenario with a small source and the models maximum velocity.
312
313    Furthermore, we know from \refSec{sec:nsstab}, that the spatial sampling
314    frequency must be at least twice this value to ensure stability. If we assume
315    the model mesh is a square equispaced grid,
316    then the sampling interval is the side length divided by the number of samples,
317    given by $\Delta x = \frac{1000.0m}{400} = 2.5m$ and the maximum sampling
318    frequency capable at this interval is
319    $f\hackscore{s}=\frac{380.0ms^{-1}}{2.5m}=152Hz$ this is just equal to the
320    required rate satisfying \refeq{eqn:samptheorem}.
321
322    \reffig{fig:ex07sampth} depicts three examples where the grid has been
323    undersampled, sampled correctly, and over sampled. The grids used had
324    200, 400 and 800 nodes per side respectively. Obviously, the oversampled grid
325    retains the best resolution of the modelled wave.
326
327    The time step required for each of these examples is simply calculated from
328    the propagation requirement. For a maximum velocity of $380.0ms^{-1}$,
329    \begin{subequations}
330
331      \Delta t \leq \frac{1000.0m}{200} \frac{1}{380.0} = 0.013s
332
333
334      \Delta t \leq \frac{1000.0m}{400} \frac{1}{380.0} = 0.0065s
335
336
337      \Delta t \leq \frac{1000.0m}{800} \frac{1}{380.0} = 0.0032s
338
339    \end{subequations}
340    We can see, that for each doubling of the number of nodes in the mesh, we halve
341    the timestep. To illustrate the impact this has, consider our model. If the
342    source is placed at the center, it is $500m$ from the nearest boundary. With a
343    velocity of $380.0ms^{-1}$ it will take $\approx1.3s$ for the wavefront to
344    reach that boundary. In each case, this equates to $100$,  $200$ and $400$ time
345    steps. This is again, only a best case scenario, for true stability these time
346    values may need to be halved and possibly havled again.
347
348    \begin{figure}[ht]
349    \centering
350    \subfigure[Undersampled Example]{
351    \includegraphics[width=3in]{figures/ex07usamp.png}
352    \label{fig:ex07usamp}
353    }
354    \subfigure[Just sampled Example]{
355    \includegraphics[width=3in]{figures/ex07jsamp.png}
356    \label{fig:ex07jsamp}
357    }
358    \subfigure[Over sampled Example]{
359    \includegraphics[width=3in]{figures/ex07nsamp.png}
360    \label{fig:ex07nsamp}
361    }
362    \label{fig:ex07sampth}
363    \caption{Sampling Theorem example for stability
364    investigation.}
365    \end{figure}
366
367
368
369
370

Legend:
 Removed from v.3003 changed lines Added in v.3054

 ViewVC Help Powered by ViewVC 1.1.26