doc/cookbook/example08.tex


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% Copyright (c) 2003-2010 by University of Queensland
% Earth Systems Science Computational Center (ESSCC)
% http://www.uq.edu.au/esscc
%
% Primary Business: Queensland, Australia
% Licensed under the Open Software License version 3.0
% http://www.opensource.org/licenses/osl-3.0.php
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

\section{Seismic Wave Propagation in Two Dimensions}
\editor{This chapter aims to expand the readers understanding of escript by
modelling the wave equations.
Challenges will include a second order differential (multiple initial
conditions). A new PDE to fit to the general form. Movement to a 3D problem
(maybe)??? }

\sslist{example08a.py}

We will now expand upon the previous chapter by introducing a vector form of
the wave equation. This means that the waves will have both a scalar magnitude,
but also a direction. This type of scenario is apparent in wave forms that
exhibit compressional and transverse particle motion. A common type of wave
that obeys this principle are seismic waves.

Wave propagation in the earth can be described by the elastic wave equation:
\begin{equation} \label{eqn:wav} \index{wave equation}
\rho \frac{\partial^{2}u\hackscore{i}}{\partial t^2} - \frac{\partial
\sigma\hackscore{ij}}{\partial x\hackscore{j}} = 0
\end{equation}
where $\sigma$ is the stress given by:
\begin{equation} \label{eqn:sigw}
 \sigma \hackscore{ij} = \lambda u\hackscore{k,k} \delta\hackscore{ij} + \mu (
u\hackscore{i,j} + u\hackscore{j,i})
\end{equation}
where $\lambda$ and $\mu$ are the Lame Coefficients. Specifically $\mu$ is the
bulk modulus. The \refEq{eqn:wav} can be written with the Einstein summation
convention as:
\begin{equation}
\rho u\hackscore{i,tt} = \sigma\hackscore{ij,j}
\end{equation}

In a similar process to the previous chapter, we will use the acceleration
solution to solve this PDE. By substituting $a$ directly for
$\frac{\partial^{2}u\hackscore{i}}{\partial t^2}$ we can derive the
displacement solution. Using $a$ \refEq{eqn:wav} becomes;
\begin{equation} \label{eqn:wava} 
\rho a\hackscore{i} - \frac{\partial
\sigma\hackscore{ij}}{\partial x\hackscore{j}} = 0
\end{equation}

\section{Vector source on the boundary}
For this particular example, we will introduce the source by applying a
displacment to the boundary during the initial time steps. The source will again
be
a radially propagating wave but due to the vector nature of the PDE used, a
direction will need to be applied to the source.

The first step is to choose an amplitude and create the source as in the
previous chapter. 
\begin{python}
 src_length = 20; print "src_length = ",src_length
# set initial values for first two time steps with source terms
y=U0*(cos(length(x-xc)*3.1415/src_length)+1)*whereNegative(length(x-xc)-src_leng
th)
\end{python}
where \verb xc  is the source point on the boundary of the model. The source
direction is then defined as an $(x,y)$ array and multiplied by the source
function. The directional array must have a magnitude of $1$ otherwise the
amplitude of the source will become modified. For this example, the source is
directed in the $-y$ direction.
\begin{python}
src_dir=numpy.array([0.,-1.]) # defines direction of point source as down
y=y*src_dir
\end{python}
The function can then be applied as a boundary condition by setting it equal to
$y$ in the general form.
\begin{python}
mypde.setValue(y=y) #set the source as a function on the boundary
\end{python}
Because we are no longer using the source to define our initial condition to
the model, we must set the model state to zero for the first two time steps.
\begin{python}
# initial value of displacement at point source is constant (U0=0.01)
# for first two time steps
u=[0.0,0.0]*whereNegative(x)
u_m1=u
\end{python}

If the source will introduce energy to the system over a period longer than one
or two time steps (ie the initial conditions), $y$ can be updated during the
iteration stage. 

\section{Time variant source}

\sslist{example08b.py}

Until this point, all of the wave propagation examples in this cookbook have
used impulsive sources which are smooth in space but not time. It is however,
advantageous to have a time smoothed source as it can reduce the temporal
frequency range and thus mitigate aliasing in the solution. 


It is quite 
simple to implement a source which is smooth in time. In addition to the
original source function the only extra requirement is a time function. For
this example it will be a decaying sinusoidal curve defined by;
\begin{python}
U0=0.1 # amplitude of point source
ls=100   # length of the source
source=np.zeros(ls,'float') # source array
g=np.log(0.01)/ls
for t in range(0,ls):
    source[t]=np.exp(g*t)*U0*np.sin(2.*np.pi*t/(0.75*ls))
\end{python}

We then build the source and the first two time steps via;
\begin{python}
# set initial values for first two time steps with source terms
y=source[0]
*(cos(length(x-xc)*3.1415/src_length)+1)*whereNegative(length(x-xc)-src_length)
src_dir=numpy.array([0.,-1.]) # defines direction of point source as down
y=y*src_dir
mypde.setValue(y=y) #set the source as a function on the boundary
# initial value of displacement at point source is constant (U0=0.01)
# for first two time steps
u=[0.0,0.0]*whereNegative(x)
u_m1=u
\end{python}

Finally, for the length of the source, we are required to update each new
solution in the itterative section of the solver. This is done via;
\begin{python}
# increment loop values
t=t+h; n=n+1
if (n < ls):
    y=source[n]**(cos(length(x-xc)*3.1415/src_length)+1)*\
                   whereNegative(length(x-xc)-src_length)
        y=y*src_dir; mypde.setValue(y=y) #set the source as a function on the
boundary
\end{python}

\section{Absorbing Boundary Conditions}
To mitigate the effect of the boundary on the model, absorbing boundary
conditions can be introduced. These conditions effectively dampen the wave
energy as they approach the bounday and thus prevent that energy from being
reflected. This type of approach is used typically when a model only represents
a small portion of the entire model, which in reality may have infinite bounds.
It is inpractical to calculate the solution for an infinite model and thus ABCs
allow us the create an approximate solution with small to zero boundary effects
on a model with a solvable size. 

To dampen the waves, the method of Cerjan(1985)
\footnote{\textit{A nonreflecting boundary condition for discrete acoustic and
elastic wave equations}, 1985, Cerjan C, Geophysics 50, doi:10.1190/1.1441945}
where the solution and the stress are multiplied by a damping function defined
on $n$ nodes of the domain adjacent to the boundary, given by;
\begin{equation}
 y=
\end{equation}
This is applied to the bounding 20-50 pts of the model using the location
specifiers of \esc;
\begin{python}
# Define where the boundary decay will be applied.
bn=30.
bleft=xstep*bn; bright=mx-(xstep*bn); bbot=my-(ystep*bn)
# btop=ystep*bn # don't apply to force boundary!!!

# locate these points in the domain
left=x[0]-bleft; right=x[0]-bright; bottom=x[1]-bbot

tgamma=0.98   # decay value for exponential function
def calc_gamma(G,npts):
    func=np.sqrt(abs(-1.*np.log(G)/(npts**2.)))
    return func

gleft  = calc_gamma(tgamma,bleft)
gright = calc_gamma(tgamma,bleft)
gbottom= calc_gamma(tgamma,ystep*bn)

print 'gamma', gleft,gright,gbottom

# calculate decay functions
def abc_bfunc(gamma,loc,x,G):
    func=exp(-1.*(gamma*abs(loc-x))**2.)
    return func

fleft=abc_bfunc(gleft,bleft,x[0],tgamma)
fright=abc_bfunc(gright,bright,x[0],tgamma)
fbottom=abc_bfunc(gbottom,bbot,x[1],tgamma)
# apply these functions only where relevant
abcleft=fleft*whereNegative(left)
abcright=fright*wherePositive(right)
abcbottom=fbottom*wherePositive(bottom)
# make sure the inside of the abc is value 1
abcleft=abcleft+whereZero(abcleft)
abcright=abcright+whereZero(abcright)
abcbottom=abcbottom+whereZero(abcbottom)
# multiply the conditions together to get a smooth result
abc=abcleft*abcright*abcbottom
\end{python}
Note that the boundary conditions are not applied to the surface, as this is
effectively a free surface where normal reflections would be experienced. The
resulting boundary damping function can be viewed in \ref{fig:abconds}

%\begin{figure}{ht}
 %\centering
 %\includegraphics[width=5in]{figures/abconds.png}
%\end{figure}


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(gt)U0np.sin(2.np.pit/(0.75ls))
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=xstepbn; bright=mx-(xstepbn); 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.(gammaabs(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=abcleftabcrightabcbottom
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