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

revision 3054 by ahallam, Wed Jun 30 02:22:25 2010 UTC revision 3308 by jfenwick, Tue Oct 26 03:24:54 2010 UTC
# Line 31  calculated via the divergence of the gra Line 31  calculated via the divergence of the gra
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 =   \nabla^2 p = \nabla \cdot \nabla p =
34      \sum\hackscore{i}^n      \sum_{i}^n
35      \frac{\partial^2 p}{\partial x^2\hackscore{i}}      \frac{\partial^2 p}{\partial x^2_{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 \autoref{eqn:laplacian}
39  becomes;  becomes;
40  \begin{equation}  \begin{equation}
41   \nabla^2 p = \frac{\partial^2 p}{\partial x^2}   \nabla^2 p = \frac{\partial^2 p}{\partial x^2}
# Line 48  intrinsic functions \textit{grad()} and Line 48  intrinsic functions \textit{grad()} and
48  For a rank 0 solution, this is of the form;  For a rank 0 solution, this is of the form;
49  \begin{equation}  \begin{equation}
50   \nabla p = \left[   \nabla p = \left[
51         \frac{\partial p}{\partial x \hackscore{0}},           \frac{\partial p}{\partial x _{0}},
52         \frac{\partial p}{\partial x \hackscore{1}}         \frac{\partial p}{\partial x _{1}}
53                    \right]                    \right]
54  \label{eqn:grad}  \label{eqn:grad}
55  \end{equation}  \end{equation}
56  Larger ranked solution objects will return gradient tensors. For example, a  Larger ranked solution objects will return gradient tensors. For example, a
57  pressure field which acts in the directions $p \hackscore{0}$ and $p pressure field which acts in the directions$p _{0}$and$p
58  \hackscore{1}$would return; _{1}$ would return;
59  \begin{equation}  \begin{equation}
60    \nabla p = \begin{bmatrix}    \nabla p = \begin{bmatrix}
61         \frac{\partial p \hackscore{0}}{\partial x \hackscore{0}} &         \frac{\partial p _{0}}{\partial x _{0}} &
62          \frac{\partial p \hackscore{1}}{\partial x \hackscore{0}} \\          \frac{\partial p _{1}}{\partial x _{0}} \\
63        \frac{\partial p \hackscore{0}}{\partial x \hackscore{1}} &        \frac{\partial p _{0}}{\partial x _{1}} &
64          \frac{\partial p \hackscore{1}}{\partial x \hackscore{1}}          \frac{\partial p _{1}}{\partial x _{1}}
65                    \end{bmatrix}                    \end{bmatrix}
66  \label{eqn:gradrank1}  \label{eqn:gradrank1}
67  \end{equation}  \end{equation}
68
69  \refEq{eqn:grad} corresponds to the Linear PDE general form value  \autoref{eqn:grad} corresponds to the Linear PDE general form value
70  $X$. Notice however that the general form contains the term $X$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}$_{i,j}$\footnote{This is the first derivative in the $j^{th}$
72  direction for the $i^{th}$ component of the solution.},  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
76  diagonal in \refeq{eqn:gradrank1}.  diagonal in \autoref{eqn:gradrank1}.
77
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{python}  \begin{python}
82   gradient=grad(p)  gradient=grad(p)
83  \end{python}  \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.
# Line 119  be significantly less than this. Of the Line 119  be significantly less than this. Of the
119
120  The wave frequency content also plays a part in numerical stability. The  The wave frequency content also plays a part in numerical stability. The
121  nyquist-sampling theorem states that a signals bandwidth content will be  nyquist-sampling theorem states that a signals bandwidth content will be
122  accurately represented when an equispaced sampling rate $f \hackscore{n}$ is  accurately represented when an equispaced sampling rate $f _{n}$ is
123  equal to or greater than twice the maximum frequency of the signal  equal to or greater than twice the maximum frequency of the signal
124  $f\hackscore{s}$, or;  $f_{s}$, or;
125  \begin{equation} \label{eqn:samptheorem}  \begin{equation} \label{eqn:samptheorem}
126   f\hackscore{n} \geqslant f\hackscore{s}   f_{n} \geqslant f_{s}
127  \end{equation}  \end{equation}
128  For example a 50Hz signal will require a sampling rate greater than 100Hz or  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,  one sample every 0.01 seconds. The wave equation relies on a spatial frequency,
# Line 148  second derivative; Line 148  second derivative;
148   f''(x) \approx \frac{f(x+h - 2f(x) + f(x-h)}{h^2}   f''(x) \approx \frac{f(x+h - 2f(x) + f(x-h)}{h^2}
149  \label{eqn:centdiff}  \label{eqn:centdiff}
150  \end{equation}  \end{equation}
151  substituting \refEq{eqn:centdiff} for $\frac{\partial ^2 p }{\partial t ^2}$  substituting \autoref{eqn:centdiff} for $\frac{\partial ^2 p }{\partial t ^2}$
152  in \refEq{eqn:acswave};  in \autoref{eqn:acswave};
153  \begin{equation}  \begin{equation}
154   \nabla ^2 p - \frac{1}{c^2h^2} \left[p\hackscore{(t+1)} - 2p\hackscore{(t)} +   \nabla ^2 p - \frac{1}{c^2h^2} \left[p_{(t+1)} - 2p_{(t)} +
155  p\hackscore{(t-1)} \right]  p_{(t-1)} \right]
156  = 0  = 0
157  \label{eqn:waveu}  \label{eqn:waveu}
158  \end{equation}  \end{equation}
159  Rearranging for $p_{(t+1)}$;  Rearranging for $p_{(t+1)}$;
160  \begin{equation}  \begin{equation}
161   p\hackscore{(t+1)} = c^2 h^2 \nabla ^2 p\hackscore{(t)} +2p\hackscore{(t)} -   p_{(t+1)} = c^2 h^2 \nabla ^2 p_{(t)} +2p_{(t)} -
162  p\hackscore{(t-1)}  p_{(t-1)}
163  \end{equation}  \end{equation}
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\hackscore{i,j}=-c^2 h^2 \nabla ^2 p_{(t)}$ and  becomes clear that $D=1$, $X_{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 others that we have created in previous  The solution script is similar to others that we have created in previous
# Line 206  initial conditions of the model to VTK, Line 206  initial conditions of the model to VTK,
206  suface in mayavi. The second is to take a cross section of the model.  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.  For the later, we will require the \textit{Locator} function.
209  First \verb Locator  must be imported;  First \verb!Locator! must be imported;
210  \begin{python}  \begin{python}
211   from esys.escript.pdetools import Locator   from esys.escript.pdetools import Locator
212  \end{python}  \end{python}
# Line 222  for i in range(ndx/2-ndx/10,ndx/2+ndx/10 Line 222  for i in range(ndx/2-ndx/10,ndx/2+ndx/10
222      cut_loc.append(xstep*i)      cut_loc.append(xstep*i)
223      src_cut.append([xstep*i,xc])      src_cut.append([xstep*i,xc])
224  \end{python}  \end{python}
225  We then submit the output to \verb Locator  and finally return the appropriate  We then submit the output to \verb!Locator! and finally return the appropriate
226  values using the \verb getValue  function.  values using the \verb!getValue! function.
227  \begin{python}  \begin{python}
228   src=Locator(mydomain,src_cut)   src=Locator(mydomain,src_cut)
229  src_cut=src.getValue(u)  src_cut=src.getValue(u)
# Line 236  pl.savefig(os.path.join(savepath,"source Line 236  pl.savefig(os.path.join(savepath,"source
236  \end{python}  \end{python}
237  \begin{figure}[h]  \begin{figure}[h]
238   \centering   \centering
239  FIXME PLEASE!      \includegraphics[width=6in]{figures/sourceline.png}
% \includegraphics[width=6in]{figures/sourceline.png}
240   \caption{Cross section of the source function.}   \caption{Cross section of the source function.}
241   \label{fig:cxsource}   \label{fig:cxsource}
242  \end{figure}  \end{figure}
# Line 255  via; Line 254  via;
254  \begin{python}  \begin{python}
255   u_rec=rec.getValue(u)   u_rec=rec.getValue(u)
256  \end{python}  \end{python}
257  For consecutive time steps one can record the values from \verb u_rec  in an  For consecutive time steps one can record the values from \verb!u_rec! in an
258  array initialised as \verb u_rec0=[]  with;  array initialised as \verb!u_rec0=[]! with;
259  \begin{python}  \begin{python}
260    u_rec0.append(rec.getValue(u))    u_rec0.append(rec.getValue(u))
261  \end{python}  \end{python}
262
263  It can be useful to monitor the value at a single or multiple individual points  It can be useful to monitor the value at a single or multiple individual points
264  in the model during the modelling process. This is done using  in the model during the modelling process. This is done using
265  the \verb Locator  function.  the \verb!Locator! function.
266
267
268  \section{Acceleration Solution}  \section{Acceleration Solution}
# Line 271  the \verb Locator  function. Line 270  the \verb Locator  function.
270
271  An alternative method is to solve for the acceleration $\frac{\partial ^2 An alternative method is to solve for the acceleration$\frac{\partial ^2
272  p}{\partial t^2}$directly, and derive the displacement solution from the p}{\partial t^2}$ directly, and derive the displacement solution from the
273  PDE solution. \refEq{eqn:waveu} is thus modified;  PDE solution. \autoref{eqn:waveu} is thus modified;
274  \begin{equation}  \begin{equation}
275    \nabla ^2 p - \frac{1}{c^2} a = 0    \nabla ^2 p - \frac{1}{c^2} a = 0
276  \label{eqn:wavea}  \label{eqn:wavea}
277  \end{equation}  \end{equation}
278  and can be solved directly with $Y=0$ and $X=-c^2 \nabla ^2 p\hackscore{(t)}$.  and can be solved directly with $Y=0$ and $X=-c^2 \nabla ^2 p_{(t)}$.
279  After each iteration the displacement is re-evaluated via;  After each iteration the displacement is re-evaluated via;
280  \begin{equation}  \begin{equation}
281   p\hackscore{(t+1)}=2p\hackscore{(t)} - p\hackscore{(t-1)} + h^2a   p_{(t+1)}=2p_{(t)} - p_{(t-1)} + h^2a
282  \end{equation}  \end{equation}
283
284  \subsection{Lumping}  \subsection{Lumping}
# Line 307  First, we let the frequency content of t Line 306  First, we let the frequency content of t
306  the source as a cosine input, than the wavlength of the input is equal to the  the source as a cosine input, than the wavlength of the input is equal to the
307  radius of the source. Let this value be 5 meters. Now, if the maximum velocity  radius of the source. Let this value be 5 meters. Now, if the maximum velocity
308  of the model is $c=380.0ms^{-1}$ then the source  of the model is $c=380.0ms^{-1}$ then the source
309  frequency is $f\hackscore{r} = \frac{380.0}{5} = 76.0 Hz$. This is a worst case  frequency is $f_{r} = \frac{380.0}{5} = 76.0 Hz$. This is a worst case
310  scenario with a small source and the models maximum velocity.  scenario with a small source and the models maximum velocity.
311
312  Furthermore, we know from \refSec{sec:nsstab}, that the spatial sampling  Furthermore, we know from \autoref{sec:nsstab}, that the spatial sampling
313  frequency must be at least twice this value to ensure stability. If we assume  frequency must be at least twice this value to ensure stability. If we assume
314  the model mesh is a square equispaced grid,  the model mesh is a square equispaced grid,
315  then the sampling interval is the side length divided by the number of samples,  then the sampling interval is the side length divided by the number of samples,
316  given by $\Delta x = \frac{1000.0m}{400} = 2.5m$ and the maximum sampling  given by $\Delta x = \frac{1000.0m}{400} = 2.5m$ and the maximum sampling
317  frequency capable at this interval is  frequency capable at this interval is
318  $f\hackscore{s}=\frac{380.0ms^{-1}}{2.5m}=152Hz$ this is just equal to the  $f_{s}=\frac{380.0ms^{-1}}{2.5m}=152Hz$ this is just equal to the
319  required rate satisfying \refeq{eqn:samptheorem}.  required rate satisfying \autoref{eqn:samptheorem}.
320
321  \reffig{fig:ex07sampth} depicts three examples where the grid has been  \autoref{fig:ex07sampth} depicts three examples where the grid has been
322  undersampled, sampled correctly, and over sampled. The grids used had  undersampled, sampled correctly, and over sampled. The grids used had
323  200, 400 and 800 nodes per side respectively. Obviously, the oversampled grid  200, 400 and 800 nodes per side respectively. Obviously, the oversampled grid
324  retains the best resolution of the modelled wave.  retains the best resolution of the modelled wave.

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

 ViewVC Help Powered by ViewVC 1.1.26