ViewVC logotype

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

Parent Directory Parent Directory | Revision Log Revision Log

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
2 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
3 %
4 % Copyright (c) 2003-2010 by University of Queensland
5 % 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 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
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)??? }
21 \sslist{example08a.py}
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.
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}
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}
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.
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}
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.

  ViewVC Help
Powered by ViewVC 1.1.26