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

Revision 3054 - (hide annotations)
Wed Jun 30 02:22:25 2010 UTC (10 years, 6 months ago) by ahallam
File MIME type: application/x-tex
File size: 8412 byte(s)
Missing Figures...

 1 ahallam 2411 2 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 3 % 4 jfenwick 2881 % Copyright (c) 2003-2010 by University of Queensland 5 ahallam 2411 % 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 ahallam 2426 \section{Seismic Wave Propagation in Two Dimensions} 15 ahallam 3029 \editor{This chapter aims to expand the readers understanding of escript by 16 modelling the wave equations. 17 Challenges will include a second order differential (multiple initial 18 conditions). A new PDE to fit to the general form. Movement to a 3D problem 19 (maybe)??? } 20 ahallam 2426 21 ahallam 3029 \sslist{example08a.py} 22 ahallam 2426 23 ahallam 3029 We will now expand upon the previous chapter by introducing a vector form of 24 the wave equation. This means that the waves will have both a scalar magnitude, 25 but also a direction. This type of scenario is apparent in wave forms that 26 exhibit compressional and transverse particle motion. A common type of wave 27 that obeys this principle are seismic waves. 28 29 Wave propagation in the earth can be described by the elastic wave equation: 30 ahallam 2434 \begin{equation} \label{eqn:wav} \index{wave equation} 31 ahallam 3029 \rho \frac{\partial^{2}u\hackscore{i}}{\partial t^2} - \frac{\partial 32 \sigma\hackscore{ij}}{\partial x\hackscore{j}} = 0 33 ahallam 2434 \end{equation} 34 ahallam 3029 where $\sigma$ is the stress given by: 35 ahallam 2434 \begin{equation} \label{eqn:sigw} 36 ahallam 3029 \sigma \hackscore{ij} = \lambda u\hackscore{k,k} \delta\hackscore{ij} + \mu ( 37 u\hackscore{i,j} + u\hackscore{j,i}) 38 ahallam 2434 \end{equation} 39 ahallam 3029 where $\lambda$ and $\mu$ are the Lame Coefficients. Specifically $\mu$ is the 40 bulk modulus. The \refEq{eqn:wav} can be written with the Einstein summation 41 convention as: 42 ahallam 2434 \begin{equation} 43 ahallam 3029 \rho u\hackscore{i,tt} = \sigma\hackscore{ij,j} 44 ahallam 2434 \end{equation} 45 46 ahallam 3029 In a similar process to the previous chapter, we will use the acceleration 47 solution to solve this PDE. By substituting $a$ directly for 48 $\frac{\partial^{2}u\hackscore{i}}{\partial t^2}$ we can derive the 49 displacement solution. Using $a$ \refEq{eqn:wav} becomes; 50 \begin{equation} \label{eqn:wava} 51 \rho a\hackscore{i} - \frac{\partial 52 \sigma\hackscore{ij}}{\partial x\hackscore{j}} = 0 53 \end{equation} 54 ahallam 2434 55 ahallam 3029 \section{Vector source on the boundary} 56 For this particular example, we will introduce the source by applying a 57 displacment to the boundary during the initial time steps. The source will again 58 be 59 a radially propagating wave but due to the vector nature of the PDE used, a 60 direction will need to be applied to the source. 61 ahallam 2460 62 ahallam 3029 The first step is to choose an amplitude and create the source as in the 63 previous chapter. 64 ahallam 3025 \begin{python} 65 ahallam 3029 src_length = 20; print "src_length = ",src_length 66 # set initial values for first two time steps with source terms 67 y=U0*(cos(length(x-xc)*3.1415/src_length)+1)*whereNegative(length(x-xc)-src_leng 68 th) 69 ahallam 3025 \end{python} 70 ahallam 3029 where \verb xc is the source point on the boundary of the model. The source 71 direction is then defined as an $(x,y)$ array and multiplied by the source 72 function. The directional array must have a magnitude of $1$ otherwise the 73 amplitude of the source will become modified. For this example, the source is 74 directed in the $-y$ direction. 75 ahallam 3025 \begin{python} 76 ahallam 3029 src_dir=numpy.array([0.,-1.]) # defines direction of point source as down 77 y=y*src_dir 78 ahallam 3025 \end{python} 79 ahallam 3029 The function can then be applied as a boundary condition by setting it equal to 80 $y$ in the general form. 81 ahallam 3025 \begin{python} 82 ahallam 3029 mypde.setValue(y=y) #set the source as a function on the boundary 83 ahallam 3025 \end{python} 84 ahallam 3029 Because we are no longer using the source to define our initial condition to 85 the model, we must set the model state to zero for the first two time steps. 86 ahallam 3025 \begin{python} 87 ahallam 3029 # initial value of displacement at point source is constant (U0=0.01) 88 # for first two time steps 89 u=[0.0,0.0]*whereNegative(x) 90 ahallam 2469 u_m1=u 91 ahallam 3025 \end{python} 92 ahallam 2434 93 ahallam 3029 If the source will introduce energy to the system over a period longer than one 94 or two time steps (ie the initial conditions), $y$ can be updated during the 95 ahallam 3054 iteration stage. 96 97 \section{Time variant source} 98 99 \sslist{example08b.py} 100 101 Until this point, all of the wave propagation examples in this cookbook have 102 used impulsive sources which are smooth in space but not time. It is however, 103 advantageous to have a time smoothed source as it can reduce the temporal 104 frequency range and thus mitigate aliasing in the solution. 105 106 107 It is quite 108 simple to implement a source which is smooth in time. In addition to the 109 original source function the only extra requirement is a time function. For 110 this example it will be a decaying sinusoidal curve defined by; 111 \begin{python} 112 U0=0.1 # amplitude of point source 113 ls=100 # length of the source 114 source=np.zeros(ls,'float') # source array 115 g=np.log(0.01)/ls 116 for t in range(0,ls): 117 source[t]=np.exp(g*t)*U0*np.sin(2.*np.pi*t/(0.75*ls)) 118 \end{python} 119 120 We then build the source and the first two time steps via; 121 \begin{python} 122 # set initial values for first two time steps with source terms 123 y=source 124 *(cos(length(x-xc)*3.1415/src_length)+1)*whereNegative(length(x-xc)-src_length) 125 src_dir=numpy.array([0.,-1.]) # defines direction of point source as down 126 y=y*src_dir 127 mypde.setValue(y=y) #set the source as a function on the boundary 128 # initial value of displacement at point source is constant (U0=0.01) 129 # for first two time steps 130 u=[0.0,0.0]*whereNegative(x) 131 u_m1=u 132 \end{python} 133 134 Finally, for the length of the source, we are required to update each new 135 solution in the itterative section of the solver. This is done via; 136 \begin{python} 137 # increment loop values 138 t=t+h; n=n+1 139 if (n < ls): 140 y=source[n]**(cos(length(x-xc)*3.1415/src_length)+1)*\ 141 whereNegative(length(x-xc)-src_length) 142 y=y*src_dir; mypde.setValue(y=y) #set the source as a function on the 143 boundary 144 \end{python} 145 146 \section{Absorbing Boundary Conditions} 147 To mitigate the effect of the boundary on the model, absorbing boundary 148 conditions can be introduced. These conditions effectively dampen the wave 149 energy as they approach the bounday and thus prevent that energy from being 150 reflected. This type of approach is used typically when a model only represents 151 a small portion of the entire model, which in reality may have infinite bounds. 152 It is inpractical to calculate the solution for an infinite model and thus ABCs 153 allow us the create an approximate solution with small to zero boundary effects 154 on a model with a solvable size. 155 156 To dampen the waves, the method of Cerjan(1985) 157 \footnote{\textit{A nonreflecting boundary condition for discrete acoustic and 158 elastic wave equations}, 1985, Cerjan C, Geophysics 50, doi:10.1190/1.1441945} 159 where the solution and the stress are multiplied by a damping function defined 160 on $n$ nodes of the domain adjacent to the boundary, given by; 161 \begin{equation} 162 y= 163 \end{equation} 164 This is applied to the bounding 20-50 pts of the model using the location 165 specifiers of \esc; 166 \begin{python} 167 # Define where the boundary decay will be applied. 168 bn=30. 169 bleft=xstep*bn; bright=mx-(xstep*bn); bbot=my-(ystep*bn) 170 # btop=ystep*bn # don't apply to force boundary!!! 171 172 # locate these points in the domain 173 left=x-bleft; right=x-bright; bottom=x-bbot 174 175 tgamma=0.98 # decay value for exponential function 176 def calc_gamma(G,npts): 177 func=np.sqrt(abs(-1.*np.log(G)/(npts**2.))) 178 return func 179 180 gleft = calc_gamma(tgamma,bleft) 181 gright = calc_gamma(tgamma,bleft) 182 gbottom= calc_gamma(tgamma,ystep*bn) 183 184 print 'gamma', gleft,gright,gbottom 185 186 # calculate decay functions 187 def abc_bfunc(gamma,loc,x,G): 188 func=exp(-1.*(gamma*abs(loc-x))**2.) 189 return func 190 191 fleft=abc_bfunc(gleft,bleft,x,tgamma) 192 fright=abc_bfunc(gright,bright,x,tgamma) 193 fbottom=abc_bfunc(gbottom,bbot,x,tgamma) 194 # apply these functions only where relevant 195 abcleft=fleft*whereNegative(left) 196 abcright=fright*wherePositive(right) 197 abcbottom=fbottom*wherePositive(bottom) 198 # make sure the inside of the abc is value 1 199 abcleft=abcleft+whereZero(abcleft) 200 abcright=abcright+whereZero(abcright) 201 abcbottom=abcbottom+whereZero(abcbottom) 202 # multiply the conditions together to get a smooth result 203 abc=abcleft*abcright*abcbottom 204 \end{python} 205 Note that the boundary conditions are not applied to the surface, as this is 206 effectively a free surface where normal reflections would be experienced. The 207 resulting boundary damping function can be viewed in \ref{fig:abconds} 208 209 %\begin{figure}{ht} 210 %\centering 211 %\includegraphics[width=5in]{figures/abconds.png} 212 %\end{figure} 213 214