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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 581 - (hide annotations)
Wed Mar 8 05:48:51 2006 UTC (13 years, 6 months ago) by lkettle
File MIME type: application/x-tex
File size: 15046 byte(s)
changes to section 2.3 on 3D wave propagation in online reference guide and added some figures

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

Properties

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

  ViewVC Help
Powered by ViewVC 1.1.26