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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 3054 - (show annotations)
Wed Jun 30 02:22:25 2010 UTC (10 years, 5 months ago) by ahallam
File MIME type: application/x-tex
File size: 8412 byte(s)
Missing Figures...
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 \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
21 \sslist{example08a.py}
22
23 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 \begin{equation} \label{eqn:wav} \index{wave equation}
31 \rho \frac{\partial^{2}u\hackscore{i}}{\partial t^2} - \frac{\partial
32 \sigma\hackscore{ij}}{\partial x\hackscore{j}} = 0
33 \end{equation}
34 where $\sigma$ is the stress given by:
35 \begin{equation} \label{eqn:sigw}
36 \sigma \hackscore{ij} = \lambda u\hackscore{k,k} \delta\hackscore{ij} + \mu (
37 u\hackscore{i,j} + u\hackscore{j,i})
38 \end{equation}
39 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 \begin{equation}
43 \rho u\hackscore{i,tt} = \sigma\hackscore{ij,j}
44 \end{equation}
45
46 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
55 \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
62 The first step is to choose an amplitude and create the source as in the
63 previous chapter.
64 \begin{python}
65 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 \end{python}
70 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 \begin{python}
76 src_dir=numpy.array([0.,-1.]) # defines direction of point source as down
77 y=y*src_dir
78 \end{python}
79 The function can then be applied as a boundary condition by setting it equal to
80 $y$ in the general form.
81 \begin{python}
82 mypde.setValue(y=y) #set the source as a function on the boundary
83 \end{python}
84 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 \begin{python}
87 # 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 u_m1=u
91 \end{python}
92
93 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 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[0]
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[0]-bleft; right=x[0]-bright; bottom=x[1]-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[0],tgamma)
192 fright=abc_bfunc(gright,bright,x[0],tgamma)
193 fbottom=abc_bfunc(gbottom,bbot,x[1],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

  ViewVC Help
Powered by ViewVC 1.1.26