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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 3054 - (hide annotations)
Wed Jun 30 02:22:25 2010 UTC (9 years, 4 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     % 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 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[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