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

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

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

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}
66    g=grad(u); stress=lam*trace(g)*kmat+mu*(g+transpose(g))
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    
102    \section{Multi-threading}
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

  ViewVC Help
Powered by ViewVC 1.1.26