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

revision 3062 by ahallam, Wed Jun 30 02:22:25 2010 UTC revision 3063 by ahallam, Thu Jul 15 02:57:46 2010 UTC
# Line 12  Line 12
12  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
13
14  \section{Seismic Wave Propagation in Two Dimensions}  \section{Seismic Wave Propagation in Two Dimensions}
\editor{This chapter aims to expand the readers understanding of escript by
modelling the wave equations.
Challenges will include a second order differential (multiple initial
conditions). A new PDE to fit to the general form. Movement to a 3D problem
(maybe)??? }
15
16  \sslist{example08a.py}  \sslist{example08a.py}
17
# Line 26  but also a direction. This type of scena Line 21  but also a direction. This type of scena
21  exhibit compressional and transverse particle motion. A common type of wave  exhibit compressional and transverse particle motion. A common type of wave
22  that obeys this principle are seismic waves.  that obeys this principle are seismic waves.
23
24  Wave propagation in the earth can be described by the elastic wave equation:  Wave propagation in the earth can be described by the elastic wave equation
25  \begin{equation} \label{eqn:wav} \index{wave equation}  \begin{equation} \label{eqn:wav} \index{wave equation}
26  \rho \frac{\partial^{2}u\hackscore{i}}{\partial t^2} - \frac{\partial  \rho \frac{\partial^{2}u\hackscore{i}}{\partial t^2} - \frac{\partial
27  \sigma\hackscore{ij}}{\partial x\hackscore{j}} = 0  \sigma\hackscore{ij}}{\partial x\hackscore{j}} = 0
28  \end{equation}  \end{equation}
29  where $\sigma$ is the stress given by:  where $\sigma$ is the stress given by
30  \begin{equation} \label{eqn:sigw}  \begin{equation} \label{eqn:sigw}
31   \sigma \hackscore{ij} = \lambda u\hackscore{k,k} \delta\hackscore{ij} + \mu (   \sigma \hackscore{ij} = \lambda u\hackscore{k,k} \delta\hackscore{ij} + \mu (
32  u\hackscore{i,j} + u\hackscore{j,i})  u\hackscore{i,j} + u\hackscore{j,i})
33  \end{equation}  \end{equation}
34  where $\lambda$ and $\mu$ are the Lame Coefficients. Specifically $\mu$ is the  with $\lambda$ and $\mu$ representing the Lame Coefficients. Specifically
35  bulk modulus. The \refEq{eqn:wav} can be written with the Einstein summation  $\mu$ is the bulk modulus.
convention as:
\begin{equation}
\rho u\hackscore{i,tt} = \sigma\hackscore{ij,j}
\end{equation}

36  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
37  solution to solve this PDE. By substituting $a$ directly for  solution to solve this PDE. By substituting $a$ directly for
38  $\frac{\partial^{2}u\hackscore{i}}{\partial t^2}$ we can derive the  $\frac{\partial^{2}u\hackscore{i}}{\partial t^2}$ we can derive the
39  displacement solution. Using $a$ \refEq{eqn:wav} becomes;  displacement solution. Using $a$ we can see that \refEq{eqn:wav} becomes
40  \begin{equation} \label{eqn:wava}  \begin{equation} \label{eqn:wava}
41  \rho a\hackscore{i} - \frac{\partial  \rho a\hackscore{i} - \frac{\partial
42  \sigma\hackscore{ij}}{\partial x\hackscore{j}} = 0  \sigma\hackscore{ij}}{\partial x\hackscore{j}} = 0
43  \end{equation}  \end{equation}
44    This is very similar to our acceleration solution in the previous chapter. Thus
45    the problem will be solved for acceleration and for displacement using the
46    backwards difference approximation.
47
48    Consider now the stress $\sigma$. One can see that the stress consists of two
49    distinct terms:
50    \begin{subequations}
51    \begin{equation} \label{eqn:sigtrace}
52    \lambda u\hackscore{k,k} \delta\hackscore{ij}
53    \end{equation}
54    \begin{equation} \label{eqn:sigtrans}
55    \mu (u\hackscore{i,j} + u\hackscore{j,i})
56    \end{equation}
57    \end{subequations}
58    One simply recognizes in \ref{eqn:sigtrace} that $u\hackscore{k,k}$ is the
59    trace of the displacement solution and that $\delta\hackscore{ij}$ is the
60    kronecker delta function with dimensions equivalent to $u$. The second term
61    \ref{eqn:sigtrans} is the sum of $u$ with its own transpose. Putting these facts
62    together we see that the spatial differential of the stress is given by the
63    gradient of $u$ and the aforementioned opperations. This value is then submitted
64    to the \esc PDE as $X$.
65    \begin{python}
67    mypde.setValue(X=-stress) # set PDE values
68    \end{python}
69    The solution is then obtained via the usual method and the displacement is
70    calculated so that the memory variables can be updated for the next time
71    iteration.
72    \begin{python}
73    accel = mypde.getSolution() #get PDE solution for accelleration
74    u_p1=(2.*u-u_m1)+h*h*accel #calculate displacement
75    u_m1=u; u=u_p1 # shift values by 1
76    \end{python}
77
78    Saving the data has been handled slightly differently in this example. The VTK
79    files generated can be quite large and take a significant amount of time to save
80    to the hard disk. To avoid doing this at every iteration a test is devised which
81    saves only at specific time intervals.
82
83    To do this there are two new parameters in our script.
84    \begin{python}
85    # data recording times
86    rtime=0.0 # first time to record
87    rtime_inc=tend/20.0 # time increment to record
88    \end{python}
89    Currently the PDE solution will be saved to file $20$ times between the start of
90    the modelling and the final time step. With these parameters set an if catch is
91    introduced to the time loop
92    \begin{python}
93    if (t >= rtime):
94        saveVTK(os.path.join(savepath,"ex08a.%05d.vtu"%n),displacement=length(u),\
95            acceleration=length(accel),tensor=stress)
96            rtime=rtime+rtime_inc #increment data save time
97    \end{python}
98    \verb!t! is the time counter. Whenever the recording time \verb!rtime! is less
99    then \verb!t! the solution is saved and \verb!rtime! is incremented. This
100    limits the number of outputs and increases the speed of the solver.
101
103    The wave equation solution can be quite demanding on cpu time. Enhancements can
104    be made by accessing multiple threads or cores on your computer. This does not
105    require any modification to the solution script and only comes into play when
106    escript is called from the shell. To use multiple threads \esc is called using
107    the \verb!-t! option with an interger argument for the number of threads
108    required. For example
109    \begin{verbatim}
110    $escript -t 4 example08a.py 111 \end{verbatim} 112 would call the script in this section and solve it using 4 threads. 113 114 \section{Vector source on the boundary} \section{Vector source on the boundary} 115 For this particular example, we will introduce the source by applying a For this particular example, we will introduce the source by applying a # Line 62 direction will need to be applied to the Line 121 direction will need to be applied to the 121 The first step is to choose an amplitude and create the source as in the The first step is to choose an amplitude and create the source as in the 122 previous chapter. previous chapter. 123 \begin{python} \begin{python} 124 src_length = 20; print "src_length = ",src_length U0=0.01 # amplitude of point source 125 # will introduce a spherical source at middle left of bottom face 126 xc=[ndx/2,0] 127 128 ############################################FIRST TIME STEPS AND SOURCE 129 # define small radius around point xc 130 src_length = 20; print "src_length = ",src_length 131 # set initial values for first two time steps with source terms # set initial values for first two time steps with source terms 132 y=U0*(cos(length(x-xc)*3.1415/src_length)+1)*whereNegative(length(x-xc)-src_leng y=U0*(cos(length(x-xc)*3.1415/src_length)+1)*whereNegative(length(x-xc)\ 133 th) -src_length) 134 \end{python} \end{python} 135 where \verb xc is the source point on the boundary of the model. The source where \verb xc is the source point on the boundary of the model. The source 136 direction is then defined as an$(x,y)$array and multiplied by the source direction is then defined as an$(x,y)$array and multiplied by the source 137 function. The directional array must have a magnitude of$1$otherwise the function. The directional array must have a magnitude of$\left| 1 \right| $138 otherwise the 139 amplitude of the source will become modified. For this example, the source is amplitude of the source will become modified. For this example, the source is 140 directed in the$-y$direction. directed in the$-y$direction. 141 \begin{python} \begin{python} # Line 81$y$in the general form. Line 147$y$in the general form. 147 \begin{python} \begin{python} 148 mypde.setValue(y=y) #set the source as a function on the boundary mypde.setValue(y=y) #set the source as a function on the boundary 149 \end{python} \end{python} 150 Because we are no longer using the source to define our initial condition to The final step is to qualify the initial conditions. Due to the fact that we are 151 the model, we must set the model state to zero for the first two time steps. no longer using the source to define our initial condition to the model, we 152 must set the model state to zero for the first two time steps. 153 \begin{python} \begin{python} 154 # initial value of displacement at point source is constant (U0=0.01) # initial value of displacement at point source is constant (U0=0.01) 155 # for first two time steps # for first two time steps # Line 90 u=[0.0,0.0]*whereNegative(x) Line 157 u=[0.0,0.0]*whereNegative(x) 157 u_m1=u u_m1=u 158 \end{python} \end{python} 159 160 If the source will introduce energy to the system over a period longer than one If the source is time progressive,$y$can be updated during the 161 or two time steps (ie the initial conditions),$y\$ can be updated during the  iteration stage. This is covered in the following section.
162  iteration stage.
163    \begin{figure}[htp]
164    \centering
165    \subfigure[Example 08a at 0.025s ]{
166    \includegraphics[width=3in]{figures/ex08pw50.png}
167    \label{fig:ex08pw50}
168    }
169    \subfigure[Example 08a at 0.175s ]{
170    \includegraphics[width=3in]{figures/ex08pw350.png}
171    \label{fig:ex08pw350}
172    } \\
173    \subfigure[Example 08a at 0.325s ]{
174    \includegraphics[width=3in]{figures/ex08pw650.png}
175    \label{fig:ex08pw650}
176    }
177    \subfigure[Example 08a at 0.475s ]{
178    \includegraphics[width=3in]{figures/ex08pw950.png}
179    \label{fig:ex08pw950}
180    }
181    \label{fig:ex08pw}
182    \caption{Results of Example 08 at various times.}
183    \end{figure}
184
185  \section{Time variant source}  \section{Time variant source}
186
187  \sslist{example08b.py}  \sslist{example08a.py}
188
189  Until this point, all of the wave propagation examples in this cookbook have  Until this point, all of the wave propagation examples in this cookbook have
190  used impulsive sources which are smooth in space but not time. It is however,  used impulsive sources which are smooth in space but not time. It is however,
# Line 204  abc=abcleft*abcright*abcbottom Line 292  abc=abcleft*abcright*abcbottom
292  \end{python}  \end{python}
293  Note that the boundary conditions are not applied to the surface, as this is  Note that the boundary conditions are not applied to the surface, as this is
294  effectively a free surface where normal reflections would be experienced. The  effectively a free surface where normal reflections would be experienced. The
295  resulting boundary damping function can be viewed in \ref{fig:abconds}  resulting boundary damping function can be viewed in Figure \ref{fig:abconds}

%\begin{figure}{ht}
%\centering
%\includegraphics[width=5in]{figures/abconds.png}
%\end{figure}

296
297    \begin{figure}[ht]
298     \centering
299     \includegraphics[width=5in]{figures/ex08babc.png}
300     \label{fig:abconds}
301     \caption{Absorbing boundary conditions for example08b.py}
302    \end{figure}
303
304    \begin{figure}[htp]
305    \centering
306    \subfigure[Example 08b at 0.03s ]{
307    \includegraphics[width=3in]{figures/ex08sw060.png}
308    \label{fig:ex08pw060}
309    }
310    \subfigure[Example 08b at 0.16s ]{
311    \includegraphics[width=3in]{figures/ex08sw320.png}
312    \label{fig:ex08pw320}
313    } \\
314    \subfigure[Example 08b at 0.33s ]{
315    \includegraphics[width=3in]{figures/ex08sw660.png}
316    \label{fig:ex08pw660}
317    }
318    \subfigure[Example 08b at 0.44s ]{
319    \includegraphics[width=3in]{figures/ex08sw880.png}
320    \label{fig:ex08pw880}
321    }
322    \label{fig:ex08pw}
323    \caption{Results of Example 08b at various times.}
324    \end{figure}

Legend:
 Removed from v.3062 changed lines Added in v.3063