# Diff of /trunk/doc/user/wave.tex

revision 1315 by gross, Wed Apr 12 23:58:02 2006 UTC revision 1316 by ksteube, Tue Sep 25 03:18:30 2007 UTC
# Line 1  Line 1
1    %
2  % $Id$  % $Id$
3  %  %
5  %               \url{http://www.access.edu.au  %
6  %         Primary Business: Queensland, Australia.  %           Copyright 2003-2007 by ACceSS MNRF
7  %   Licensed under the Open Software License version 3.0  %       Copyright 2007 by University of Queensland
9    %                http://esscc.uq.edu.au
10    %        Primary Business: Queensland, Australia
13  %  %
14    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
15    %
16
17  \section{3D Wave Propagation}  \section{3D Wave Propagation}
18  \label{WAVE CHAP}  \label{WAVE CHAP}
19
20  We want to calculate the displacement field $u\hackscore{i}$ for any time $t>0$ by solving the wave equation  In this next example we want to calculate the displacement field $u\hackscore{i}$ for any time $t>0$ by solving the wave equation
21  \index{wave equation}:  \index{wave equation}:
22  \begin{eqnarray}\label{WAVE general problem}  \begin{eqnarray}\label{WAVE general problem}
23  \rho u\hackscore{i,tt} - \sigma\hackscore{ij,j}=0  \rho u\hackscore{i,tt} - \sigma\hackscore{ij,j}=0
# Line 40  Here we are modelling a point source at Line 48  Here we are modelling a point source at
48   set the initial displacement to be $U\hackscore 0$ in a sphere of small radius around the point   set the initial displacement to be $U\hackscore 0$ in a sphere of small radius around the point
49   $x\hackscore C$.     $x\hackscore C$.
50
51  We use an explicit time-integration scheme to calculate the displacement field $u$ at  We use an explicit time integration scheme to calculate the displacement field $u$ at
52  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$  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$
53  which is defined by  which is defined by
54  \begin{eqnarray} \label{WAVE dyn 2}  \begin{eqnarray} \label{WAVE dyn 2}
# Line 68  Now we have converted our problem for di Line 76  Now we have converted our problem for di
76  acceleration, $a^(n)$, which now depends  acceleration, $a^(n)$, which now depends
77  on the solution at the previous two time steps, $u^{(n-1)}$  and $u^{(n-2)}$.  on the solution at the previous two time steps, $u^{(n-1)}$  and $u^{(n-2)}$.
78
79  In each time step we have to solve this problem to get the acceleration $a^{(n)}$ and we will  In each time step we have to solve this problem to get the acceleration $a^{(n)}$, and we will
80  use the \LinearPDE class of the \linearPDEs to do so. The general form of the PDE defined through  use the \LinearPDE class of the \linearPDEs to do so. The general form of the PDE defined through
81  the \LinearPDE class is discussed in \Sec{SEC LinearPDE}. The form which is relevant here are    the \LinearPDE class is discussed in \Sec{SEC LinearPDE}. The form which is relevant here is
82  \label{WAVE dyn 100}  \label{WAVE dyn 100}
83  D\hackscore{ij} a^{(n)}\hackscore{j} = - X\hackscore{ij,j}\; .  D\hackscore{ij} a^{(n)}\hackscore{j} = - X\hackscore{ij,j}\; .
84
85  Implicitly, the natural boundary condition  The natural boundary condition
86  \label{WAVE dyn 101}  \label{WAVE dyn 101}
87  n\hackscore{j}X\hackscore{ij} =0  n\hackscore{j}X\hackscore{ij} =0
88
89  is assumed.  is used.
90
91  With $u=a^{(n)}$ we can identify the values to be assigned to $D$, $X$:  With $u=a^{(n)}$ we can identify the values to be assigned to $D$ and $X$:
92  \label{WAVE  dyn 6}  \label{WAVE  dyn 6}
93  \begin{array}{ll}  \begin{array}{ll}
94  D\hackscore{ij}=\rho \delta\hackscore{ij}&  D\hackscore{ij}=\rho \delta\hackscore{ij}&
# Line 95  X\hackscore{ij}=-\sigma^{(n-1)}\hackscor Line 103  X\hackscore{ij}=-\sigma^{(n-1)}\hackscor
103  The following script defines a the function \function{wavePropagation} which  The following script defines a the function \function{wavePropagation} which
104  implements the Verlet scheme to solve our wave propagation problem.  implements the Verlet scheme to solve our wave propagation problem.
105  The argument \var{domain} which is a \Domain class object  The argument \var{domain} which is a \Domain class object
106  defines the domain of the problem. \var{h} and \var{tend} are the step size  defines the domain of the problem. \var{h} and \var{tend} are the time step size
107  and the time mark after which the simulation will be terminated respectively. \var{lam}, \var{mu} and  and the end time of the simulation. \var{lam}, \var{mu} and
108  \var{rho} are material properties.  \var{rho} are material properties.
109  \begin{python}  \begin{python}
110  from esys.linearPDEs import LinearPDE  from esys.linearPDEs import LinearPDE
# Line 164  def wavePropagation(domain,h,tend,lam,mu Line 172  def wavePropagation(domain,h,tend,lam,mu
172     u_pc_data.close()     u_pc_data.close()
173  \end{python}  \end{python}
174  Notice that  Notice that
175  all coefficients of the PDE which are independent from time $t$ are set outside the \code{while}  all coefficients of the PDE which are independent of time $t$ are set outside the \code{while}
176  loop. This allows the \LinearPDE class to reuse information as much as possible  loop. This is very efficient since it allows the \LinearPDE class to reuse information as much as possible
177  when iterating over time.    when iterating over time.
178
179  We have seen most of the functions in previous examples but there are some new functions here:  there are a few new \escript functions in this example:
180  \code{grad(u)} returns the gradient $u\hackscore{i,j}$ of $u$ (in fact \var{grad(g)[i,j]}=$=u\hackscore{i,j}$).  \code{grad(u)} returns the gradient $u\hackscore{i,j}$ of $u$ (in fact \var{grad(g)[i,j]}=$=u\hackscore{i,j}$).
181  It is pointed out that in general there are restrictions on the argument of the \function{grad} function, for instance  There are restrictions on the argument of the \function{grad} function, for instance
183  \code{trace(g)} returns the sum of the main diagonal elements \var{g[k,k]} of \var{g}  \code{trace(g)} returns the sum of the main diagonal elements \var{g[k,k]} of \var{g}
184  and \code{transpose(g)} returns the transposed of the matrix \var{g} (ie.  and \code{transpose(g)} returns the matrix transpose of \var{g} (ie. $\var{transpose(g)[i,j]}=\var{g[j,i]}$).
$\var{transpose(g)[i,j]}=\var{g[j,i]}$).
185
186  We initialize the values of \code{u} and \code{u_last} to be $U\hackscore 0$ for a small  We initialize the values of \code{u} and \code{u_last} to be $U\hackscore 0$ for a small
187  sphere of radius, \code{src_radius} around the point source located at $x\hackscore C$.  sphere of radius, \code{src_radius} around the point source located at $x\hackscore C$.
# Line 192  plotted easily using any plotting packag Line 199  plotted easily using any plotting packag
199  %\end{python}  %\end{python}
200  %declares \var{u} as a vector-valued function of its coordinates  %declares \var{u} as a vector-valued function of its coordinates
201  %in the domain (i.e. with two or three components in the case of  %in the domain (i.e. with two or three components in the case of
202  %a two or three dimensional domain, repectively). The first argument defines the value of the function which is equal  %a two or three dimensional domain, respectively). The first argument defines the value of the function which is equal
203  %to $0$ for all vectot components. The second argument \var{ContinuousFunction(domain)}  %to $0$ for all vectot components. The second argument \var{ContinuousFunction(domain)}
204  %specifies the smoothness of the function. In our case we want the displacement field to be  %specifies the smoothness of the function. In our case we want the displacement field to be
205  %a continuous function over the \Domain \var{domain}. On other option would be  %a continuous function over the \Domain \var{domain}. On other option would be
206  %\var{Function(domain)} in which case continuity is not garanteed. In fact the  %\var{Function(domain)} in which case continuity is not guaranteed. In fact the
207  %argument of the \function{grad} function has to be of the typ \var{ContinuousFunction(domain)} while  %argument of the \function{grad} function has to be of the typ \var{ContinuousFunction(domain)} while
208  %returned function is of the var{Function(domain)} type. For more details see \Sec{SEC ESCRIPT DATA}.  %returned function is of the var{Function(domain)} type. For more details see \Sec{SEC ESCRIPT DATA}.
209
210  The statement \code{myPDE.setLumpingOn()}  The statement \code{myPDE.setLumpingOn()}
211  switches on the usage of an aggressive approximation of the PDE operator, in this case  switches on the use of an aggressive approximation of the PDE operator as a diagonal matrix
212  formed by the coefficient \var{D} (actually the discrete  formed from the coefficient \var{D}.
213  version of the PDE operator is approximated by a diagonal matrix). The approximation allows, at the costs of  The approximation allows, at the cost of
215  solving of the PDE when the solution is requested. In the general case, the presence of \var{A}, \var{B} or \var  solution of the PDE. When using \code{setLumpingOn} the presence of \var{A}, \var{B} or \var
216  {C} \code{setLumpingOn} will produce wrong results.  {C} will produce wrong results.
217
218
219  One of the big advantages of the Verlet scheme is the fact that the problem to be solved  One of the big advantages of the Verlet scheme is the fact that the problem to be solved
220  in each time step is very simple and does not involve any spatial derivatives (which actually allows us to use lumping).  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).
221  The problem becomes so simple because we use the stress from the last time step rather then the stress which is  The problem becomes so simple because we use the stress from the last time step rather then the stress which is
222  actually present at the current time step. Schemes using this approach are called an explicit time integration  actually present at the current time step. Schemes using this approach are called an explicit time integration
223  scheme \index{explicit scheme} \index{time integration!explicit}. The  schemes \index{explicit scheme} \index{time integration!explicit}. The
224  backward Euler scheme we have used in \Chap{DIFFUSION CHAP} is  backward Euler scheme we have used in \Chap{DIFFUSION CHAP} is
225  an example of an implicit schemes  an example of an implicit scheme
226  \index{implicit scheme} \index{time integration!implicit}. In this case one uses the actual status of  \index{implicit scheme} \index{time integration!implicit}. In this case one uses the actual status of
227  each variable at a particular time rather then values from previous time steps. This will lead to a problem which is  each variable at a particular time rather then values from previous time steps. This will lead to a problem which is
228  more expensive to solve, in particular for non-linear problems.  more expensive to solve, in particular for non-linear problems.

Legend:
 Removed from v.1315 changed lines Added in v.1316