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

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 5 ahallam 2411 % 6 % Primary Business: Queensland, Australia 7 % Licensed under the Open Software License version 3.0 8 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