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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 3232 - (show annotations)
Fri Oct 1 02:08:38 2010 UTC (8 years, 11 months ago) by ahallam
File MIME type: application/x-tex
File size: 16368 byte(s)
Updates to cookbook, new chapter DC Resis and IP. Using new packages hyperref, natbib
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\hackscore{i}}{\partial t^2} - \frac{\partial
26 \sigma\hackscore{ij}}{\partial x\hackscore{j}} = 0
27 \end{equation}
28 where $\sigma$ is the stress given by
29 \begin{equation} \label{eqn:sigw}
30 \sigma \hackscore{ij} = \lambda u\hackscore{k,k} \delta\hackscore{ij} + \mu (
31 u\hackscore{i,j} + u\hackscore{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\hackscore{i}}{\partial t^2}$ we can derive the
38 displacement solution. Using $a$ we can see that \autoref{eqn:wav} becomes
39 \begin{equation} \label{eqn:wava}
40 \rho a\hackscore{i} - \frac{\partial
41 \sigma\hackscore{ij}}{\partial x\hackscore{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.
45
46 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 \lambda u\hackscore{k,k} \delta\hackscore{ij}
51 \end{equation}
52 \begin{equation} \label{eqn:sigtrans}
53 \mu (u\hackscore{i,j} + u\hackscore{j,i})
54 \end{equation}
55 \end{subequations}
56 One simply recognizes in \autoref{eqn:sigtrace} that $u\hackscore{k,k}$ is the
57 trace of the displacement solution and that $\delta\hackscore{ij}$ is the
58 kronecker delta function with dimensions equivalent to $u$. The second term
59 \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 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 the modelling and the final time step. With these parameters set, an if
89 statement is introduced to the time loop
90 \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 The computation times on an increasing number of cores is outlines in
113 \autoref{tab:wpcores}.
114
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 \section{Vector source on the boundary}
136 \sslist{example08b.py}
137 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
143 The first step is to choose an amplitude and create the source as in the
144 previous chapter.
145 \begin{python}
146 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 src_length = 40; print "src_length = ",src_length
153 # set initial values for first two time steps with source terms
154 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 \end{python}
160 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 \begin{python}
169 src_dir=numpy.array([0.,-1.]) # defines direction of point source as down
170 y=y*src_dir
171 \end{python}
172 The function can then be applied as a boundary condition by setting it equal to
173 $y$ in the general form.
174 \begin{python}
175 mypde.setValue(y=y) #set the source as a function on the boundary
176 \end{python}
177 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 \begin{python}
181 # initial value of displacement at point source is constant (U0=0.01)
182 # for first two time steps
183 u=[0.0,0.0]*wherePositive(x)
184 u_m1=u
185 \end{python}
186
187 If the source is time progressive, $y$ can be updated during the
188 iteration stage. This is covered in the following section.
189
190 \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 \clearpage
212
213 \section{Time variant source}
214
215 \sslist{example08b.py}
216 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 this example the time variant source will be the derivative of a gausian curve
225 defined by the required dominant frequency (\autoref{fig:tvsource}).
226 \begin{python}
227 #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 source=np.zeros(ls,'float') # source array
235 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 \end{python}
245 \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
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 To dampen the waves, the method of \citet{Cerjan1985}
289 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 \gamma =\sqrt{\frac{| -log( \gamma \hackscore{b} ) |}{n^2}}
293 \end{equation}
294 \begin{equation}
295 y=e^{-(\gamma x)^2}
296 \end{equation}
297 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 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 resulting boundary damping function can be viewed in
342 \autoref{fig:abconds}.
343
344 \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 \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
366 \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 \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 in \autoref{CHAP HEAT 2}. The domain contains a syncline with two set of
394 material properties on either side of the interface.
395
396 \begin{figure}[ht]
397 \begin{center}
398 \includegraphics[width=5in]{figures/gmsh-example08c.png}
399 \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