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

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 7 % 8 % Primary Business: Queensland, Australia 9 % Licensed under the Open Software License version 3.0 10 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*(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 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-bleft; right=x-bright; bottom=x-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,tgamma) 328 fright=abc_bfunc(gright,bright,x,tgamma) 329 fbottom=abc_bfunc(gbottom,bbot,x,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