/[escript]/trunk/doc/user/wave.tex
ViewVC logotype

Annotation of /trunk/doc/user/wave.tex

Parent Directory Parent Directory | Revision Log Revision Log


Revision 2104 - (hide annotations)
Thu Nov 27 05:36:47 2008 UTC (13 years, 8 months ago) by jfenwick
File MIME type: application/x-tex
File size: 14972 byte(s)
Fixed non-compilage.
Removed mysterious t and t! from output.
1 ksteube 1811
2     %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
3 ksteube 1316 %
4 ksteube 1811 % Copyright (c) 2003-2008 by University of Queensland
5     % Earth Systems Science Computational Center (ESSCC)
6     % http://www.uq.edu.au/esscc
7 gross 625 %
8 ksteube 1811 % 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 gross 625 %
12 ksteube 1811 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
13 ksteube 1316
14 ksteube 1811
15 ksteube 1318 \section{3-D Wave Propagation}
16 jgs 107 \label{WAVE CHAP}
17    
18 artak 1971 In this next example we want to calculate the displacement field $u\hackscore{i}$ for any time $t>0$ by solving the wave equation:
19     \index{wave equation}
20 jgs 107 \begin{eqnarray}\label{WAVE general problem}
21     \rho u\hackscore{i,tt} - \sigma\hackscore{ij,j}=0
22     \end{eqnarray}
23     in a three dimensional block of length $L$ in $x\hackscore{0}$
24     and $x\hackscore{1}$ direction and height $H$
25     in $x\hackscore{2}$ direction. $\rho$ is the known density which may be a function of its location.
26     $\sigma\hackscore{ij}$ is the stress field \index{stress} which in case of an isotropic, linear elastic material is given by
27     \begin{eqnarray} \label{WAVE stress}
28     \sigma\hackscore{ij} & = & \lambda u\hackscore{k,k} \delta\hackscore{ij} + \mu ( u\hackscore{i,j} + u\hackscore{j,i})
29     \end{eqnarray}
30     where $\lambda$ and $\mu$ are the Lame coefficients
31     \index{Lame coefficients} and $\delta\hackscore{ij}$ denotes the Kronecker symbol\index{Kronecker symbol}.
32     On the boundary the normal stress is given by
33     \begin{eqnarray} \label{WAVE natural}
34     \sigma\hackscore{ij}n\hackscore{j}=0
35     \end{eqnarray}
36 lkettle 581 for all time $t>0$.
37    
38     At initial time $t=0$ the displacement
39     $u\hackscore{i}$ and the velocity $u\hackscore{i,t}$ are given:
40 jgs 107 \begin{eqnarray} \label{WAVE initial}
41 lkettle 581 u\hackscore{i}(0,x)= \left\{\begin{array}{cc} U\hackscore{0} & \mbox{for $x$ at point charge, $x\hackscore{C}$} \\ 0 & \mbox{elsewhere}\end{array} \right. & \mbox{ and } & u\hackscore{i,t}(0,x)=0 \;
42 jgs 107 \end{eqnarray}
43     for all $x$ in the domain.
44    
45 lkettle 581 Here we are modelling a point source at the point $x\hackscore C$, in the numerical solution we
46     set the initial displacement to be $U\hackscore 0$ in a sphere of small radius around the point
47     $x\hackscore C$.
48 jgs 107
49 ksteube 1316 We use an explicit time integration scheme to calculate the displacement field $u$ at
50 jgs 107 certain time marks $t^{(n)}$ where $t^{(n)}=t^{(n-1)}+h$ with time step size $h>0$. In the following the upper index ${(n)}$ refers to values at time $t^{(n)}$. We use the Verlet scheme \index{Verlet scheme} with constant time step size $h$
51     which is defined by
52     \begin{eqnarray} \label{WAVE dyn 2}
53 lkettle 575 u^{(n)}=2u^{(n-1)}-u^{(n-2)} + h^2 a^{(n)} \\
54 jgs 107 \end{eqnarray}
55 lkettle 581 for all $n=2,3,\ldots$. It is designed to solve a system of equations of the form
56 jgs 107 \begin{eqnarray} \label{WAVE dyn 2b}
57     u\hackscore{,tt}=G(u)
58     \end{eqnarray}
59 gross 660 where one sets $a^{(n)}=G(u^{(n-1)})$.
60 jgs 107
61     In our case $a^{(n)}$ is given by
62     \begin{eqnarray}\label{WAVE dyn 3}
63     \rho a^{(n)}\hackscore{i}=\sigma^{(n-1)}\hackscore{ij,j}
64     \end{eqnarray}
65     and boundary conditions
66 gross 593 \begin{eqnarray} \label{WAVE natural at n}
67 jgs 107 \sigma^{(n-1)}\hackscore{ij}n\hackscore{j}=0
68     \end{eqnarray}
69     derived from \eqn{WAVE natural} where
70     \begin{eqnarray} \label{WAVE dyn 3a}
71 lkettle 581 \sigma\hackscore{ij}^{(n-1)} & = & \lambda u^{(n-1)}\hackscore{k,k} \delta\hackscore{ij} + \mu ( u^{(n-1)}\hackscore{i,j} + u^{(n-1)}\hackscore{j,i}).
72 jgs 107 \end{eqnarray}
73 lkettle 581 Now we have converted our problem for displacement, $u^{(n)}$, into a problem for
74     acceleration, $a^(n)$, which now depends
75     on the solution at the previous two time steps, $u^{(n-1)}$ and $u^{(n-2)}$.
76 jgs 107
77 ksteube 1316 In each time step we have to solve this problem to get the acceleration $a^{(n)}$, and we will
78 gross 565 use the \LinearPDE class of the \linearPDEs to do so. The general form of the PDE defined through
79 ksteube 1316 the \LinearPDE class is discussed in \Sec{SEC LinearPDE}. The form which is relevant here is
80 jgs 107 \begin{equation}\label{WAVE dyn 100}
81 lkettle 581 D\hackscore{ij} a^{(n)}\hackscore{j} = - X\hackscore{ij,j}\; .
82 jgs 107 \end{equation}
83 ksteube 1316 The natural boundary condition
84 jgs 107 \begin{equation}\label{WAVE dyn 101}
85     n\hackscore{j}X\hackscore{ij} =0
86     \end{equation}
87 ksteube 1316 is used.
88 jgs 107
89 ksteube 1316 With $u=a^{(n)}$ we can identify the values to be assigned to $D$ and $X$:
90 jgs 107 \begin{equation}\label{WAVE dyn 6}
91     \begin{array}{ll}
92     D\hackscore{ij}=\rho \delta\hackscore{ij}&
93     X\hackscore{ij}=-\sigma^{(n-1)}\hackscore{ij} \\
94     \end{array}
95     \end{equation}
96 lkettle 581
97    
98     %Remember that \var{y.whereZero()} returns a function which is $1$ for those points $x$ where $y$ is zero
99     %(up to a certain tolerance, see \Sec{SEC ESCRIPT DATA}) and $0$ for those points where $y$ is not equal to zero.
100 jgs 107
101 lkettle 581 The following script defines a the function \function{wavePropagation} which
102 jgs 107 implements the Verlet scheme to solve our wave propagation problem.
103     The argument \var{domain} which is a \Domain class object
104 ksteube 1316 defines the domain of the problem. \var{h} and \var{tend} are the time step size
105     and the end time of the simulation. \var{lam}, \var{mu} and
106 lkettle 581 \var{rho} are material properties.
107 jgs 107 \begin{python}
108     from esys.linearPDEs import LinearPDE
109 lkettle 581 from numarray import identity,zeros,ones
110 jgs 107 from esys.escript import *
111 lkettle 581 from esys.escript.pdetools import Locator
112     def wavePropagation(domain,h,tend,lam,mu,rho,U0):
113 jgs 107 x=domain.getX()
114     # ... open new PDE ...
115     mypde=LinearPDE(domain)
116 lkettle 581 mypde.setSolverMethod(LinearPDE.LUMPING)
117 jgs 110 kronecker=identity(mypde.getDim())
118 lkettle 581 # spherical source at middle of bottom face
119     xc=[width/2.,width/2.,0.]
120     # define small radius around point xc
121     # Lsup(x) returns the maximum value of the argument x
122     src_radius = 0.1*Lsup(domain.getSize())
123     print "src_radius = ",src_radius
124     dunit=numarray.array([1.,0.,0.]) # defines direction of point source
125     mypde.setValue(D=kronecker*rho)
126 jgs 107 # ... set initial values ....
127 jgs 110 n=0
128 lkettle 581 # initial value of displacement at point source is constant (U0=0.01)
129     # for first two time steps
130     u=U0*whereNegative(length(x-xc)-src_radius)*dunit
131     u_last=U0*whereNegative(length(x-xc)-src_radius)*dunit
132 jgs 107 t=0
133 lkettle 581 # define the location of the point source
134     L=Locator(domain,xc)
135     # find potential at point source
136     u_pc=L.getValue(u)
137     print "u at point charge=",u_pc
138     u_pc_x = u_pc[0]
139     u_pc_y = u_pc[1]
140     u_pc_z = u_pc[2]
141     # open file to save displacement at point source
142     u_pc_data=open('./data/U_pc.out','w')
143     u_pc_data.write("%f %f %f %f\n"%(t,u_pc_x,u_pc_y,u_pc_z))
144 jgs 107 while t<tend:
145     # ... get current stress ....
146     g=grad(u)
147 jgs 110 stress=lam*trace(g)*kronecker+mu*(g+transpose(g))
148     # ... get new acceleration ....
149 lkettle 581 mypde.setValue(X=-stress)
150 jgs 107 a=mypde.getSolution()
151     # ... get new displacement ...
152     u_new=2*u-u_last+h**2*a
153     # ... shift displacements ....
154     u_last=u
155     u=u_new
156     t+=h
157 jgs 110 n+=1
158     print n,"-th time step t ",t
159 lkettle 581 L=Locator(domain,xc)
160     u_pc=L.getValue(u)
161     print "u at point charge=",u_pc
162     u_pc_x=u_pc[0]
163     u_pc_y=u_pc[1]
164     u_pc_z=u_pc[2]
165     # save displacements at point source to file for t > 0
166     u_pc_data.write("%f %f %f %f\n"%(t,u_pc_x,u_pc_y,u_pc_z))
167     # ... save current acceleration in units of gravity and displacements
168     if n==1 or n%10==0: saveVTK("./data/usoln.%i.vtu"%(n/10),acceleration=length(a)/9.81,
169     displacement = length(u), Ux = u[0] )
170     u_pc_data.close()
171 jgs 107 \end{python}
172 jgs 110 Notice that
173 ksteube 1316 all coefficients of the PDE which are independent of time $t$ are set outside the \code{while}
174     loop. This is very efficient since it allows the \LinearPDE class to reuse information as much as possible
175 jgs 110 when iterating over time.
176    
177 ksteube 1316 there are a few new \escript functions in this example:
178 jgs 107 \code{grad(u)} returns the gradient $u\hackscore{i,j}$ of $u$ (in fact \var{grad(g)[i,j]}=$=u\hackscore{i,j}$).
179 ksteube 1316 There are restrictions on the argument of the \function{grad} function, for instance
180     the statement \code{grad(grad(u))} will raise an exception.
181 jgs 107 \code{trace(g)} returns the sum of the main diagonal elements \var{g[k,k]} of \var{g}
182 ksteube 1316 and \code{transpose(g)} returns the matrix transpose of \var{g} (ie. $\var{transpose(g)[i,j]}=\var{g[j,i]}$).
183 jgs 107
184 lkettle 581 We initialize the values of \code{u} and \code{u_last} to be $U\hackscore 0$ for a small
185     sphere of radius, \code{src_radius} around the point source located at $x\hackscore C$.
186     These satisfy the initial conditions defined in \eqn{WAVE initial}.
187 jgs 107
188 lkettle 581 The output \code{U_pc.out} and \code{vtu} files are saved in a directory called \code{data}.
189     The \code{U_pc.out} stores 4 columns of data: $t,u\hackscore x,u\hackscore y,u\hackscore z$
190     respectively, where $u\hackscore x,u\hackscore y,u\hackscore z$ are the $x\hackscore 0,x\hackscore 1,x\hackscore 2$ components of
191     the displacement vector $u$ at the point source. These can be
192     plotted easily using any plotting package. In gnuplot the command:
193     \code{plot 'U_pc.out' u 1:2 title 'U_x' w l lw 2, 'U_pc.out' u 1:3 title 'U_y' w l lw 2,
194     'U_pc.out' u 1:4 title 'U_z' w l lw 2} will reproduce Figure \ref{WAVE FIG 1}.
195     %\begin{python}
196     %u=Vector(0.,ContinuousFunction(domain))
197     %\end{python}
198     %declares \var{u} as a vector-valued function of its coordinates
199     %in the domain (i.e. with two or three components in the case of
200 ksteube 1316 %a two or three dimensional domain, respectively). The first argument defines the value of the function which is equal
201 lkettle 581 %to $0$ for all vectot components. The second argument \var{ContinuousFunction(domain)}
202     %specifies the `smoothness` of the function. In our case we want the displacement field to be
203     %a continuous function over the \Domain \var{domain}. On other option would be
204 ksteube 1316 %\var{Function(domain)} in which case continuity is not guaranteed. In fact the
205 lkettle 581 %argument of the \function{grad} function has to be of the typ \var{ContinuousFunction(domain)} while
206     %returned function is of the var{Function(domain)} type. For more details see \Sec{SEC ESCRIPT DATA}.
207    
208 jgs 110 The statement \code{myPDE.setLumpingOn()}
209 ksteube 1316 switches on the use of an aggressive approximation of the PDE operator as a diagonal matrix
210     formed from the coefficient \var{D}.
211     The approximation allows, at the cost of
212     additional error, very fast
213     solution of the PDE. When using \code{setLumpingOn} the presence of \var{A}, \var{B} or \var
214     {C} will produce wrong results.
215 jgs 110
216 jgs 107
217 jgs 110 One of the big advantages of the Verlet scheme is the fact that the problem to be solved
218 ksteube 1316 in each time step is very simple and does not involve any spatial derivatives (which is what allows us to use lumping in this simulation).
219 lkettle 581 The problem becomes so simple because we use the stress from the last time step rather then the stress which is
220     actually present at the current time step. Schemes using this approach are called an explicit time integration
221 ksteube 1316 schemes \index{explicit scheme} \index{time integration!explicit}. The
222 jgs 110 backward Euler scheme we have used in \Chap{DIFFUSION CHAP} is
223 ksteube 1316 an example of an implicit scheme
224 jgs 110 \index{implicit scheme} \index{time integration!implicit}. In this case one uses the actual status of
225 lkettle 581 each variable at a particular time rather then values from previous time steps. This will lead to a problem which is
226 jgs 110 more expensive to solve, in particular for non-linear problems.
227     Although
228 lkettle 581 explicit time integration schemes are cheap to finalize a single time step, they need significantly smaller time
229     steps then implicit schemes and can suffer from stability problems. Therefore they need a
230 jgs 110 very careful selection of the time step size $h$.
231 jgs 107
232 jgs 110 An easy, heuristic way of choosing an appropriate
233 jgs 107 time step size is the Courant condition \index{Courant condition} \index{explicit scheme!Courant condition}
234     which says that within a time step a information should not travel further than a cell used in the
235     discretization scheme. In the case of the wave equation the velocity of a (p-) wave is given as
236     $\sqrt{\frac{\lambda+2\mu}{\rho}}$ so one should choose $h$ from
237     \begin{eqnarray}\label{WAVE dyn 66}
238     h= \frac{1}{5} \sqrt{\frac{\rho}{\lambda+2\mu}} \Delta x
239     \end{eqnarray}
240 lkettle 581 where $\Delta x$ is the cell diameter. The factor $\frac{1}{5}$ is a safety factor considering the heuristics of
241 jgs 107 the formula.
242    
243 lkettle 581 The following script uses the \code{wavePropagation} function to solve the
244     wave equation for a point source located at the bottom face of a block. The width of the block in
245     each direction on the bottom face is $10\mbox{km}$ ($x\hackscore 0$ and $x\hackscore 1$ directions ie. \code{l0} and \code{l1}).
246 artak 1971 The \code{ne} gives the number of elements in $x\hackscore{0}$ and $x\hackscore 1$ directions.
247 lkettle 581 The depth of the block is aligned with the $x\hackscore{2}$-direction.
248     The depth (\code{l2}) of the block in
249     the $x\hackscore{2}$-direction is chosen so that there are 10 elements and the magnitude of the
250     of the depth is chosen such that the elements
251     become cubic. We chose 10 for the number of elements in $x\hackscore{2}$-direction so that the
252     computation would be faster. \code{Brick($n\hackscore 0, n\hackscore 1, n\hackscore 2,l\hackscore 0,l\hackscore 1,l\hackscore 2$)} is an \finley function which creates a rectangular mesh
253     with $n\hackscore 0 \times n\hackscore 1 \times n\hackscore2$ elements over the brick $[0,l\hackscore 0] \times [0,l\hackscore 1] \times [0,l\hackscore 2]$.
254 jgs 107 \begin{python}
255 lkettle 581 from esys.finley import Brick
256     ne=32 # number of cells in x_0 and x_1 directions
257     width=10000. # length in x_0 and x_1 directions
258 jgs 110 lam=3.462e9
259     mu=3.462e9
260 jgs 107 rho=1154.
261 jgs 110 tend=60
262 lkettle 581 h=(1./5.)*sqrt(rho/(lam+2*mu))*(width/ne)
263     print "time step size = ",h
264     U0=0.01 # amplitude of point source
265 jgs 110
266 lkettle 581 mydomain=Brick(ne,ne,10,l0=width,l1=width,l2=10.*width/32.)
267     wavePropagation(mydomain,h,tend,lam,mu,rho,U0)
268 jgs 107 \end{python}
269 lkettle 581 The script is available as \file{wave.py} in the \ExampleDirectory \index{scripts!\file{wave.py}}. Before running the code make sure you have created a directory called \code{data} in the current
270     working directory.
271     To visualize the results from the data directory:
272     \begin{verbatim} mayavi -d usoln.1.vtu -m SurfaceMap &\end{verbatim} You can rotate this figure by clicking on it with the mouse and moving it around.
273     Again use \code{Configure Data} to move backwards and forward in time, and
274     also to choose the results (acceleration, displacement or $u\hackscore x$) by using \code{Select Scalar}. Figure \ref{WAVE FIG 2} shows the results for the displacement at various time steps.
275 jgs 107
276 jfenwick 2104 \begin{figure}[t!]
277 gross 599 \centerline{\includegraphics[width=3.in,angle=270]{figures/WavePC.eps}}
278 gross 593 \caption{Amplitude at Point source}
279     \label{WAVE FIG 1}
280     \end{figure}
281 jgs 107
282 jfenwick 2104 \begin{figure}[t]
283 lkettle 581 \begin{center}
284 gross 599 \includegraphics[width=2in,angle=270]{figures/Wavet0.eps}
285     \includegraphics[width=2in,angle=270]{figures/Wavet1.eps}
286     \includegraphics[width=2in,angle=270]{figures/Wavet3.eps}
287     \includegraphics[width=2in,angle=270]{figures/Wavet10.eps}
288     \includegraphics[width=2in,angle=270]{figures/Wavet30.eps}
289     \includegraphics[width=2in,angle=270]{figures/Wavet288.eps}
290 lkettle 581 \end{center}
291     \caption{Selected time steps ($n = 0,1,30,100,300$ and 2880) of a wave propagation over a $10\mbox{km} \times 10\mbox{km} \times 3.125\mbox{km}$ block
292     from a point source initially at $(5\mbox{km},5\mbox{km},0)$
293     with time step size $h=0.02083$. Color represents the displacement.
294 gross 707 Here the view is oriented onto the bottom face.
295     \label{WAVE FIG 2}}
296 jgs 110 \end{figure}

Properties

Name Value
svn:eol-style native
svn:keywords Author Date Id Revision

  ViewVC Help
Powered by ViewVC 1.1.26