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

revision 3053 by ahallam, Fri May 21 02:01:37 2010 UTC revision 3054 by ahallam, Wed Jun 30 02:22:25 2010 UTC
# Line 92  u_m1=u Line 92  u_m1=u
92
93  If the source will introduce energy to the system over a period longer than one  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  or two time steps (ie the initial conditions), $y$ can be updated during the
iteration stage.
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
162     y=
163
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

Legend:
 Removed from v.3053 changed lines Added in v.3054