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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 3989 - (show annotations)
Tue Sep 25 02:21:54 2012 UTC (7 years ago) by jfenwick
File MIME type: application/x-tex
File size: 16506 byte(s)
More copyright fixes.
pyvisi traces removed.
Some install doco
1
2 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
3 % Copyright (c) 2003-2012 by University of Queensland
4 % http://www.uq.edu.au
5 %
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 % Development until 2012 by Earth Systems Science Computational Center (ESSCC)
11 % Development since 2012 by School of Earth Sciences
12 %
13 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
14
15 \section{Seismic Wave Propagation in Two Dimensions}
16
17 \sslist{example08a.py}
18 We will now expand upon the previous chapter by introducing a vector form of
19 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
24 Wave propagation in the earth can be described by the elastic wave equation
25 \begin{equation} \label{eqn:wav} \index{wave equation}
26 \rho \frac{\partial^{2}u_{i}}{\partial t^2} - \frac{\partial
27 \sigma_{ij}}{\partial x_{j}} = 0
28 \end{equation}
29 where $\sigma$ is the stress given by
30 \begin{equation} \label{eqn:sigw}
31 \sigma _{ij} = \lambda u_{k,k} \delta_{ij} + \mu (
32 u_{i,j} + u_{j,i})
33 \end{equation}
34 and $\lambda$ and $\mu$ represent Lame's parameters. Specifically for seismic
35 waves, $\mu$ is the propagation materials shear modulus.
36 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 $\frac{\partial^{2}u_{i}}{\partial t^2}$ we can derive the
39 acceleration solution. Using $a$ we can see that \autoref{eqn:wav} becomes
40 \begin{equation} \label{eqn:wava}
41 \rho a_{i} - \frac{\partial
42 \sigma_{ij}}{\partial x_{j}} = 0
43 \end{equation}
44 Thus the problem will be solved for acceleration and then converted to
45 displacement using the backwards difference approximation as for the acoustic
46 example in the previous chapter.
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_{k,k} \delta_{ij}
53 \end{equation}
54 \begin{equation} \label{eqn:sigtrans}
55 \mu (u_{i,j} + u_{j,i})
56 \end{equation}
57 \end{subequations}
58 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 kronecker delta function with dimensions equivalent to $u$. The second term
61 \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 gradient of $u$ and the aforementioned operations. 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 acceleration
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
91 statement is 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 \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 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 The computation times on an increasing number of cores is outlined in
115 \autoref{tab:wpcores}.
116
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 \section{Vector source on the boundary}
138 \sslist{example08b.py}
139 For this particular example, we will introduce the source by applying a
140 displacement to the boundary during the initial time steps. The source will
141 again be
142 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
145 The first step is to choose an amplitude and create the source as in the
146 previous chapter.
147 \begin{python}
148 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 src_length = 40; print "src_length = ",src_length
155 # set initial values for first two time steps with source terms
156 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 \end{python}
162 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 \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 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 \begin{python}
171 src_dir=numpy.array([0.,-1.]) # defines direction of point source as down
172 y=y*src_dir
173 \end{python}
174 The function can then be applied as a boundary condition by setting it equal to
175 $y$ in the general form.
176 \begin{python}
177 mypde.setValue(y=y) #set the source as a function on the boundary
178 \end{python}
179 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 \begin{python}
183 # initial value of displacement at point source is constant (U0=0.01)
184 # for first two time steps
185 u=[0.0,0.0]*wherePositive(x)
186 u_m1=u
187 \end{python}
188
189 If the source is time progressive, $y$ can be updated during the
190 iteration stage. This is covered in the following section.
191
192 \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 \clearpage
214
215 \section{Time variant source}
216
217 \sslist{example08b.py}
218 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 this example the time variant source will be the derivative of a Gaussian curve
227 defined by the required dominant frequency (\autoref{fig:tvsource}).
228 \begin{python}
229 #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 source=np.zeros(ls,'float') # source array
237 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 \end{python}
247 \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
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 solution in the iterative section of the solver. This is done via;
270 \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 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
292 To dampen the waves, the method of \citet{Cerjan1985}
293 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 \gamma =\sqrt{\frac{| -log( \gamma _{b} ) |}{n^2}}
297 \end{equation}
298 \begin{equation}
299 y=e^{-(\gamma x)^2}
300 \end{equation}
301 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 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 resulting boundary damping function can be viewed in
346 \autoref{fig:abconds}.
347
348 \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 function accepts an \verb!order! keyword argument.
353 \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 \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
370 \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 \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 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
399 \begin{figure}[ht]
400 \begin{center}
401 \includegraphics[width=5in]{figures/gmsh-example08c.png}
402 \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 accept this domain. Multiple material parameters must now be deined and assigned
409 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 Don't forget that the source boundary must also be tagged and added so it can
422 be referenced
423 \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 It is now possible to solve the script as in the previous examples
429 (\autoref{fig:ex08cres}).
430
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