# Diff of /branches/stage3.0/doc/user/wave.tex

revision 2583 by jfenwick, Fri Jul 31 05:39:36 2009 UTC revision 2584 by jfenwick, Wed Aug 5 02:44:51 2009 UTC
# Line 48  for all time $t>0$. Line 48  for all time $t>0$.
48  \end{figure}  \end{figure}
49
50  Here we are modelling a point source at the point $x\hackscore C$ in the $x\hackscore{0}$-direction  Here we are modelling a point source at the point $x\hackscore C$ in the $x\hackscore{0}$-direction
51  which is a negative puls of amplitude $U\hackscore{0}$ followd by the same  which is a negative pulse of amplitude $U\hackscore{0}$ followed by the same
52  positive puls. In mathematical terms we use  positive pulse. In mathematical terms we use
53  \begin{eqnarray} \label{WAVE source}  \begin{eqnarray} \label{WAVE source}
54  u\hackscore{0}(x\hackscore C,t)= U\hackscore{0} \; \sqrt{2}  \frac{(t-t\hackscore{0})}{\alpha} e^{\frac{1}{2}-\frac{(t-t\hackscore{0})^2}{\alpha^2}} \  u\hackscore{0}(x\hackscore C,t)= U\hackscore{0} \; \sqrt{2}  \frac{(t-t\hackscore{0})}{\alpha} e^{\frac{1}{2}-\frac{(t-t\hackscore{0})^2}{\alpha^2}} \
55  \end{eqnarray}  \end{eqnarray}
56  for all $t\ge0$ where $\alpha$ is the width of the puls and $t\hackscore{0}$ is the time when  for all $t\ge0$ where $\alpha$ is the width of the pulse and $t\hackscore{0}$ is the time when
57  the puls changes from negative to positive. In the simulations we will choose $\alpha=0.3$ and $t\hackscore{0}=2$, see Figure~\ref{WAVE FIG 1.2}  the pulse changes from negative to positive. In the simulations we will choose $\alpha=0.3$ and $t\hackscore{0}=2$, see Figure~\ref{WAVE FIG 1.2}
58  and will apply the source as a constraint in a in a sphere of small radius around the point  and will apply the source as a constraint in a in a sphere of small radius around the point
59  $x\hackscore C$.    $x\hackscore C$.
60
# Line 129  and the end time of the simulation. \var Line 129  and the end time of the simulation. \var
129  \var{rho} are material properties.  \var{rho} are material properties.
130  \begin{python}  \begin{python}
132       # lists to collect displacement at point source which is returned to the caller
133       ts, u_pc0,u_pc1,u_pc2=[], [], [], []
134
135     x=domain.getX()     x=domain.getX()
136     # ... open new PDE ...     # ... open new PDE ...
137     mypde=LinearPDE(domain)     mypde=LinearPDE(domain)
138     mypde.getSolverOptions().setSolverMethod(mypde.getSolverOptions().LUMPING)     mypde.getSolverOptions().setSolverMethod(mypde.getSolverOptions().LUMPING)
139     kronecker=identity(mypde.getDim())     kronecker=identity(mypde.getDim())

140     dunit=numpy.array([1.,0.,0.]) # defines direction of point source     dunit=numpy.array([1.,0.,0.]) # defines direction of point source

142     # ... set initial values ....     # ... set initial values ....
143     n=0     n=0
# Line 146  def wavePropagation(domain,h,tend,lam,mu Line 145  def wavePropagation(domain,h,tend,lam,mu
145     u=Vector(0.,Solution(domain))     u=Vector(0.,Solution(domain))
146     u_last=Vector(0.,Solution(domain))     u_last=Vector(0.,Solution(domain))
147     t=0     t=0

148     # define the location of the point source     # define the location of the point source
149     L=Locator(domain,xc)     L=Locator(domain,xc)
150     # find potential at point source     # find potential at point source
# Line 154  def wavePropagation(domain,h,tend,lam,mu Line 152  def wavePropagation(domain,h,tend,lam,mu
152     print "u at point charge=",u_pc     print "u at point charge=",u_pc
153     # open file to save displacement at point source     # open file to save displacement at point source
154     u_pc_data=FileWriter('./data/U_pc.out')     u_pc_data=FileWriter('./data/U_pc.out')
155     u_pc_data.write("%f %f %f %f\n"%(t,u_pc[0],u_pc[1],u_pc[2]))     ts.append(t); u_pc0.append(u_pc[0]), u_pc1.append(u_pc[1]), u_pc2.append(u_pc[2])
156
157     while t<tend:     while t<tend:
158       t+=h       t+=h
# Line 162  def wavePropagation(domain,h,tend,lam,mu Line 160  def wavePropagation(domain,h,tend,lam,mu
161       stress=lam*trace(g)*kronecker+mu*(g+transpose(g))       stress=lam*trace(g)*kronecker+mu*(g+transpose(g))
162       # ... get new acceleration ....       # ... get new acceleration ....
163       amplitude=U0*(4*(t-t0)**3/alpha**3-6*(t-t0)/alpha)*sqrt(2.)/alpha**2*exp(1./2.-(t-t0)**2/alpha**2)       amplitude=U0*(4*(t-t0)**3/alpha**3-6*(t-t0)/alpha)*sqrt(2.)/alpha**2 \
164                                                 *exp(1./2.-(t-t0)**2/alpha**2)
165       mypde.setValue(X=-stress, r=dunit*amplitude)       mypde.setValue(X=-stress, r=dunit*amplitude)
166       a=mypde.getSolution()       a=mypde.getSolution()
167       # ... get new displacement ...       # ... get new displacement ...
# Line 175  def wavePropagation(domain,h,tend,lam,mu Line 174  def wavePropagation(domain,h,tend,lam,mu
174       u_pc=L.getValue(u)       u_pc=L.getValue(u)
175       print "u at point charge=",u_pc       print "u at point charge=",u_pc
176       # save displacements at point source to file for t > 0       # save displacements at point source to file for t > 0
177       u_pc_data.write("%f %f %f %f\n"%(t,u_pc[0],u_pc[1],u_pc[2]))       ts.append(t); u_pc0.append(u_pc[0]), u_pc1.append(u_pc[1]), \
178                                                        u_pc2.append(u_pc[2])
179
180       # ... save current acceleration in units of gravity and displacements       # ... save current acceleration in units of gravity and displacements
181       if n==1 or n%10==0: saveVTK("./data/usoln.%i.vtu"%(n/10),acceleration=length(a)/9.81,       if n==1 or n%10==0: saveVTK("./data/usoln.%i.vtu"%(n/10), \
182       displacement = length(u), tensor = stress, Ux = u[0] )                                   acceleration=length(a)/9.81,
183                                     displacement = length(u), \
184                     tensor = stress, Ux = u[0] )
185
186     u_pc_data.close()     return ts, u_pc0, u_pc1, u_pc2
187  \end{python}  \end{python}
188  Notice that  Notice that
189  all coefficients of the PDE which are independent of 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}
# Line 214  L=Locator(domain,xc) Line 216  L=Locator(domain,xc)
216  u=...  u=...
217  u_pc=L.getValue(u)  u_pc=L.getValue(u)
218  \end{python}  \end{python}
219  The return value \code{u_pc} is the value of \code{u} at the location \code{xc}\footnote{In fact the finite element node which is closest to the given position. The usage of  \class{Locator} is MPI save.}.  The return value \code{u_pc} is the value of \code{u} at the location \code{xc}\footnote{In fact the finite element node which is closest to the given position. The usage of  \class{Locator} is MPI save.}. The values
220    are collected in the lists \var{u_pc0}, \var{u_pc1} and \var{u_pc2} together with the
221    corresponding time marker in \var{ts}. The values are handed back to the caller. Later we will show to ways to
222  The output \code{U_pc.out} and \code{vtu} files are saved in a directory called \code{data}.  access these data.
The \code{U_pc.out} stores 4 columns of data: $t,u\hackscore x,u\hackscore y,u\hackscore z$
respectively, where $u\hackscore x,u\hackscore y,u\hackscore z$ are the $x\hackscore 0,x\hackscore 1,x\hackscore 2$ components of
the displacement vector $u$ at the point source. These can be
plotted easily using any plotting package. In gnuplot the command:
\begin{verbatim}
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,
'U_pc.out' u 1:4 title 'U_z' w l lw 2
\end{verbatim}
will reproduce Figure~\ref{WAVE FIG 1} (As expected this is identical to the input signal shown in Figure~\ref{WAVE FIG 1.2})2. It is pointed out that we are not using the
standart \PYTHON \function{open} to create the file \code{U_pc.out} as it is not
when running \escript under MPI, see chapter~\ref{EXECUTION} for more details.
223
224  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
225  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).  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).
# Line 256  h= \frac{1}{5} \sqrt{\frac{\rho}{\lambda Line 247  h= \frac{1}{5} \sqrt{\frac{\rho}{\lambda
247  where $\Delta x$ is the cell diameter. The factor $\frac{1}{5}$ is a safety factor considering the heuristics of  where $\Delta x$ is the cell diameter. The factor $\frac{1}{5}$ is a safety factor considering the heuristics of
248  the formula.  the formula.
249
250    \begin{figure}[t]
251    \begin{center}
252    \includegraphics[width=2in]{figures/Wave11}
253    \includegraphics[width=2in]{figures/Wave22}
254    \includegraphics[width=2in]{figures/Wave28}
255    \includegraphics[width=2in]{figures/Wave32}
256    \includegraphics[width=2in]{figures/Wave36}
257    \end{center}
258    \caption{Selected time steps ($n = 11, 22, 32, 36$) of a wave propagation over a $10\mbox{km} \times 10\mbox{km} \times 3.125\mbox{km}$ block
259    from a point source initially at $(5\mbox{km},5\mbox{km},0)$
260    with time step size $h=0.02083$. Color represents the displacement.
261    Here the view is oriented onto the bottom face.
262    \label{WAVE FIG 2}}
263    \end{figure}
264
265  The following script uses the \code{wavePropagation} function to solve the  The following script uses the \code{wavePropagation} function to solve the
266  wave equation for a point source located at the bottom face of a block. The width of the block in  wave equation for a point source located at the bottom face of a block. The width of the block in
267  each direction on the bottom face is $10\mbox{km}$ ($x\hackscore 0$ and $x\hackscore 1$ directions ie. \code{l0} and \code{l1}).  each direction on the bottom face is $10\mbox{km}$ ($x\hackscore 0$ and $x\hackscore 1$ directions ie. \code{l0} and \code{l1}).
# Line 281  xc=[width/2.,width/2.,0.] Line 287  xc=[width/2.,width/2.,0.]
287  # define small radius around point xc  # define small radius around point xc

290  mydomain=Brick(ne,ne,10,l0=width,l1=width,l2=10.*width/32.)  mydomain=Brick(ne,ne,10,l0=width,l1=width,l2=10.*width/32.)
291  h=(1./5.)*inf(sqrt(rho/(lam+2*mu))*inf(domain.getSize())  h=(1./5.)*inf(sqrt(rho/(lam+2*mu))*inf(domain.getSize())
292  print "time step size = ",h  print "time step size = ",h
293  wavePropagation(mydomain,h,tend,lam,mu,rho,xc, src_radius, U0)  ts, u_pc0, u_pc1, u_pc2 =   \
295  \end{python}  \end{python}
296  The \function{domain.getSize()} return the local element size $\Delta x$. The  The \function{domain.getSize()} return the local element size $\Delta x$. The
297  \function{inf} makes sure that the Courant condition~\ref{WAVE dyn 66} olds everywhere in the domain.  \function{inf} makes sure that the Courant condition~\ref{WAVE dyn 66} olds everywhere in the domain.
298
299  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}}.
300  To visualize the results from the data directory:  To visualize the results from the data directory:
301  \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.  \begin{verbatim}
302    mayavi2 -d usoln.1.vtu -m SurfaceMap &
303    \end{verbatim}
304    You can rotate this figure by clicking on it with the mouse and moving it around.
305  Again use \code{Configure Data} to move backwards and forward in time, and  Again use \code{Configure Data} to move backwards and forward in time, and
306  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.  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.
307
# Line 302  also to choose the results (acceleration Line 311  also to choose the results (acceleration
311  \label{WAVE FIG 1}  \label{WAVE FIG 1}
312  \end{figure}  \end{figure}
313
314  \begin{figure}[t]  It remains to show some possibilities to inspect the collected data \var{u_pc0}, \var{u_pc1} and \var{u_pc2}.
315  \begin{center}  One way is to write the data to a file and then use an external package such as \gnuplot, excel or OpenOffice to read the data for further analysis. The following code shows one possible form to write the data to the
316  \includegraphics[width=2in]{figures/Wave11}  file \file{./data/U_pc.out}:
317  \includegraphics[width=2in]{figures/Wave22}  \begin{python}
318  \includegraphics[width=2in]{figures/Wave28}  u_pc_data=FileWriter('./data/U_pc.out')
319  \includegraphics[width=2in]{figures/Wave32}  for i in xrange(len(ts)) :
320  \includegraphics[width=2in]{figures/Wave36}      u_pc_data.write("%f %f %f %f\n"%(ts[i],u_pc0[i],u_pc1[i],u_pc2[i]))
321  \end{center}  u_pc_data.close()
322  \caption{Selected time steps ($n = 11, 22, 32, 36$) of a wave propagation over a $10\mbox{km} \times 10\mbox{km} \times 3.125\mbox{km}$ block  \end{python}
323  from a point source initially at $(5\mbox{km},5\mbox{km},0)$  The \code{U_pc.out} stores 4 columns of data: $t,u\hackscore x,u\hackscore y,u\hackscore z$
324  with time step size $h=0.02083$. Color represents the displacement.  respectively, where $u\hackscore x,u\hackscore y,u\hackscore z$ are the $x\hackscore 0,x\hackscore 1,x\hackscore 2$ components of the displacement vector $u$ at the point source. These can be
325  Here the view is oriented onto the bottom face.  plotted easily using any plotting package. In \gnuplot the command:
326  \label{WAVE FIG 2}}  \begin{verbatim}
327  \end{figure}   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,
328    'U_pc.out' u 1:4 title 'U_z' w l lw 2
329    \end{verbatim}
330    will reproduce Figure~\ref{WAVE FIG 1} (As expected this is identical to the input signal shown in Figure~\ref{WAVE FIG 1.2})2. It is pointed out that we are not using the
331    standard \PYTHON \function{open} to write to the file \code{U_pc.out} as it is not safe
332    when running \escript under MPI, see chapter~\ref{EXECUTION} for more details.
333
334    Alternatively, one can implement plotting the results at run time rather than in a post-processing step. This avoids
335    the generation of an intermediate data file. In {\it escript} the preferred way of creating 2D plots of
336    time dependent data is \MATPLOTLIB. The following script creates the plot and writes it into the
337    file \file{u_pc.png} in the PNG image format:
338    \begin{python}
339    import matplotlib.pyplot as plt
340    if getMPIRankWorld() == 0:
341        plt.title("Displacement at Point Source")
342        plt.plot(ts, u_pc0, '-', label="x_0", linewidth=1)
343        plt.plot(ts, u_pc1, '-', label="x_1", linewidth=1)
344        plt.plot(ts, u_pc2, '-', label="x_2", linewidth=1)
345        plt.xlabel('time')
346        plt.ylabel('displacement')
347        plt.legend()
348        plt.savefig('u_pc.png', format='png')
349    \end{python}
350    You can add the \function{plt.show()} to create a interactive browser window. Please not that
351    through the \code{getMPIRankWorld() == 0} statement the plot is generated on one processor only (in this case the rank 0 processor) when run under MPI.
352
353    Both options for processing the point source data are include in the example file \file{wave.py}. There other options available to process these data in particular through the \SCIPY package , eg Fourier transformations. It is beyond the scope of this users guide to document the usage of \SCIPY for time series analysis but is highly recommended that users  use \SCIPY to any further data processing.
354

Legend:
 Removed from v.2583 changed lines Added in v.2584