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

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

Parent Directory Parent Directory | Revision Log Revision Log | View Patch Patch

revision 3004 by ahallam, Wed Apr 14 05:27:30 2010 UTC revision 3054 by ahallam, Wed Jun 30 02:22:25 2010 UTC
# Line 30  The Laplacian opperator which can be wri Line 30  The Laplacian opperator which can be wri
30  calculated via the divergence of the gradient of the object, which is in this  calculated via the divergence of the gradient of the object, which is in this
31  example $p$. Thus we can write;  example $p$. Thus we can write;
32  \begin{equation}  \begin{equation}
33   \nabla^2 p = \nabla \cdot \nabla p = \frac{\partial^2 p}{\partial   \nabla^2 p = \nabla \cdot \nabla p =
34  x^2\hackscore{i}}      \sum\hackscore{i}^n
35        \frac{\partial^2 p}{\partial x^2\hackscore{i}}
36   \label{eqn:laplacian}   \label{eqn:laplacian}
37  \end{equation}  \end{equation}
38  For the two dimensional case in Cartesian coordinates \refEq{eqn:laplacian}  For the two dimensional case in Cartesian coordinates \refEq{eqn:laplacian}
# Line 66  pressure field which acts in the directi Line 67  pressure field which acts in the directi
67  \end{equation}  \end{equation}
68    
69  \refEq{eqn:grad} corresponds to the Linear PDE general form value  \refEq{eqn:grad} corresponds to the Linear PDE general form value
70  $X$. Notice however that the gernal form contains the term $X \hackscore{i,j}$,  $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  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  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  object, it is nesscary to calculate the trace also. This is the sum of the
# Line 75  diagonal in \refeq{eqn:gradrank1}. Line 78  diagonal in \refeq{eqn:gradrank1}.
78  Thus when solving for equations containing the Laplacian one of two things must  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  be completed. If the object \verb p   is less than rank 1 the gradient is
80  calculated via;  calculated via;
81  \begin{verbatim}  \begin{python}
82   gradient=grad(p)   gradient=grad(p)
83  \end{verbatim}  \end{python}
84  and if the object is greater thank or equal to a rank 1 tensor, the trace of  and if the object is greater thank or equal to a rank 1 tensor, the trace of
85  the gradient is calculated.  the gradient is calculated.
86  \begin{verbatim}  \begin{python}
87   gradient=trace(grad(p))   gradient=trace(grad(p))
88  \end{verbatim}  \end{python}
   
89  These valuse can then be submitted to the PDE solver via the general form term  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  $X$. The Laplacian is then computed in the solution process by taking the
91  divergence of $X$.  divergence of $X$.
92    
93  \section{Numerical Solution Stability}  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  Unfortunately, the wave equation belongs to a class of equations called
102  \textbf{stiff} PDEs. These types of equations can be difficult to solve  \textbf{stiff} PDEs. These types of equations can be difficult to solve
103  numerically as they tend to oscilate about the exact solution and can  numerically as they tend to oscilate about the exact solution which can
104  eventually fail. To counter this problem, explicitly stable schemes like  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  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.  considered for stability when numerically trying to solve the wave equation.
108    For linear media, the two variables are related via;
109  \begin{equation} \label{eqn:freqvel}  \begin{equation} \label{eqn:freqvel}
110  f=\frac{v}{\lambda}  f=\frac{v}{\lambda}
111  \end{equation}  \end{equation}
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  Velocity is one of these variables. For stability the  is able to, and in general, needs to be much slower than the numerical wave.
 analytical wave must not propagate faster than the numerical wave 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
# Line 129  problems. Line 137  problems.
137  This requirement makes the wave equation arduous to  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 frequencies will be the worst  solution. Models with very high velocities and frequencies will be the worst
140  effected by this problem.  affected by this problem.
141    
142  \section{Displacement Solution}  \section{Displacement Solution}
143  \sslist{example07a.py}  \sslist{example07a.py}
# Line 157  this can be compared with the general fo Line 165  this can be compared with the general fo
165  becomes clear that $D=1$, $X\hackscore{i,j}=-c^2 h^2 \nabla ^2 p_{(t)}$ and  becomes clear that $D=1$, $X\hackscore{i,j}=-c^2 h^2 \nabla ^2 p_{(t)}$ and
166  $Y=2p_{(t)} - p_{(t-1)}$.  $Y=2p_{(t)} - p_{(t-1)}$.
167    
168  The solution script is similar to other that we have created in previous  The solution script is similar to others that we have created in previous
169  chapters. The general steps are;  chapters. The general steps are;
170  \begin{enumerate}  \begin{enumerate}
171   \item The necessary libraries must be imported.   \item The necessary libraries must be imported.
# Line 180  solution. This model has the source loca Line 188  solution. This model has the source loca
188  be smooth and cover a number of samples to satisfy the frequency stability  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  criterion. Small sources will generate high frequency signals. Here, the source
190  is defined by a cosine function.  is defined by a cosine function.
191  \begin{verbatim}  \begin{python}
192  U0=0.01 # amplitude of point source  U0=0.01 # amplitude of point source
193  xc=[500,500] #location of point source  xc=[500,500] #location of point source
194  # define small radius around point xc  # define small radius around point xc
195  src_radius = 30  src_radius = 30
196  # for first two time steps  # for first two time steps
197  u=U0*(cos(length(x-xc)*3.1415/src_radius)+1)*whereNegative(length(x-xc)-src_radi  u=U0*(cos(length(x-xc)*3.1415/src_radius)+1)*\
198  us)      whereNegative(length(x-xc)-src_radius)
199  u_m1=u  u_m1=u
200  \end{verbatim}  \end{python}
201  When using a rectangular domain  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  \begin{equation}  \begin{equation}
276    \nabla ^2 p - \frac{1}{c^2} a = 0    \nabla ^2 p - \frac{1}{c^2} a = 0
# Line 208  After each iteration the displacement is Line 282  After each iteration the displacement is
282   p\hackscore{(t+1)}=2p\hackscore{(t)} - p\hackscore{(t-1)} + h^2a   p\hackscore{(t+1)}=2p\hackscore{(t)} - p\hackscore{(t-1)} + h^2a
283  \end{equation}  \end{equation}
284    
285    \subsection{Lumping}
286  For \esc, the acceleration solution is prefered as it allows the use of matrix  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  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  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  diagonal. The use of Lumping is motivaed by the simplicity of diagonal matrix
290   inversion. As a result, Lumping can significantly reduce the computational   inversion. As a result, Lumping can significantly reduce the computational
291  requirements of a problem.  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;  To turn lumping on in \esc one can use the command;
296  \begin{verbatim}  \begin{python}
297   mypde.getSolverOptions().setSolverMethod(mypde.getSolverOptions().LUMPING)   mypde.getSolverOptions().setSolverMethod(mypde.getSolverOptions().LUMPING)
298  \end{verbatim}  \end{python}
299  It is also possible to check if lumping is set using;  It is also possible to check if lumping is set using;
300  \begin{verbatim}  \begin{python}
301    print mypde.isUsingLumping()    print mypde.isUsingLumping()
302  \end{verbatim}  \end{python}
303    
304  \section{Stability Investigation}  \section{Stability Investigation}
305  It is now prudent to investigate the stability limitations of this problem.  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 the radius  First, we let the frequency content of the source be very small. If we define
307  of the source which equals the wavelength is 5 meters, than the frequency is  the source as a cosine input, than the wavlength of the input is equal to the
308  the inverse of the wavelength  The velocity is $c=380.0ms^{-1}$ thus the source  radius of the source. Let this value be 5 meters. Now, if the maximum velocity
309  frequency is $f\hackscore{r} = \frac{380.0}{5} = 76.0 Hz$. The sampling  of the model is $c=380.0ms^{-1}$ then the source
310  frequency must be at least twice this. Assuming a rectangular equispaced grid,  frequency is $f\hackscore{r} = \frac{380.0}{5} = 76.0 Hz$. This is a worst case
311  the sampling interval is $\Delta x = \frac{1000.0}{400} = 2.5$ and the sampling  scenario with a small source and the models maximum velocity.
312  frequency $f\hackscore{s}=\frac{380.0}{2.5}=152$ this is just equal to the  
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}.  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     \begin{equation}
331      \Delta t \leq \frac{1000.0m}{200} \frac{1}{380.0} = 0.013s
332     \end{equation}
333     \begin{equation}
334      \Delta t \leq \frac{1000.0m}{400} \frac{1}{380.0} = 0.0065s
335     \end{equation}
336     \begin{equation}
337      \Delta t \leq \frac{1000.0m}{800} \frac{1}{380.0} = 0.0032s
338     \end{equation}
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.3004  
changed lines
  Added in v.3054
  ViewVC Help
Powered by ViewVC 1.1.26