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

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

Parent Directory Parent Directory | Revision Log Revision Log | View Patch Patch

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  \begin{equation} \label{eqn:wav} \index{wave equation}  \begin{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  \end{equation}  \end{equation}
34  where $\sigma$ represents stress and is given by:  where $\sigma$ is the stress given by:
35  \begin{equation} \label{eqn:sigw}  \begin{equation} \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  \end{equation}  \end{equation}
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  \begin{equation}  \begin{equation}
43   \rho u\hackscore{i,tt} = \sigma\hackscore{ij,j}  \rho u\hackscore{i,tt} = \sigma\hackscore{ij,j}
44  \end{equation}  \end{equation}
 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}  
  src_radius = 50  
 \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=U0*(cos(length(x-xc)*3.1415/src_radius)+1)*\  
               whereNegative(length(x-xc)-src_radius)*dunit  
  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}  \begin{equation} \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}  
 g=grad(u)  
 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;  
 \begin{equation}  
 u^n = 2u^{n-1}-u^{n-2}+h^2 \biggl(\frac{\partial ^2 u}{\partial t^2}\biggr)^n  
 \end{equation}  
 Substituting for the double time derivative we see;  
 \begin{equation}  
 \rho u^n = 2\rho u^{n-1}- \rho u^{n-2}+h^2 \sigma \hackscore{ij,j} ^n  
53  \end{equation}  \end{equation}
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

  ViewVC Help
Powered by ViewVC 1.1.26