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

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

Parent Directory Parent Directory | Revision Log Revision Log


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

  ViewVC Help
Powered by ViewVC 1.1.26