/[escript]/trunk/doc/cookbook/example07.tex
ViewVC logotype

Annotation of /trunk/doc/cookbook/example07.tex

Parent Directory Parent Directory | Revision Log Revision Log


Revision 3025 - (hide annotations)
Thu May 6 01:20:46 2010 UTC (10 years, 2 months ago) by ahallam
File MIME type: application/x-tex
File size: 14624 byte(s)
Updates to wave equation examples. Pressure wave examples should be completed now. Working on seismic wave examples. Cookbook chapter for pressure wave examples also added.
1 ahallam 3003
2     %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
3     %
4     % Copyright (c) 2003-2010 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     %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
13    
14    
15    
16     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
18 ahallam 3004 or longitudinal style wave motion is not possible. An obvious example is sound
19 ahallam 3003 waves.
20    
21 ahallam 3004 The acoustic wave equation is defined as;
22 ahallam 3003 \begin{equation}
23     \nabla ^2 p - \frac{1}{c^2} \frac{\partial ^2 p}{\partial t^2} = 0
24     \label{eqn:acswave}
25     \end{equation}
26 ahallam 3004 where $p$ is the pressure, $t$ is the time and $c$ is the wave velocity.
27 ahallam 3003
28 ahallam 3004 \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     \begin{equation}
33     \nabla^2 p = \nabla \cdot \nabla p = \frac{\partial^2 p}{\partial
34     x^2\hackscore{i}}
35     \label{eqn:laplacian}
36     \end{equation}
37     For the two dimensional case in Cartesian coordinates \refEq{eqn:laplacian}
38     becomes;
39     \begin{equation}
40     \nabla^2 p = \frac{\partial^2 p}{\partial x^2}
41     + \frac{\partial^2 p}{\partial y^2}
42     \end{equation}
43 ahallam 3003
44 ahallam 3004 In \esc the Laplacian is calculated using the divergence representation and the
45     intrinsic functions \textit{grad()} and \textit{trace()}. The fucntion
46     \textit{grad{}} will return the spatial gradients of an object.
47     For a rank 0 solution, this is of the form;
48     \begin{equation}
49     \nabla p = \left[
50     \frac{\partial p}{\partial x \hackscore{0}},
51     \frac{\partial p}{\partial x \hackscore{1}}
52     \right]
53     \label{eqn:grad}
54     \end{equation}
55     Larger ranked solution objects will return gradient tensors. For example, a
56     pressure field which acts in the directions $p \hackscore{0}$ and $p
57     \hackscore{1}$ would return;
58     \begin{equation}
59     \nabla p = \begin{bmatrix}
60     \frac{\partial p \hackscore{0}}{\partial x \hackscore{0}} &
61     \frac{\partial p \hackscore{1}}{\partial x \hackscore{0}} \\
62     \frac{\partial p \hackscore{0}}{\partial x \hackscore{1}} &
63     \frac{\partial p \hackscore{1}}{\partial x \hackscore{1}}
64     \end{bmatrix}
65     \label{eqn:gradrank1}
66     \end{equation}
67 ahallam 3003
68 ahallam 3004 \refEq{eqn:grad} corresponds to the Linear PDE general form value
69     $X$. Notice however that the gernal form contains the term $X \hackscore{i,j}$,
70     hence for a rank 0 object there is no need to do more than calculate the
71     gradient and submit it to the solver. In the case of the rank 1 or greater
72     object, it is nesscary to calculate the trace also. This is the sum of the
73     diagonal in \refeq{eqn:gradrank1}.
74    
75     Thus when solving for equations containing the Laplacian one of two things must
76     be completed. If the object \verb p is less than rank 1 the gradient is
77     calculated via;
78 ahallam 3025 \begin{python}
79 ahallam 3004 gradient=grad(p)
80 ahallam 3025 \end{python}
81 ahallam 3004 and if the object is greater thank or equal to a rank 1 tensor, the trace of
82     the gradient is calculated.
83 ahallam 3025 \begin{python}
84 ahallam 3004 gradient=trace(grad(p))
85 ahallam 3025 \end{python}
86 ahallam 3004 These valuse can then be submitted to the PDE solver via the general form term
87     $X$. The Laplacian is then computed in the solution process by taking the
88     divergence of $X$.
89    
90 ahallam 3025 Note, if you are unsure about the rank of your tensor, the \textit{getRank}
91     command will return the rank of the PDE object.
92     \begin{python}
93     rank = p.getRank()
94     \end{python}
95    
96    
97     \section{Numerical Solution Stability} \label{sec:nsstab}
98 ahallam 3004 Unfortunately, the wave equation belongs to a class of equations called
99     \textbf{stiff} PDEs. These types of equations can be difficult to solve
100 ahallam 3025 numerically as they tend to oscilate about the exact solution which can
101     eventually lead to a catastrophic failure in the solution. To counter this
102     problem, explicitly stable schemes like
103 ahallam 3004 the backwards Euler method are required. There are two variables which must be
104     considered for stability when numerically trying to solve the wave equation.
105 ahallam 3025 For linear media, the two variables are related via;
106 ahallam 3004 \begin{equation} \label{eqn:freqvel}
107     f=\frac{v}{\lambda}
108     \end{equation}
109 ahallam 3025 The velocity $v$ that a wave travels in a medium is an important variable. For
110     stability the analytical wave must not propagate faster than the numerical wave
111     is able to, and in general, needs to be much slower than the numerical wave.
112 ahallam 3003 For example, a line 100m long is discretised into 1m intervals or 101 nodes. If
113     a wave enters with a propagation velocity of 100m/s then the travel time for
114     the wave between each node will be 0.01 seconds. The time step, must therefore
115     be significantly less than this. Of the order $10E-4$ would be appropriate.
116    
117 ahallam 3004 The wave frequency content also plays a part in numerical stability. The
118     nyquist-sampling theorem states that a signals bandwidth content will be
119     accurately represented when an equispaced sampling rate $f \hackscore{n}$ is
120     equal to or greater than twice the maximum frequency of the signal
121     $f\hackscore{s}$, or;
122     \begin{equation} \label{eqn:samptheorem}
123     f\hackscore{n} \geqslant f\hackscore{s}
124     \end{equation}
125     For example a 50Hz signal will require a sampling rate greater than 100Hz or
126     one sample every 0.01 seconds. The wave equation relies on a spatial frequency,
127     thus the sampling theorem in this case applies to the solution mesh spacing. In
128     this way, the frequency content of the input signal directly affects the time
129     discretisation of the problem.
130    
131     To accurately model the wave equation with high resolutions and velocities
132     means that very fine spatial and time discretisation is necessary for most
133     problems.
134     This requirement makes the wave equation arduous to
135 ahallam 3003 solve numerically due to the large number of time iterations required in each
136 ahallam 3004 solution. Models with very high velocities and frequencies will be the worst
137     effected by this problem.
138 ahallam 3003
139     \section{Displacement Solution}
140     \sslist{example07a.py}
141    
142     We begin the solution to this PDE with the centred difference formula for the
143     second derivative;
144     \begin{equation}
145     f''(x) \approx \frac{f(x+h - 2f(x) + f(x-h)}{h^2}
146     \label{eqn:centdiff}
147     \end{equation}
148     substituting \refEq{eqn:centdiff} for $\frac{\partial ^2 p }{\partial t ^2}$
149     in \refEq{eqn:acswave};
150     \begin{equation}
151 ahallam 3004 \nabla ^2 p - \frac{1}{c^2h^2} \left[p\hackscore{(t+1)} - 2p\hackscore{(t)} +
152     p\hackscore{(t-1)} \right]
153 ahallam 3003 = 0
154     \label{eqn:waveu}
155     \end{equation}
156     Rearranging for $p_{(t+1)}$;
157     \begin{equation}
158 ahallam 3004 p\hackscore{(t+1)} = c^2 h^2 \nabla ^2 p\hackscore{(t)} +2p\hackscore{(t)} -
159     p\hackscore{(t-1)}
160 ahallam 3003 \end{equation}
161     this can be compared with the general form of the \modLPDE module and it
162 ahallam 3004 becomes clear that $D=1$, $X\hackscore{i,j}=-c^2 h^2 \nabla ^2 p_{(t)}$ and
163     $Y=2p_{(t)} - p_{(t-1)}$.
164 ahallam 3003
165 ahallam 3025 The solution script is similar to others that we have created in previous
166 ahallam 3004 chapters. The general steps are;
167     \begin{enumerate}
168     \item The necessary libraries must be imported.
169     \item The domain needs to be defined.
170     \item The time iteration and control parameters need to be defined.
171     \item The PDE is initialised with source and boundary conditions.
172     \item The time loop is started and the PDE is solved at consecutive time steps.
173     \item All or select solutions are saved to file for visualisation lated on.
174     \end{enumerate}
175    
176     Parts of the script which warrant more attention are the definition of the
177     source, visualising the source, the solution time loop and the VTK data export.
178    
179     \subsection{Pressure Sources}
180     As the pressure is a scalar, one need only define the pressure for two
181     time steps prior to the start of the solution loop. Two known solutions are
182     required because the wave equation contains a double partial derivative with
183     respect to time. This is often a good opportunity to introduce a source to the
184     solution. This model has the source located at it's centre. The source should
185     be smooth and cover a number of samples to satisfy the frequency stability
186     criterion. Small sources will generate high frequency signals. Here, the source
187     is defined by a cosine function.
188 ahallam 3025 \begin{python}
189 ahallam 3004 U0=0.01 # amplitude of point source
190     xc=[500,500] #location of point source
191     # define small radius around point xc
192     src_radius = 30
193     # for first two time steps
194 ahallam 3025 u=U0*(cos(length(x-xc)*3.1415/src_radius)+1)*\
195     whereNegative(length(x-xc)-src_radius)
196 ahallam 3004 u_m1=u
197 ahallam 3025 \end{python}
198 ahallam 3004 When using a rectangular domain
199    
200 ahallam 3025 \subsection{Visualising the Source}
201     There are two options for visualising the source. The first is to export the
202     initial conditions of the model to VTK, which can be interpreted as a scalar
203     suface in mayavi. The second is to take a cross section of the model.
204    
205     For the later, we will require the \textit{Locator} function.
206     First \verb Locator must be imported;
207     \begin{python}
208     from esys.escript.pdetools import Locator
209     \end{python}
210     The function can then be used on the domain to locate the nearest domain node
211     to the point or points of interest.
212    
213     It is now necessary to build a list of $(x,y)$ locations that specify where are
214     model slice will go. This is easily implemeted with a loop;
215     \begin{python}
216     cut_loc=[]
217     src_cut=[]
218     for i in range(ndx/2-ndx/10,ndx/2+ndx/10):
219     cut_loc.append(xstep*i)
220     src_cut.append([xstep*i,xc[1]])
221     \end{python}
222     We then submit the output to \verb Locator and finally return the appropriate
223     values using the \verb getValue function.
224     \begin{python}
225     src=Locator(mydomain,src_cut)
226     src_cut=src.getValue(u)
227     \end{python}
228     It is then a trivial task to plot and save the output using \mpl.
229     \begin{python}
230     pl.plot(cut_loc,src_cut)
231     pl.axis([xc[0]-src_radius*3,xc[0]+src_radius*3,0.,2*U0])
232     pl.savefig(os.path.join(savepath,"source_line.png"))
233     \end{python}
234     \begin{figure}[h]
235     \includegraphics[width=5in]{figures/sourceline.png}
236     \caption{Cross section of the source function.}
237     \label{fig:cxsource}
238     \end{figure}
239    
240    
241     \subsection{Point Monitoring}
242     In the more general case where the solution mesh is irregular or specific
243     locations need to be monitored, it is simple enough to use the \textit{Locator}
244     function.
245     \begin{python}
246     rec=Locator(mydomain,[250.,250.])
247     \end{python}
248     When the solution \verb u is update we can extract the value at that point
249     via;
250     \begin{python}
251     u_rec=rec.getValue(u)
252     \end{python}
253     For consecutive time steps one can record the values from \verb u_rec in an
254     array initialised as \verb u_rec0=[] with;
255     \begin{python}
256     u_rec0.append(rec.getValue(u))
257     \end{python}
258    
259     It can be useful to monitor the value at a single or multiple individual points
260     in the model during the modelling process. This is done using
261     the \verb Locator function.
262    
263    
264 ahallam 3003 \section{Acceleration Solution}
265     \sslist{example07b.py}
266    
267     An alternative method is to solve for the acceleration $\frac{\partial ^2
268     p}{\partial t^2}$ directly, and derive the the displacement solution from the
269     PDE solution. \refEq{eqn:waveu} is thus modified;
270     \begin{equation}
271     \nabla ^2 p - \frac{1}{c^2} a = 0
272     \label{eqn:wavea}
273     \end{equation}
274 ahallam 3004 and can be solved directly with $Y=0$ and $X=-c^2 \nabla ^2 p\hackscore{(t)}$.
275 ahallam 3003 After each iteration the displacement is re-evaluated via;
276     \begin{equation}
277 ahallam 3004 p\hackscore{(t+1)}=2p\hackscore{(t)} - p\hackscore{(t-1)} + h^2a
278 ahallam 3003 \end{equation}
279    
280 ahallam 3004 For \esc, the acceleration solution is prefered as it allows the use of matrix
281     lumping. Lumping or mass lumping as it is sometimes known, is the process of
282     aggressively approximating the density elements of a mass matrix into the main
283     diagonal. The use of Lumping is motivaed by the simplicity of diagonal matrix
284     inversion. As a result, Lumping can significantly reduce the computational
285     requirements of a problem.
286 ahallam 3003
287 ahallam 3004 To turn lumping on in \esc one can use the command;
288 ahallam 3025 \begin{python}
289 ahallam 3004 mypde.getSolverOptions().setSolverMethod(mypde.getSolverOptions().LUMPING)
290 ahallam 3025 \end{python}
291 ahallam 3004 It is also possible to check if lumping is set using;
292 ahallam 3025 \begin{python}
293 ahallam 3004 print mypde.isUsingLumping()
294 ahallam 3025 \end{python}
295 ahallam 3004
296     \section{Stability Investigation}
297     It is now prudent to investigate the stability limitations of this problem.
298 ahallam 3025 First, we let the frequency content of the source be very small. If we define
299     the source as a cosine input, than the wavlength of the input is equal to the
300     radius of the source. Let this value be 5 meters. Now, if the maximum velocity
301     of the model is $c=380.0ms^{-1}$ then the source
302     frequency is $f\hackscore{r} = \frac{380.0}{5} = 76.0 Hz$. This is a worst case
303     scenario with a small source and the models maximum velocity.
304    
305     Furthermore, we know from \refSec{sec:nsstab}, that the spatial sampling
306     frequency must be at least twice this value to ensure stability. If we assume
307     the model mesh is a square equispaced grid,
308     then the sampling interval is the side length divided by the number of samples,
309     given by $\Delta x = \frac{1000.0m}{400} = 2.5m$ and the maximum sampling
310     frequency capable at this interval is
311     $f\hackscore{s}=\frac{380.0ms^{-1}}{2.5m}=152Hz$ this is just equal to the
312 ahallam 3004 required rate satisfying \refeq{eqn:samptheorem}.
313    
314 ahallam 3025 \reffig{fig:ex07sampth} depicts three examples where the grid has been
315     undersampled, sampled correctly, and over sampled. The grids used had
316     200, 400 and 800 nodes per side respectively. Obviously, the oversampled grid
317     retains the best resolution of the modelled wave.
318 ahallam 3004
319 ahallam 3025 The time step required for each of these examples is simply calculated from
320     the propagation requirement. For a maximum velocity of $380.0ms^{-1}$,
321     \begin{subequations}
322     \begin{equation}
323     \Delta t \leq \frac{1000.0m}{200} \frac{1}{380.0} = 0.013s
324     \end{equation}
325     \begin{equation}
326     \Delta t \leq \frac{1000.0m}{400} \frac{1}{380.0} = 0.0065s
327     \end{equation}
328     \begin{equation}
329     \Delta t \leq \frac{1000.0m}{800} \frac{1}{380.0} = 0.0032s
330     \end{equation}
331     \end{subequations}
332     We can see, that for each doubling of the number of nodes in the mesh, we halve
333     the timestep. To illustrate the impact this has, consider our model. If the
334     source is placed at the center, it is $500m$ from the nearest boundary. With a
335     velocity of $380.0ms^{-1}$ it will take $\approx1.3s$ for the wavefront to
336     reach that boundary. In each case, this equates to $100$, $200$ and $400$ time
337     steps. This is again, only a best case scenario, for true stability these time
338     values may need to be halved and possibly havled again.
339 ahallam 3004
340 ahallam 3025 \begin{figure}[ht]
341     \centering
342     \subfigure[Undersampled Example]{
343     \includegraphics[width=3in]{figures/ex07usamp.png}
344     \label{fig:ex07usamp}
345     }
346     \subfigure[Just sampled Example]{
347     \includegraphics[width=3in]{figures/ex07jsamp.png}
348     \label{fig:ex07jsamp}
349     }
350     \subfigure[Over sampled Example]{
351     \includegraphics[width=3in]{figures/ex07nsamp.png}
352     \label{fig:ex07nsamp}
353     }
354     \label{fig:ex07sampth}
355     \caption{Sampling Theorem example for stability
356     investigation.}
357     \end{figure}
358 ahallam 3004
359 ahallam 3025
360    
361    
362    

  ViewVC Help
Powered by ViewVC 1.1.26