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

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

Parent Directory Parent Directory | Revision Log Revision Log


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

  ViewVC Help
Powered by ViewVC 1.1.26