# Annotation of /trunk/doc/cookbook/twodswp001.tex

Revision 2434 - (hide annotations)
Thu May 21 07:59:53 2009 UTC (11 years, 2 months ago) by ahallam
Original Path: trunk/doc/cookbook/TEXT/twodswp001.tex
File MIME type: application/x-tex
File size: 7043 byte(s)
Week 8: Still working on seismic wave propagation. Need to come up with a smooth source function next week. All TEXT files up to date.

 1 ahallam 2411 2 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 3 % 4 % Copyright (c) 2003-2009 by University of Queensland 5 % Earth Systems Science Computational Center (ESSCC) 6 7 % 8 % Primary Business: Queensland, Australia 9 % Licensed under the Open Software License version 3.0 10 11 % 12 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 13 14 ahallam 2426 \section{Seismic Wave Propagation in Two Dimensions} 15 \editor{This chapter aims to expand the readers understanding of escript by modelling the wave equations. 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)??? } 17 18 ahallam 2434 \verb|wavesolver2d.py| 19 ahallam 2426 20 ahallam 2434 Wave propagation in the earth can be described by the wave equation: 21 \begin{equation} \label{eqn:wav} \index{wave equation} 22 \rho \frac{\partial^{2}u\hackscore{i}}{\partial t^2} - \frac{\partial \sigma\hackscore{ij}}{\partial\hackscore{j}} = 0 23 \end{equation} 24 where $\sigma$ represents stress and is given by: 25 \begin{equation} \label{eqn:sigw} 26 \sigma \hackscore{ij} = \lambda u\hackscore{k,k} \delta\hackscore{ij} + \mu ( u\hackscore{i,j} + u\hackscore{j,i}) 27 \end{equation} 28 $\lambda$ and $\mu$ are the Lame Coefficients. Specifically $\mu$ is the bulk modulus. Equation \ref{eqn:wav} can be written with the Einstein summation convention as: 29 \begin{equation} 30 \rho u\hackscore{i,tt} = \sigma\hackscore{ij,j} 31 \end{equation} 32 For this model we will specify the boundary conditions such that the normal of the stress from the boundary is zero. 33 \begin{eqnarray} \label{eqn:bdw} 34 \sigma \hackscore{ij}n\hackscore{j}=0 35 \end{eqnarray} 36 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. 37 38 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; 39 \begin{verbatim} 40 domain : domain to solve over 41 h : time step 42 tend : end time 43 lam, mu : lames constants for domain 44 rho : density of domain 45 U0 : magnitude of source 46 xc : source location in domain (Vector) 47 savepath: where to output the data files 48 \end{verbatim} 49 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; 50 \begin{verbatim} 51 mypde.setSolverMethod(LinearPDE.LUMPING) 52 \end{verbatim} 53 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 with however as this function can only be used when the $A$, $B$ and$C$ coefficients of the general form are zero. 54 55 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 specified individually. It should be noted here that if multiple time steps are known or understood in the begining of a modle, they can be added to the simulation manually. The solver is then able to continue the model from where the data ends. Alternatively, if the source motion is understood, its position can be corrected for each itteration to create a more accurate recreation of an event. The source in this example will induce a radially propagating wave. To do this we will introduce a small displacement 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 is defined by 0.1 times the maxium (\verb Lsup() ) size of the domain; 56 \begin{verbatim} 57 src_radius = 0.1*Lsup(domain.getSize()) 58 \end{verbatim} 59 60 61 \esc can be used to model the propgation of waves through a medium. 62 63 64 65 66 67 In this example we will see how the wave equation can be implemented using \esc and solved for in two dimensions. Our domain is defined by a thin sheet that has dimensions $x$ and $y$ and to model waves we will introduce a point source displacement at time zero. The affects of this displacement should propagate radially from the source and eventually be reflected from the boundaries of the model. 68 69 70 71 72 73 74 ahallam 2426 The code described in this section can be found in \fileex{wavesolver2d001.py} 75 In a similar manner to the previous chapter the first step to creating our script is to import the necessary modules and functions. 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. This example includes an acknowledgement clause 76 \begin{verbatim} 77 #Check to make sure number of time steps is not too large. 78 print "Time step size= ",h, "Expected number of outputs= ",tend/h 79 proceeder = raw_input("Is this ok?(y/n)") 80 #Exit if user thinks too many outputs. 81 if proceeder == "n": 82 sys.exit() 83 \end{verbatim} 84 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. 85 86 To solve this PDE we are going to introduce the concept of a python library. A library is useful as it allows a user to store defined functions that can be called to solve generic problems. The 2D wave equation satisfies this criteria. The first step is to create a new python file which we have called \verb cblib.py within this file we can set all of the necessary includes to make things easier in the future. Other advantages of libraries include a reduction in the duplication of code and the ability to modularise functions and variables. 87 88 89