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

revision 107 by jgs, Thu Jan 27 06:21:48 2005 UTC revision 110 by jgs, Mon Feb 14 04:14:42 2005 UTC
# Line 71  a^{(n)}\hackscore{2}(x)= s\hackscore{,tt Line 71  a^{(n)}\hackscore{2}(x)= s\hackscore{,tt
71  \end{eqnarray}  \end{eqnarray}
72  derived from \eqn{WAVE impact}.  derived from \eqn{WAVE impact}.
73
74
75  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
76  use the \LinearPDE class of the \linearPDEsPack to do so. The general form of the PDE defined through  use the \LinearPDE class of the \linearPDEsPack to do so. The general form of the PDE defined through
77  the \LinearPDE class is discussed in \Sec{SEC LinearPDE}. The form which are relevant here are    the \LinearPDE class is discussed in \Sec{SEC LinearPDE}. The form which are relevant here are
78  \label{WAVE dyn 100}  \label{WAVE dyn 100}
79  D\hacksco  D\hackscore{ij} u\hackscore{j} = - X\hackscore{ij,j}\; .
re{ij} u\hackscore{j} = - X\hackscore{ij,j}\; .
80
81  Impliciltly, the natural boundary condition  Implicitly, the natural boundary condition
82  \label{WAVE dyn 101}  \label{WAVE dyn 101}
83  n\hackscore{j}X\hackscore{ij} =0  n\hackscore{j}X\hackscore{ij} =0
84
85  is assumed. The \LinearPDE object allows defining constriants of the form  is assumed. The \LinearPDE object allows defining constraints of the form
86  \label{WAVE dyn 101}  \label{WAVE dyn 101}
87  u\hackscore{i}(x)=r\hackscore{i}(x) \mbox{ for } x \mbox{ with } q\hackscore{i}>0  u\hackscore{i}(x)=r\hackscore{i}(x) \mbox{ for } x \mbox{ with } q\hackscore{i}>0
88
# Line 104  implements the Verlet scheme to solve ou Line 104  implements the Verlet scheme to solve ou
104  The argument \var{domain} which is a \Domain class object  The argument \var{domain} which is a \Domain class object
105  defines the domain of the problem. \var{h} and \var{tend} is the step size  defines the domain of the problem. \var{h} and \var{tend} is the step size
106  and the time mark after which the simulation will be terminated. \var{lam}, \var{mu} and  and the time mark after which the simulation will be terminated. \var{lam}, \var{mu} and
107  \var{rho} are material properties. \var{s_tt} is a fucntion which returns $s\hackscore{,tt}$ at a given time:  \var{rho} are material properties. \var{s_tt} is a function which returns $s\hackscore{,tt}$ at a given time:
108  \begin{python}  \begin{python}
109  from esys.linearPDEs import LinearPDE  from esys.linearPDEs import LinearPDE
110  from numarray import identity  from numarray import identity
111  from esys.escript import *  from esys.escript import *
112  def wavePropagation(domain,h,tend,lam,mu,rho,s_tt):  def wavePropagation(domain,h,tend,lam,mu,rho,s_tt):
113     x=domain.getX()     x=domain.getX()
ndim=domain.getDim()
114     # ... open new PDE ...     # ... open new PDE ...
115     mypde=LinearPDE(domain)     mypde=LinearPDE(domain)
116     mypde.setLumpingOn()     mypde.setLumpingOn()
117     kronecker=identity(ndim)     kronecker=identity(mypde.getDim())
118     mypde.setValue(D=kronecker*rho, \     mypde.setValue(D=kronecker*rho, \
119                    q=x[0].whereZero()*kronecker[ndim-1,:])                    q=x[0].whereZero()*kronecker[1,:])
120     # ... set initial values ....     # ... set initial values ....
121       n=0
122     u=Vector(0,ContinuousFunction(domain))     u=Vector(0,ContinuousFunction(domain))
123     u_last=Vector(0,ContinuousFunction(domain))     u_last=Vector(0,ContinuousFunction(domain))
124     t=0     t=0
125     while t<tend:     while t<tend:
126       # ... get current stress ....       # ... get current stress ....
128       stress=trace(g)*lam*kronecker+mu*(g+transpose(g))       stress=lam*trace(g)*kronecker+mu*(g+transpose(g))
129       # ... get new acceleration ....           # ... get new acceleration ....
130       mypde.setValue(X=-stress,r=s_tt(t+h)*kronecker[ndim-1,:])       mypde.setValue(X=-stress,r=s_tt(t+h)*kronecker[1,:])
131       a=mypde.getSolution()       a=mypde.getSolution()
132       # ... get new displacement ...       # ... get new displacement ...
133       u_new=2*u-u_last+h**2*a       u_new=2*u-u_last+h**2*a
134       # ... shift displacements ....       # ... shift displacements ....
135       u_last=u       u_last=u
136       u=u_new       u=u_new
# ... next time step ...
137       t+=h       t+=h
138         n+=1
139         print n,"-th time step t ",t
140         print "a=",inf(a),sup(a)
141         print "u=",inf(u),sup(u)
142         # ... save current acceleration in units of gravity
143         if n%10==0 : (length(a)/9.81).saveDX("/tmp/res/u.%i.dx"%(n/10))
144  \end{python}  \end{python}
145  We have seen most of the functions in previous examples. However,    Notice that
146  there are some new functions here:  all coefficients of the PDE which are independent from time $t$ are set outside the \code{while}
147    loop. This allows the \LinearPDE class to ruse information as much as possible
148    when iterating over time.
149
150    We have seen most of the functions in previous examples but there are some new functions here:
151  \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}$).
152  It is pointed out that in general there are restrictions on the argument of the \function{grad} function, for instance  It is pointed out that in general there are restrictions on the argument of the \function{grad} function, for instance
# Line 148  $\var{transpose(g)[i,j]}=\var{g[j,i]}$). Line 157  $\var{transpose(g)[i,j]}=\var{g[j,i]}$).
157
158  The initialization of \var{u} and \var{u_last} needs some explanation. The statement  The initialization of \var{u} and \var{u_last} needs some explanation. The statement
159  \begin{python}  \begin{python}
160  u=Vector(0,ContinuousFunction(domain))  u=Vector(0.,ContinuousFunction(domain))
161  \end{python}  \end{python}
162  creates a new object which is of class \Data, see \Sec{SEC ESCRIPT DATA}. The object represnts a vector-valued  declares \var{u} as a vector-valued function of its coordinates
163  function of the location coordinates. Vector-value means that it has two or three components on a two or three dimensional domain. respectively. The second argument The first argument (here $=0$) is the initial value  in the domain (i.e. with two or three components in the case of
164  assigned to any spatial point.  a two or three dimensional domain, repectively). The first argument defines the value of the function which is equal
165    to $0$ for all vectot components. The second argument \var{ContinuousFunction(domain)}
166    specifies the smoothness of the function. In our case we want the displacement field to be
167    a continuous function over the \Domain \var{domain}. On other option would be
168    \var{Function(domain)} in which case continuity is not garanteed. In fact the
169    argument of the \function{grad} function has to be of the typ \var{ContinuousFunction(domain)} while
170    returned function is of the var{Function(domain)} type. For more details see \Sec{SEC ESCRIPT DATA}.
171
172    The statement \code{myPDE.setLumpingOn()}
All coefficients of the PDE which are independent from time are set outside the \code{while}
loop. This allows the \LinearPDE class to ruse information as much as possible
when iterating over time.  The statement \code{myPDE.setLumpingOn()}
173  switches on the usage of an aggressive approximation of the PDE operator, in this case  switches on the usage of an aggressive approximation of the PDE operator, in this case
174  formed by the coefficient \var{D} (actually the discrete  formed by the coefficient \var{D} (actually the discrete
175  version of the PDE operator is approximated by a diagonal matrix). The approximation allows, at the costs of  version of the PDE operator is approximated by a diagonal matrix). The approximation allows, at the costs of
177  solving of the PDE when the solution is requested. In the general case of the presents of \var{A}, \var{B} or \var  solving of the PDE when the solution is requested. In the general case of the presents of \var{A}, \var{B} or \var
178  {C} \code{setLumpingOn} will produce wrong results.  {C} \code{setLumpingOn} will produce wrong results.
179  =========
180  It is pointed out that the value of the constraint \var{r} for the acceleration \var{a} at time $t=0$ is consistent  It is pointed out that the value of the constraint \var{r} for the acceleration \var{a} at time $t=0$ is consistent
181  with the initial value of the acceleration which is identical zero. Otherwise, a discontinuity of \var{a} in time  with the initial value of the acceleration which is identical zero. Otherwise, a discontinuity of \var{a} in time
182  direction would lead to heavy oscillations in the displacement. It is also pointed out, that at time $t=0$ the  direction would lead to heavy oscillations in the displacement. It is also pointed out, that at time $t=0$ the
# Line 173  constraint defined by \eqn{WAVE impact} Line 184  constraint defined by \eqn{WAVE impact}
184  the initial conditions in \eqn{WAVE initial}.  the initial conditions in \eqn{WAVE initial}.
185
186
187  ================  One of the big advantages of the Verlet scheme is the fact that the problem to be solved
188  The advantage of the Verlat scheme is that the PDE we have to solve in each time step is very  in each time step is very simple and does not involve any spatial derivatives (which actually allows us to use lumping).
189  simple as we are using the stress based on the displacements of last time step. One could, similar to the  The problem becomes so simple because we use the stress from the last time step rather then the stress which
190  backward Euler scheme we have used in \Chap{DIFFUSION CHAP}, use the stress calculated with current  actually is present at the current time step. Schemes using this approach are called an explicit time integration
191  displacement which would lead to a more complex PDE involving a coefficient \var{A}. These so-called implicit schemes  scheme \index{explicit scheme} \index{time integration!explicit}. The
192  \index{implicit scheme} \index{time integration!implicit} is more expensive in each time then the Verlat scheme which  backward Euler scheme we have used in \Chap{DIFFUSION CHAP} is
193  is a so called explicit scheme \index{explicit scheme} \index{time integration!explicit}.      an example of an implicit schemes
194    \index{implicit scheme} \index{time integration!implicit}. In this case one uses the actual status of
195    each variable at a particular time rather then values from previous time steps. This will lead to problem which is
196    more expensive to solve, in particular for non-linear problems.
197    Although
198    explicit time integration schemes are cheep to finalize a single time step, they need a significant smaller time
199    step then implicit schemes and can suffer from stability problems. Therefor they need a
200    very careful selection of the time step size $h$.
201
202  Explicit time integration schemes need a very careful selection of the time step size $h$ otherwise  An easy, heuristic way of choosing an appropriate
if $h$ is too big the scheme can be unstable. An easy, heuristic way of choosing an appropriate
203  time step size is the Courant condition \index{Courant condition} \index{explicit scheme!Courant condition}  time step size is the Courant condition \index{Courant condition} \index{explicit scheme!Courant condition}
204  which says that within a time step a information should not travel further than a cell used in the  which says that within a time step a information should not travel further than a cell used in the
205  discretization scheme. In the case of the wave equation the velocity of a (p-) wave is given as  discretization scheme. In the case of the wave equation the velocity of a (p-) wave is given as
# Line 194  where $\Delta x$ is the cell diameter. T Line 211  where $\Delta x$ is the cell diameter. T
211  the formula.  the formula.
212
213  The following script uses the \code{wavePropagtion} function to solve the  The following script uses the \code{wavePropagtion} function to solve the
214  wave equation over a block in crust of the earth. The  wave equation over a block in crust of the earth. The depth of the block is
215  depth is $l\hackscore{0}=20km$ and the length is $l\hackscore{1}=l\hackscore{2}=200km$.  $10km$ and the width in each direction is $100km$. The depth of the block is alligned
216  The time of the impact is about $2sec$ and maximum displacement is $u^{max}=15m$ which corresponds to an heavy  with the $x\hackscore{0}$-direction. \var{ne} gives the number of elements in $x\hackscore{0}$-direction. The number
217  earthquake. We are using a $6 \times 60 \times 60$ mesh:  of elements in $x\hackscore{1}$- and $x\hackscore{2}$-direction is chosen such that the elements
218    become cubic. The impact is $2m$ acting within a time of about $10sec$:
219  \begin{python}  \begin{python}
220  ne=6  ne=10           # number of cells in x_0-direction
221  lame_lambda=3.462e9  depth=10000.   # length in x_0-direction
222  lame_mu=3.462e9  width=100000.  # length in x_1 and x_2 direction
223    lam=3.462e9
224    mu=3.462e9
225  rho=1154.  rho=1154.
226  tau=2.  tau=10.
227  umax=15.  umax=2.
228  xc=[0,1000,1000]  tend=60
229  r=1.  h=1./5.*sqrt(rho/(lam+2*mu))*(depth/ne)
230  tend=10.  def s_tt(t): return umax/tau**2*(6*t/tau-9*(t/tau)**4)*exp(-(t/tau)**3)
231  dt=1./5.*sqrt(rho/(lame_lambda+2*lame_mu))(20000./ne)
232  print "step size = ",dt  print "time step size = ",h
233  mydomain=finely.Brick(ne,10*ne,10*ne,l0=20000,l1=200000,l2=200000)  mydomain=Brick(ne,int(width/depth)*ne,int(width/depth)*ne,l0=depth,l1=width,l2=width)
234  wavePropagation(mydomain,dt,tend,lame_lambda,lame_mu,rho,xc,r,tau,umax)  wavePropagation(mydomain,h,tend,lam,mu,rho,s_tt)
235  \end{python}  \end{python}
236  The script is available as \file{wave.py} in the \ExampleDirectory \index{scripts!\file{wave.py}}.  The script is available as \file{wave.py} in the \ExampleDirectory \index{scripts!\file{wave.py}}.
237
# Line 222  The script is available as \file{wave.py Line 242  The script is available as \file{wave.py
242  % \label{WAVE FIG 1}  % \label{WAVE FIG 1}
243  % \end{figure}  % \end{figure}
244
245  % \begin{figure}  \begin{figure}
246  % \centerline{\includegraphics[width=\figwidth]{DiffusionRes1}}  % \centerline{\includegraphics[width=\figwidth]{DiffusionRes1}}
247  % \centerline{\includegraphics[width=\figwidth]{DiffusionRes16}}  % \centerline{\includegraphics[width=\figwidth]{DiffusionRes16}}
248  % \centerline{\includegraphics[width=\figwidth]{DiffusionRes32}}  % \centerline{\includegraphics[width=\figwidth]{DiffusionRes32}}
249  %\centerline{\includegraphics[width=\figwidth]{DiffusionRes48}}  %\centerline{\includegraphics[width=\figwidth]{DiffusionRes48}}
250  % \caption{Results of the Temperture Diffusion Problem for Time Steps $1$ $16$, $32$ and $48$.}  \caption{Selected time steps of a wave propagation over a $10km \times 100km \times 100km$
251  % \label{WAVE FIG 2}  with time step size $h=0.0666$. Color represents the acceleration.}
252  %\end{figure}  \label{WAVE FIG 2}
253    \end{figure}
254
255

Legend:
 Removed from v.107 changed lines Added in v.110