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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 4286 - (hide annotations)
Thu Mar 7 04:28:11 2013 UTC (6 years, 6 months ago) by caltinay
File MIME type: application/x-tex
File size: 16507 byte(s)
Assorted spelling fixes.

1 ahallam 2411
2 jfenwick 3989 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
3 jfenwick 4154 % Copyright (c) 2003-2013 by University of Queensland
4 jfenwick 3989 % http://www.uq.edu.au
5 ahallam 2411 %
6     % Primary Business: Queensland, Australia
7     % Licensed under the Open Software License version 3.0
8     % http://www.opensource.org/licenses/osl-3.0.php
9     %
10 jfenwick 3989 % Development until 2012 by Earth Systems Science Computational Center (ESSCC)
11     % Development since 2012 by School of Earth Sciences
12     %
13     %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
14 ahallam 2411
15 ahallam 2426 \section{Seismic Wave Propagation in Two Dimensions}
16    
17 ahallam 3029 \sslist{example08a.py}
18     We will now expand upon the previous chapter by introducing a vector form of
19 ahallam 3153 the wave equation. This means that the waves will have not only a scalar
20     magnitude as for the pressure wave solution, but also a direction. This type of
21     scenario is apparent in wave types that exhibit compressional and transverse
22     particle motion. An example of this would be seismic waves.
23 ahallam 3029
24 ahallam 3063 Wave propagation in the earth can be described by the elastic wave equation
25 ahallam 2434 \begin{equation} \label{eqn:wav} \index{wave equation}
26 jfenwick 3308 \rho \frac{\partial^{2}u_{i}}{\partial t^2} - \frac{\partial
27     \sigma_{ij}}{\partial x_{j}} = 0
28 ahallam 2434 \end{equation}
29 ahallam 3063 where $\sigma$ is the stress given by
30 ahallam 2434 \begin{equation} \label{eqn:sigw}
31 jfenwick 3308 \sigma _{ij} = \lambda u_{k,k} \delta_{ij} + \mu (
32     u_{i,j} + u_{j,i})
33 ahallam 2434 \end{equation}
34 ahallam 3153 and $\lambda$ and $\mu$ represent Lame's parameters. Specifically for seismic
35     waves, $\mu$ is the propagation materials shear modulus.
36 ahallam 3029 In a similar process to the previous chapter, we will use the acceleration
37     solution to solve this PDE. By substituting $a$ directly for
38 jfenwick 3308 $\frac{\partial^{2}u_{i}}{\partial t^2}$ we can derive the
39 ahallam 3370 acceleration solution. Using $a$ we can see that \autoref{eqn:wav} becomes
40 ahallam 3029 \begin{equation} \label{eqn:wava}
41 jfenwick 3308 \rho a_{i} - \frac{\partial
42     \sigma_{ij}}{\partial x_{j}} = 0
43 ahallam 3029 \end{equation}
44 ahallam 3153 Thus the problem will be solved for acceleration and then converted to
45 ahallam 3370 displacement using the backwards difference approximation as for the acoustic
46     example in the previous chapter.
47 ahallam 2434
48 ahallam 3063 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 jfenwick 3308 \lambda u_{k,k} \delta_{ij}
53 ahallam 3063 \end{equation}
54     \begin{equation} \label{eqn:sigtrans}
55 jfenwick 3308 \mu (u_{i,j} + u_{j,i})
56 ahallam 3063 \end{equation}
57     \end{subequations}
58 jfenwick 3308 One simply recognizes in \autoref{eqn:sigtrace} that $u_{k,k}$ is the
59     trace of the displacement solution and that $\delta_{ij}$ is the
60 ahallam 3063 kronecker delta function with dimensions equivalent to $u$. The second term
61 ahallam 3232 \autoref{eqn:sigtrans} is the sum of $u$ with its own transpose. Putting these
62     facts together we see that the spatial differential of the stress is given by the
63 ahallam 3370 gradient of $u$ and the aforementioned operations. This value is then submitted
64 ahallam 3063 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 ahallam 3370 accel = mypde.getSolution() #get PDE solution for acceleration
74 ahallam 3063 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 ahallam 3153 the modelling and the final time step. With these parameters set, an if
91     statement is introduced to the time loop
92 ahallam 3063 \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 ahallam 3370 The wave equation solution can be quite demanding on CPU time. Enhancements can
104 ahallam 3063 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 ahallam 3370 \esc is called from the shell. To use multiple threads \esc is called using
107     the \verb!-t! option with an integer argument for the number of threads
108 ahallam 3063 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 ahallam 3370 The computation times on an increasing number of cores is outlined in
115 ahallam 3232 \autoref{tab:wpcores}.
116 ahallam 3153
117     \begin{table}[ht]
118     \begin{center}
119     \caption{Computation times for an increasing number of cores.}
120     \label{tab:wpcores}
121     \begin{tabular}{| c | c |}
122     \hline
123     Number of Cores & Time (s) \\
124     \hline
125     1 & 691.0 \\
126     2 & 400.0 \\
127     3 & 305.0 \\
128     4 & 328.0 \\
129     5 & 323.0 \\
130     6 & 292.0 \\
131     7 & 282.0 \\
132     8 & 445.0 \\ \hline
133     \end{tabular}
134     \end{center}
135     \end{table}
136    
137 ahallam 3029 \section{Vector source on the boundary}
138 ahallam 3153 \sslist{example08b.py}
139 ahallam 3029 For this particular example, we will introduce the source by applying a
140 ahallam 3370 displacement to the boundary during the initial time steps. The source will
141     again be
142 ahallam 3029 a radially propagating wave but due to the vector nature of the PDE used, a
143     direction will need to be applied to the source.
144 ahallam 2460
145 ahallam 3029 The first step is to choose an amplitude and create the source as in the
146     previous chapter.
147 ahallam 3025 \begin{python}
148 ahallam 3063 U0=0.01 # amplitude of point source
149     # will introduce a spherical source at middle left of bottom face
150     xc=[ndx/2,0]
151    
152     ############################################FIRST TIME STEPS AND SOURCE
153     # define small radius around point xc
154 ahallam 3153 src_length = 40; print "src_length = ",src_length
155 ahallam 3029 # set initial values for first two time steps with source terms
156 ahallam 3153 xb=FunctionOnBoundary(domain).getX()
157     y=source[0]*(cos(length(x-xc)*3.1415/src_length)+1)*\
158     whereNegative(length(xb-src_length))
159     src_dir=numpy.array([0.,1.]) # defines direction of point source as down
160     y=y*src_dir
161 ahallam 3025 \end{python}
162 ahallam 3153 where \verb xc is the source point on the boundary of the model. Note that
163     because the source is specifically located on the boundary, we have used the
164 ahallam 3370 \verb!FunctionOnBoundary! call to ensure the nodes used to define the source
165     are also located upon the boundary. These boundary nodes are passed to
166     source as \verb!xb!. The source direction is then defined as an $(x,y)$ array and multiplied by the
167 ahallam 3153 source function. The directional array must have a magnitude of $\left| 1
168     \right| $ otherwise the amplitude of the source will become modified. For this
169     example, the source is directed in the $-y$ direction.
170 ahallam 3025 \begin{python}
171 ahallam 3029 src_dir=numpy.array([0.,-1.]) # defines direction of point source as down
172     y=y*src_dir
173 ahallam 3025 \end{python}
174 ahallam 3029 The function can then be applied as a boundary condition by setting it equal to
175     $y$ in the general form.
176 ahallam 3025 \begin{python}
177 ahallam 3029 mypde.setValue(y=y) #set the source as a function on the boundary
178 ahallam 3025 \end{python}
179 ahallam 3063 The final step is to qualify the initial conditions. Due to the fact that we are
180     no longer using the source to define our initial condition to the model, we
181     must set the model state to zero for the first two time steps.
182 ahallam 3025 \begin{python}
183 ahallam 3029 # initial value of displacement at point source is constant (U0=0.01)
184     # for first two time steps
185 ahallam 3153 u=[0.0,0.0]*wherePositive(x)
186 ahallam 2469 u_m1=u
187 ahallam 3025 \end{python}
188 ahallam 2434
189 ahallam 3063 If the source is time progressive, $y$ can be updated during the
190     iteration stage. This is covered in the following section.
191 ahallam 3054
192 ahallam 3063 \begin{figure}[htp]
193     \centering
194     \subfigure[Example 08a at 0.025s ]{
195     \includegraphics[width=3in]{figures/ex08pw50.png}
196     \label{fig:ex08pw50}
197     }
198     \subfigure[Example 08a at 0.175s ]{
199     \includegraphics[width=3in]{figures/ex08pw350.png}
200     \label{fig:ex08pw350}
201     } \\
202     \subfigure[Example 08a at 0.325s ]{
203     \includegraphics[width=3in]{figures/ex08pw650.png}
204     \label{fig:ex08pw650}
205     }
206     \subfigure[Example 08a at 0.475s ]{
207     \includegraphics[width=3in]{figures/ex08pw950.png}
208     \label{fig:ex08pw950}
209     }
210     \label{fig:ex08pw}
211     \caption{Results of Example 08 at various times.}
212     \end{figure}
213 ahallam 3153 \clearpage
214 ahallam 3063
215 ahallam 3054 \section{Time variant source}
216    
217 ahallam 3153 \sslist{example08b.py}
218 ahallam 3054 Until this point, all of the wave propagation examples in this cookbook have
219     used impulsive sources which are smooth in space but not time. It is however,
220     advantageous to have a time smoothed source as it can reduce the temporal
221     frequency range and thus mitigate aliasing in the solution.
222    
223     It is quite
224     simple to implement a source which is smooth in time. In addition to the
225     original source function the only extra requirement is a time function. For
226 ahallam 3370 this example the time variant source will be the derivative of a Gaussian curve
227 ahallam 3232 defined by the required dominant frequency (\autoref{fig:tvsource}).
228 ahallam 3054 \begin{python}
229 ahallam 3153 #Creating the time function of the source.
230     dfeq=50 #Dominant Frequency
231     a = 2.0 * (np.pi * dfeq)**2.0
232     t0 = 5.0 / (2.0 * np.pi * dfeq)
233     srclength = 5. * t0
234     ls = int(srclength/h)
235     print 'source length',ls
236 ahallam 3054 source=np.zeros(ls,'float') # source array
237 ahallam 3153 ampmax=0
238     for it in range(0,ls):
239     t = it*h
240     tt = t-t0
241     dum1 = np.exp(-a * tt * tt)
242     source[it] = -2. * a * tt * dum1
243     if (abs(source[it]) > ampmax):
244     ampmax = abs(source[it])
245     time[t]=t*h
246 ahallam 3054 \end{python}
247 ahallam 3153 \begin{figure}[ht]
248     \centering
249     \includegraphics[width=3in]{figures/source.png}
250     \caption{Time variant source with a dominant frequency of 50Hz.}
251     \label{fig:tvsource}
252     \end{figure}
253 ahallam 3054
254     We then build the source and the first two time steps via;
255     \begin{python}
256     # set initial values for first two time steps with source terms
257     y=source[0]
258     *(cos(length(x-xc)*3.1415/src_length)+1)*whereNegative(length(x-xc)-src_length)
259     src_dir=numpy.array([0.,-1.]) # defines direction of point source as down
260     y=y*src_dir
261     mypde.setValue(y=y) #set the source as a function on the boundary
262     # initial value of displacement at point source is constant (U0=0.01)
263     # for first two time steps
264     u=[0.0,0.0]*whereNegative(x)
265     u_m1=u
266     \end{python}
267    
268     Finally, for the length of the source, we are required to update each new
269 ahallam 3370 solution in the iterative section of the solver. This is done via;
270 ahallam 3054 \begin{python}
271     # increment loop values
272     t=t+h; n=n+1
273     if (n < ls):
274     y=source[n]**(cos(length(x-xc)*3.1415/src_length)+1)*\
275     whereNegative(length(x-xc)-src_length)
276     y=y*src_dir; mypde.setValue(y=y) #set the source as a function on the
277     boundary
278     \end{python}
279    
280     \section{Absorbing Boundary Conditions}
281     To mitigate the effect of the boundary on the model, absorbing boundary
282 ahallam 3370 conditions can be introduced. These conditions effectively dampen the wave energy
283     as they approach the boundary and thus prevent that energy from being reflected.
284     This type of approach is typically used when a model is shrunk to decrease
285     computational requirements. In practise this applies to almost all models,
286     especially in earth sciences where the entire planet or a large enough
287     portional of it cannot be modelled efficiently when considering small scale
288     problems. It is impractical to calculate the solution for an infinite model and thus ABCs allow
289     us the create an approximate solution with small to zero boundary effects on a
290     model with a solvable size.
291 ahallam 3054
292 ahallam 3232 To dampen the waves, the method of \citet{Cerjan1985}
293 ahallam 3054 where the solution and the stress are multiplied by a damping function defined
294     on $n$ nodes of the domain adjacent to the boundary, given by;
295     \begin{equation}
296 jfenwick 3308 \gamma =\sqrt{\frac{| -log( \gamma _{b} ) |}{n^2}}
297 ahallam 3054 \end{equation}
298 ahallam 3153 \begin{equation}
299     y=e^{-(\gamma x)^2}
300     \end{equation}
301 ahallam 3054 This is applied to the bounding 20-50 pts of the model using the location
302     specifiers of \esc;
303     \begin{python}
304     # Define where the boundary decay will be applied.
305     bn=30.
306     bleft=xstep*bn; bright=mx-(xstep*bn); bbot=my-(ystep*bn)
307     # btop=ystep*bn # don't apply to force boundary!!!
308    
309     # locate these points in the domain
310     left=x[0]-bleft; right=x[0]-bright; bottom=x[1]-bbot
311    
312     tgamma=0.98 # decay value for exponential function
313     def calc_gamma(G,npts):
314     func=np.sqrt(abs(-1.*np.log(G)/(npts**2.)))
315     return func
316    
317     gleft = calc_gamma(tgamma,bleft)
318     gright = calc_gamma(tgamma,bleft)
319     gbottom= calc_gamma(tgamma,ystep*bn)
320    
321     print 'gamma', gleft,gright,gbottom
322    
323     # calculate decay functions
324     def abc_bfunc(gamma,loc,x,G):
325     func=exp(-1.*(gamma*abs(loc-x))**2.)
326     return func
327    
328     fleft=abc_bfunc(gleft,bleft,x[0],tgamma)
329     fright=abc_bfunc(gright,bright,x[0],tgamma)
330     fbottom=abc_bfunc(gbottom,bbot,x[1],tgamma)
331     # apply these functions only where relevant
332     abcleft=fleft*whereNegative(left)
333     abcright=fright*wherePositive(right)
334     abcbottom=fbottom*wherePositive(bottom)
335     # make sure the inside of the abc is value 1
336     abcleft=abcleft+whereZero(abcleft)
337     abcright=abcright+whereZero(abcright)
338     abcbottom=abcbottom+whereZero(abcbottom)
339     # multiply the conditions together to get a smooth result
340     abc=abcleft*abcright*abcbottom
341     \end{python}
342     Note that the boundary conditions are not applied to the surface, as this is
343 ahallam 3153 effectively a free surface where normal reflections would be experienced.
344     Special conditions can be introduced at this surface if they are known. The
345 ahallam 3232 resulting boundary damping function can be viewed in
346     \autoref{fig:abconds}.
347 ahallam 3054
348 ahallam 3153 \section{Second order Meshing}
349     For stiff problems like the wave equation it is often prudent to implement
350     second order meshing. This creates a more accurate mesh approximation with some
351     increased processing cost. To turn second order meshing on, the \verb!rectangle!
352 ahallam 3370 function accepts an \verb!order! keyword argument.
353 ahallam 3153 \begin{python}
354     domain=Rectangle(l0=mx,l1=my,n0=ndx, n1=ndy,order=2) # create the domain
355     \end{python}
356     Other pycad functions and objects have similar keyword arguments for higher
357     order meshing.
358    
359     Note that when implementing second order meshing, a smaller timestep is required
360     then for first order meshes as the second order essentially reduces the size of
361     the mesh by half.
362    
363 ahallam 3063 \begin{figure}[ht]
364     \centering
365     \includegraphics[width=5in]{figures/ex08babc.png}
366     \label{fig:abconds}
367     \caption{Absorbing boundary conditions for example08b.py}
368     \end{figure}
369 ahallam 3054
370 ahallam 3063 \begin{figure}[htp]
371     \centering
372     \subfigure[Example 08b at 0.03s ]{
373     \includegraphics[width=3in]{figures/ex08sw060.png}
374     \label{fig:ex08pw060}
375     }
376     \subfigure[Example 08b at 0.16s ]{
377     \includegraphics[width=3in]{figures/ex08sw320.png}
378     \label{fig:ex08pw320}
379     } \\
380     \subfigure[Example 08b at 0.33s ]{
381     \includegraphics[width=3in]{figures/ex08sw660.png}
382     \label{fig:ex08pw660}
383     }
384     \subfigure[Example 08b at 0.44s ]{
385     \includegraphics[width=3in]{figures/ex08sw880.png}
386     \label{fig:ex08pw880}
387     }
388     \label{fig:ex08pw}
389     \caption{Results of Example 08b at various times.}
390     \end{figure}
391 ahallam 3153 \clearpage
392    
393     \section{Pycad example}
394     \sslist{example08c.py}
395     To make the problem more interesting we will now introduce an interface to the
396 ahallam 3370 middle of the domain. In fact we will use the same domain as we did fora
397     different set of material properties on either side of the interface.
398 ahallam 3153
399     \begin{figure}[ht]
400     \begin{center}
401 ahallam 3168 \includegraphics[width=5in]{figures/gmsh-example08c.png}
402 ahallam 3153 \caption{Domain geometry for example08c.py showing line tangents.}
403     \label{fig:ex08cgeo}
404     \end{center}
405     \end{figure}
406    
407     It is simple enough to slightly modify the scripts of the previous sections to
408 caltinay 4286 accept this domain. Multiple material parameters must now be defined and assigned
409 ahallam 3153 to specific tagged areas. Again this is done via
410     \begin{python}
411     lam=Scalar(0,Function(domain))
412     lam.setTaggedValue("top",lam1)
413     lam.setTaggedValue("bottom",lam2)
414     mu=Scalar(0,Function(domain))
415     mu.setTaggedValue("top",mu1)
416     mu.setTaggedValue("bottom",mu2)
417     rho=Scalar(0,Function(domain))
418     rho.setTaggedValue("top",rho1)
419     rho.setTaggedValue("bottom",rho2)
420     \end{python}
421 ahallam 3370 Don't forget that the source boundary must also be tagged and added so it can
422     be referenced
423 ahallam 3153 \begin{python}
424     # Add the subdomains and flux boundaries.
425     d.addItems(PropertySet("top",tblock),PropertySet("bottom",bblock),\
426     PropertySet("linetop",l30))
427     \end{python}
428 ahallam 3370 It is now possible to solve the script as in the previous examples
429     (\autoref{fig:ex08cres}).
430 ahallam 3153
431     \begin{figure}[ht]
432     \centering
433     \includegraphics[width=4in]{figures/ex08c2601.png}
434     \caption{Modelling results of example08c.py at 0.2601 seconds. Notice the
435     refraction of the wave front about the boundary between the two layers.}
436     \label{fig:ex08cres}
437     \end{figure}
438    
439    

  ViewVC Help
Powered by ViewVC 1.1.26