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

trunk/doc/cookbook/twodswp001.tex revision 3025 by ahallam, Thu May 6 01:20:46 2010 UTC trunk/doc/cookbook/example08.tex revision 3029 by ahallam, Fri May 21 02:01:37 2010 UTC
# Line 12  Line 12
12  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
13
14  \section{Seismic Wave Propagation in Two Dimensions}  \section{Seismic Wave Propagation in Two Dimensions}
15  \editor{This chapter aims to expand the readers understanding of escript by modelling the wave equations.  \editor{This chapter aims to expand the readers understanding of escript by
16  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)??? }  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  \verb|wavesolver2d.py|  Wave propagation in the earth can be described by the elastic wave equation:

Wave propagation in the earth can be described by the wave equation:
30   \label{eqn:wav} \index{wave equation}   \label{eqn:wav} \index{wave equation}
31  \rho \frac{\partial^{2}u\hackscore{i}}{\partial t^2} - \frac{\partial \sigma\hackscore{ij}}{\partial\hackscore{j}} = 0  \rho \frac{\partial^{2}u\hackscore{i}}{\partial t^2} - \frac{\partial
32    \sigma\hackscore{ij}}{\partial x\hackscore{j}} = 0
33
34  where $\sigma$ represents stress and is given by:  where $\sigma$ is the stress given by:
35   \label{eqn:sigw}   \label{eqn:sigw}
36   \sigma \hackscore{ij} = \lambda u\hackscore{k,k} \delta\hackscore{ij} + \mu ( u\hackscore{i,j} + u\hackscore{j,i})   \sigma \hackscore{ij} = \lambda u\hackscore{k,k} \delta\hackscore{ij} + \mu (
37    u\hackscore{i,j} + u\hackscore{j,i})
38
39  $\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:  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
43   \rho u\hackscore{i,tt} = \sigma\hackscore{ij,j}  \rho u\hackscore{i,tt} = \sigma\hackscore{ij,j}
44
For this model we will specify the boundary conditions such that the normal of the stress from the boundary is zero.
\begin{eqnarray} \label{eqn:bdw}
\sigma \hackscore{ij}n\hackscore{j}=0
\end{eqnarray}
To solve this PDE we are going to write a more generic solution routine than we have in the previous chapters. There are a number of advantages to this approach. Firstly, by writing a subroutine that will solve a 2D wave propagation problem it reduces the amount of code duplication that may occur. When errors arrise one need only ammend the subroutine rather than all iterations of it that may have been created. This saves time and effort in the long run.

To create the our subroutine we will need to import all our necessary libraries again. This is as for previous examples. Then we can define our script and the variables it will take. Our subroutine is located in \verb|/examples/cblib/wavesolver2d.py|  . The arguments of the subroutine are;
\begin{python}
domain  : domain to solve over
h       : time step
tend    : end time
lam, mu : lames constants for domain
rho : density of domain
U0  : magnitude of source
xc  : source location in domain (Vector)
savepath: where to output the data files
\end{python}
There are a few differences between the wave equation and the heat diffusion
problem of the previous chapter. While the nodes are defined the same way with
the function \verb getX  and the PDE is still linear; one must consider the
solution method. Without the line;
\begin{python}
mypde.setSolverMethod(LinearPDE.LUMPING)
\end{python}
the PDE would take a significant amount of time to solve.
The \verb LUMPING functionality implements an aggressive approximation for the
$D$ coefficient matrix of the \esc linear PDE general form.
While \verb LUMPING introduces additional error to the solution it can
significantly reduce the solution time. Care should be taken however, as this
function can only be used when the $A$, $B$ and $C$ coefficients of the general
form are zero.

As the wave equation has a double time derivative, it is not sufficient to only
stipulate the initial conditions for one time step. Two time steps must be
specified so that the equation can be solved. For this example $u$ (\verb u )
and $u(t-1)$ (\verb u_m1 ) will be the same but if both of these condititions
are known, they can be specified individually. It should be noted here that if
multiple time steps are known at the begining of a model, they can be added to
the simulation manually. The solver can then continue the model from the end of
the known data. Alternatively, if the source motion is understood, its position
can be corrected for each iteration to create a more accurate recreation of an
event.

The source in this example will induce a radially propagating wave. A small
displacement will be applied to the medium about a singularity which we have
called \verb xc  , this is the source location. We start by giving the source
some spatial magnitude by defining a small radius about \verb xc  which is
affected. The \verb src_radius  needs to cover a significant portion of grid
nodes, otherwise the waves generated will suffer from dispersion due to an
inadequate grid step size. If the source is small, the grid stepping must
reflect the size of the source for more accurate results. Our radius will be;
\begin{python}
\end{python}
Now that the extent of the source has been allocated it needs two more things; a
direction and a magnitude. We can choose a direction based on the 360 degrees
that exist in a full circle. If we take $\theta=0$ to be the x-axis and move
counter clockwise then we can create a directional vector $U=[dx,dy]$  where
$tan(\theta) = dy/dx$. It is also necessary to ensure that our directional
vector of unit length such that $|U|=1$; which implies $\sqrt{dx^2+dy^2}=1$. By
doing this we ensure that no accidental scaling is introduced to our source
term. Here are three examples of different directions which satisfy the above
conditions;
\begin{enumerate}
\item Along the x-axis: $U=[dx=1,dy=0]$
\item Along the y-axis: $U=[dx=0,dy=1]$
\item At 45deg: $U=[dx=\frac{1}{\sqrt2},dy=\frac{1}{\sqrt2}]$
\end{enumerate}
There are limitations to specifying the source in this manner. Realistically we
would not expect a 2D surface source to move form side to side as an isotropic
source makes more sense. \editor{I am not sure here how to create an isotropic
source function.}. In the 3D case things are not quite so bad. Normally we are
interested in the p-waves that are directed dowwards and thus we need not have
any x or y component in our source directionality. This still introduced
assumptions and removes realistic wave motions both P and S from the model.
For our example we will use;
\begin{python}
dunit=numarray.array([1.,0.])
\end{python}
Next we must define the values of our entire domain for the first and second
time step. For the purposes of this example it is sufficient to have these two
timesteps as equal. Setting the source is similar to earlier problems where we
can use \esc functions to set specific areas of the domain to certain values. We
must also smooth our sourse to its surrounds to prevent ?diffusion? errors. This
is acheived by using a cosine taper. Our source terms then become;
\begin{python}
u_m1=u
\end{python}
45
46  It is often useful to know the values of PDE at certain locations in the model.  In a similar process to the previous chapter, we will use the acceleration
47  To acheive this we are going to use a new generic function  solution to solve this PDE. By substituting $a$ directly for
48  called \verb cbphones  which allows us to specify receiver locations to record  $\frac{\partial^{2}u\hackscore{i}}{\partial t^2}$ we can derive the
49  the PDE values at those points. The function \verb cbphones as the arguments;  displacement solution. Using $a$ \refEq{eqn:wav} becomes;
50  \begin{python}   \label{eqn:wava}
51  #   domain  : domain of model  \rho a\hackscore{i} - \frac{\partial
52  #   U       : Current time state displacement solution.  \sigma\hackscore{ij}}{\partial x\hackscore{j}} = 0
#   phones  : Geophone Locations
#   dim     : model dimesions
#   savepath: where to output the data files local is default
\end{python}
\editor{not generic as of yet but may move to make cbphones and the phones
positioning a serious part of wavesolver 2d}
We have chosen to have three receivers and they are called using;
\begin{python}
u_pot = cbphones(domain,u,[[0,500],[250,500],[400,500]],2)
\end{python}
This places the receivers on the surface at the source location and two
locations further along the top of the model. The output \verb u_pot  can then
be split and saved to file using the following command;
\begin{python}
u_pc_data=open(os.path.join(savepath,'U_pc.out'),'w')
u_pc_data.write("%f %f %f %f %f %f
%f\n"%(t,u_pc_x1,u_pc_y1,u_pc_x2,u_pc_y2,u_pc_x3,u_pc_y3))
\end{python}
Convieniently this saves the time, x direction displacement and y direction
displacement values for these locations. Now that the initial conditions have
been defined we can tackle the task of solving the wave equation for the number
of required time steps. To do this we require a while loop and form of the wave
equation which fits our general linear PDE form. We start with the form of the
equation for stress \refEq{eqn:sigw}. We can define the kronecker matrix using
the domain and take the derivative of \verb u  via the function \verb|grad(u)|
. As $\lambda$ and $\mu$ are constants we can now define $\sigma$;
\begin{python}
stress=lam*trace(g)*kmat+mu*(g+transpose(g))
\end{python}
Solving for the double time derivative of u on the LHS of \refEq{eqn:wav}
required us to use the centred difference forumlua which returns;

u^n = 2u^{n-1}-u^{n-2}+h^2 \biggl(\frac{\partial ^2 u}{\partial t^2}\biggr)^n

Substituting for the double time derivative we see;

\rho u^n = 2\rho u^{n-1}- \rho u^{n-2}+h^2 \sigma \hackscore{ij,j} ^n
53
54  This fits the general form $Du=-X \hackscore{j,j} + Y$ where $D=rho$; $Y=2\rho 55 u^{n-1}- \rho u^{n-2}$ and $X=-h^2 \sigma \hackscore{ij,j} ^n$. \verb D does not  \section{Vector source on the boundary}
56  vary between time steps can be defined before our iteration loop via;  For this particular example, we will introduce the source by applying a
57  \begin{python}  displacment to the boundary during the initial time steps. The source will again
58  mypde.setValue(D=kmat*rho)  be
59  \end{python}  a radially propagating wave but due to the vector nature of the PDE used, a
60  The values for \verb u  must be refreshed after each iteration and are thus  direction will need to be applied to the source.
61  defined within our while loop via;
62  \begin{python}  The first step is to choose an amplitude and create the source as in the
63  mypde.setValue(X=-stress*(h*h),Y=(rho*2*u-rho*u_m1))  previous chapter.
64  \end{python}  \begin{python}
65  With each iteration we update \verb u  and migrate our answers into the correct   src_length = 20; print "src_length = ",src_length
66  variables. The iterative values must also be updated as well as the response  # set initial values for first two time steps with source terms
67  from our receiver locations. All this is acheived via;  y=U0*(cos(length(x-xc)*3.1415/src_length)+1)*whereNegative(length(x-xc)-src_leng
68  \begin{python}  th)
69  u_p1 = mypde.getSolution()  \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  u_m1=u
u=u_p1
t+=h
n+=1
u_pot = cbphones(domain,u,[[125.,250.],[250.,250.],[250.,375.]],2)
# save displacements at point source to file for t > 0
u_pc_data.write("%f %f %f %f %f %f %f\n"%(t,u_pc_x1,u_pc_y1,u_pc_x2,u_pc_y2,u_pc_x3,u_pc_y3))
\end{python}
With an appropriate file saving output we now have a working generic solver for
our 2D wave equation problem. We have also included our two new generic programs
in the \verb cblib  library so they can be more simply imported along with their
own dependancies to our test script.

Writing a test program will allow us to more easily pass the variables required
to the solver to generate an output solution. Our testing code described in this
section can be found in \fileex{wavesolver2d001.py}. In a similar manner to the
previous chapter the first step to creating our script is to import the
necessary modules and functions including our new library file. Following this
the PDE and control variables must be defined. This includes the domain
dimensions and type, the time scale and the time step. To ensure stability the
time step can be calcuated such that it satisfies the Courant stability criteria
\editor{MORE HERE ONCE METHOD FINALISED}. Considering the complexity of the
computational solution to the wave equation it is proudant to consider how many
steps will need to be solved. Our test script thus includes an acknowledgement
clause;
\begin{python}
#Check to make sure number of time steps is not too large.
print "Time step size= ",h, "Expected number of outputs= ",tend/h
proceeder = raw_input("Is this ok?(y/n)")
#Exit if user thinks too many outputs.
if proceeder == "n":
sys.exit()
91  \end{python}  \end{python}
92  This requires that the user knows the number of itterations that will be required to solve the model for the time period \verb 0  to \verb tend . The command \verb sys.exit()  is used here to halt the script if the input to preceeder is \verb n  and thus prevent a forced crash of the script should its projected solve time be too large.
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.

Legend:
 Removed from v.3025 changed lines Added in v.3029