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

Revision 3029 - (show annotations)
Fri May 21 02:01:37 2010 UTC (10 years, 1 month ago) by ahallam
File MIME type: application/x-tex
File size: 4008 byte(s)
small updates and work on seismic code

 1 2 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 3 % 4 % Copyright (c) 2003-2010 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 \section{Seismic Wave Propagation in Two Dimensions} 15 \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 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 Wave propagation in the earth can be described by the elastic wave equation: 30 \begin{equation} \label{eqn:wav} \index{wave equation} 31 \rho \frac{\partial^{2}u\hackscore{i}}{\partial t^2} - \frac{\partial 32 \sigma\hackscore{ij}}{\partial x\hackscore{j}} = 0 33 \end{equation} 34 where $\sigma$ is the stress given by: 35 \begin{equation} \label{eqn:sigw} 36 \sigma \hackscore{ij} = \lambda u\hackscore{k,k} \delta\hackscore{ij} + \mu ( 37 u\hackscore{i,j} + u\hackscore{j,i}) 38 \end{equation} 39 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} 43 \rho u\hackscore{i,tt} = \sigma\hackscore{ij,j} 44 \end{equation} 45 46 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 55 \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 62 The first step is to choose an amplitude and create the source as in the 63 previous chapter. 64 \begin{python} 65 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 \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 91 \end{python} 92 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.