 # Contents of /branches/doubleplusgood/doc/cookbook/example08.tex

Revision 4345 - (show annotations)
Fri Mar 29 07:09:41 2013 UTC (6 years, 9 months ago) by jfenwick
File MIME type: application/x-tex
File size: 16507 byte(s)
Spelling fixes

 1 2 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 3 % Copyright (c) 2003-2013 by University of Queensland 4 5 % 6 % Primary Business: Queensland, Australia 7 % Licensed under the Open Software License version 3.0 8 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*(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 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-bleft; right=x-bright; bottom=x-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,tgamma) 329 fright=abc_bfunc(gright,bright,x,tgamma) 330 fbottom=abc_bfunc(gbottom,bbot,x,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 defined 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