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

revision 3369 by jfenwick, Tue Oct 26 03:24:54 2010 UTC revision 3370 by ahallam, Sun Nov 21 23:22:25 2010 UTC
# Line 35  waves, $\mu$ is the propagation material Line 35  waves, $\mu$ is the propagation material
35  In a similar process to the previous chapter, we will use the acceleration  In a similar process to the previous chapter, we will use the acceleration
36  solution to solve this PDE. By substituting $a$ directly for  solution to solve this PDE. By substituting $a$ directly for
37  $\frac{\partial^{2}u_{i}}{\partial t^2}$ we can derive the  $\frac{\partial^{2}u_{i}}{\partial t^2}$ we can derive the
38  displacement solution. Using $a$ we can see that \autoref{eqn:wav} becomes  acceleration solution. Using $a$ we can see that \autoref{eqn:wav} becomes
39  \begin{equation} \label{eqn:wava}  \begin{equation} \label{eqn:wava}
40  \rho a_{i} - \frac{\partial  \rho a_{i} - \frac{\partial
41  \sigma_{ij}}{\partial x_{j}} = 0  \sigma_{ij}}{\partial x_{j}} = 0
42  \end{equation}  \end{equation}
43  Thus the problem will be solved for acceleration and then converted to  Thus the problem will be solved for acceleration and then converted to
44  displacement using the backwards difference approximation.  displacement using the backwards difference approximation as for the acoustic
45    example in the previous chapter.
46
47  Consider now the stress $\sigma$. One can see that the stress consists of two  Consider now the stress $\sigma$. One can see that the stress consists of two
48  distinct terms:  distinct terms:
# Line 58  trace of the displacement solution and t Line 59  trace of the displacement solution and t
59  kronecker delta function with dimensions equivalent to $u$. The second term  kronecker delta function with dimensions equivalent to $u$. The second term
60  \autoref{eqn:sigtrans} is the sum of $u$ with its own transpose. Putting these  \autoref{eqn:sigtrans} is the sum of $u$ with its own transpose. Putting these
61  facts together we see that the spatial differential of the stress is given by the  facts together we see that the spatial differential of the stress is given by the
62  gradient of $u$ and the aforementioned opperations. This value is then submitted  gradient of $u$ and the aforementioned operations. This value is then submitted
63  to the \esc PDE as $X$.  to the \esc PDE as $X$.
64  \begin{python}  \begin{python}
# Line 68  The solution is then obtained via the us Line 69  The solution is then obtained via the us
69  calculated so that the memory variables can be updated for the next time  calculated so that the memory variables can be updated for the next time
70  iteration.  iteration.
71  \begin{python}  \begin{python}
72  accel = mypde.getSolution() #get PDE solution for accelleration  accel = mypde.getSolution() #get PDE solution for acceleration
73  u_p1=(2.*u-u_m1)+h*h*accel #calculate displacement  u_p1=(2.*u-u_m1)+h*h*accel #calculate displacement
74  u_m1=u; u=u_p1 # shift values by 1  u_m1=u; u=u_p1 # shift values by 1
75  \end{python}  \end{python}
# Line 98  then \verb!t! the solution is saved and Line 99  then \verb!t! the solution is saved and
99  limits the number of outputs and increases the speed of the solver.  limits the number of outputs and increases the speed of the solver.
100
102  The wave equation solution can be quite demanding on cpu time. Enhancements can  The wave equation solution can be quite demanding on CPU time. Enhancements can
104  require any modification to the solution script and only comes into play when  require any modification to the solution script and only comes into play when
105  escript is called from the shell. To use multiple threads \esc is called using  \esc is called from the shell. To use multiple threads \esc is called using
106  the \verb!-t! option with an interger argument for the number of threads  the \verb!-t! option with an integer argument for the number of threads
107  required. For example  required. For example
108  \begin{verbatim}  \begin{verbatim}
109  $escript -t 4 example08a.py$escript -t 4 example08a.py
110  \end{verbatim}  \end{verbatim}
111  would call the script in this section and solve it using 4 threads.  would call the script in this section and solve it using 4 threads.
112
113  The computation times on an increasing number of cores is outlines in  The computation times on an increasing number of cores is outlined in
114  \autoref{tab:wpcores}.  \autoref{tab:wpcores}.
115
116  \begin{table}[ht]  \begin{table}[ht]
# Line 135  Number of Cores & Time (s) \\ Line 136  Number of Cores & Time (s) \\
136  \section{Vector source on the boundary}  \section{Vector source on the boundary}
137  \sslist{example08b.py}  \sslist{example08b.py}
138  For this particular example, we will introduce the source by applying a  For this particular example, we will introduce the source by applying a
139  displacment to the boundary during the initial time steps. The source will again  displacement to the boundary during the initial time steps. The source will
140  be  again be
141  a radially propagating wave but due to the vector nature of the PDE used, a  a radially propagating wave but due to the vector nature of the PDE used, a
142  direction will need to be applied to the source.  direction will need to be applied to the source.
143
# Line 159  y=y*src_dir Line 160  y=y*src_dir
160  \end{python}  \end{python}
161  where \verb xc  is the source point on the boundary of the model. Note that  where \verb xc  is the source point on the boundary of the model. Note that
162  because the source is specifically located on the boundary, we have used the  because the source is specifically located on the boundary, we have used the
163  \verb!FunctionOnBoundary! call to ensure the nodes are located upon the  \verb!FunctionOnBoundary! call to ensure the nodes used to define the source
164  boundary only. These boundary nodes are passed to source as \verb!xb!. The  are also located upon the boundary. These boundary nodes are passed to
165  source direction is then defined as an $(x,y)$ array and multiplied by the  source as \verb!xb!. The source direction is then defined as an $(x,y)$ array and multiplied by the
166  source function. The directional array must have a magnitude of $\left| 1 source function. The directional array must have a magnitude of$\left| 1
167  \right| $otherwise the amplitude of the source will become modified. For this \right|$ otherwise the amplitude of the source will become modified. For this
168  example, the source is directed in the $-y$ direction.  example, the source is directed in the $-y$ direction.
# Line 221  frequency range and thus mitigate aliasi Line 222  frequency range and thus mitigate aliasi
222  It is quite  It is quite
223  simple to implement a source which is smooth in time. In addition to the  simple to implement a source which is smooth in time. In addition to the
224  original source function the only extra requirement is a time function. For  original source function the only extra requirement is a time function. For
225  this example the time variant source will be the derivative of a gausian curve  this example the time variant source will be the derivative of a Gaussian curve
226  defined by the required dominant frequency (\autoref{fig:tvsource}).  defined by the required dominant frequency (\autoref{fig:tvsource}).
227  \begin{python}  \begin{python}
228  #Creating the time function of the source.  #Creating the time function of the source.
# Line 264  u_m1=u Line 265  u_m1=u
265  \end{python}  \end{python}
266
267  Finally, for the length of the source, we are required to update each new  Finally, for the length of the source, we are required to update each new
268  solution in the itterative section of the solver. This is done via;  solution in the iterative section of the solver. This is done via;
269  \begin{python}  \begin{python}
270  # increment loop values  # increment loop values
271  t=t+h; n=n+1  t=t+h; n=n+1
# Line 277  boundary Line 278  boundary
278
279  \section{Absorbing Boundary Conditions}  \section{Absorbing Boundary Conditions}
280  To mitigate the effect of the boundary on the model, absorbing boundary  To mitigate the effect of the boundary on the model, absorbing boundary
281  conditions can be introduced. These conditions effectively dampen the wave  conditions can be introduced. These conditions effectively dampen the wave energy
282  energy as they approach the bounday and thus prevent that energy from being  as they approach the boundary and thus prevent that energy from being reflected.
283  reflected. This type of approach is used typically when a model only represents  This type of approach is typically used when a model is shrunk to decrease
284  a small portion of the entire model, which in reality may have infinite bounds.  computational requirements. In practise this applies to almost all models,
285  It is inpractical to calculate the solution for an infinite model and thus ABCs  especially in earth sciences where the entire planet or a large enough
286  allow us the create an approximate solution with small to zero boundary effects  portional of it cannot be modelled efficiently when considering small scale
287  on a model with a solvable size.  problems. It is impractical to calculate the solution for an infinite model and thus ABCs allow
288    us the create an approximate solution with small to zero boundary effects on a
289    model with a solvable size.
290
291  To dampen the waves, the method of \citet{Cerjan1985}  To dampen the waves, the method of \citet{Cerjan1985}
292  where the solution and the stress are multiplied by a damping function defined  where the solution and the stress are multiplied by a damping function defined
# Line 345  resulting boundary damping function can Line 348  resulting boundary damping function can
348  For stiff problems like the wave equation it is often prudent to implement  For stiff problems like the wave equation it is often prudent to implement
349  second order meshing. This creates a more accurate mesh approximation with some  second order meshing. This creates a more accurate mesh approximation with some
350  increased processing cost. To turn second order meshing on, the \verb!rectangle!  increased processing cost. To turn second order meshing on, the \verb!rectangle!
351  function accpets an \verb!order! keyword argument.  function accepts an \verb!order! keyword argument.
352  \begin{python}  \begin{python}
353  domain=Rectangle(l0=mx,l1=my,n0=ndx, n1=ndy,order=2) # create the domain  domain=Rectangle(l0=mx,l1=my,n0=ndx, n1=ndy,order=2) # create the domain
354  \end{python}  \end{python}
# Line 389  the mesh by half. Line 392  the mesh by half.
393  \sslist{example08c.py}  \sslist{example08c.py}
394  To make the problem more interesting we will now introduce an interface to the  To make the problem more interesting we will now introduce an interface to the
395  middle of the domain. Infact we will use the same domain as we did for heat flux  middle of the domain. In fact we will use the same domain as we did fora
396  in \autoref{CHAP HEAT 2}. The domain contains a syncline with two set of  different set of material properties on either side of the interface.
material properties on either side of the interface.
397
398  \begin{figure}[ht]  \begin{figure}[ht]
399  \begin{center}  \begin{center}
# Line 415  rho=Scalar(0,Function(domain)) Line 417  rho=Scalar(0,Function(domain))
417  rho.setTaggedValue("top",rho1)  rho.setTaggedValue("top",rho1)
418  rho.setTaggedValue("bottom",rho2)  rho.setTaggedValue("bottom",rho2)
419  \end{python}  \end{python}
420  Don't forget that teh source boudnary must also be tagged and added so it can be  Don't forget that the source boundary must also be tagged and added so it can
421  referenced  be referenced
422  \begin{python}  \begin{python}
423  # Add the subdomains and flux boundaries.  # Add the subdomains and flux boundaries.