| 1 |
|
| 2 |
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% |
| 3 |
% |
| 4 |
% Copyright (c) 2003-2009 by University of Queensland |
| 5 |
% Earth Systems Science Computational Center (ESSCC) |
| 6 |
% http://www.uq.edu.au/esscc |
| 7 |
% |
| 8 |
% 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 |
% |
| 12 |
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% |
| 13 |
|
| 14 |
|
| 15 |
\section{3-D Wave Propagation} |
| 16 |
\label{WAVE CHAP} |
| 17 |
|
| 18 |
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 |
\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 |
for all time $t>0$. |
| 37 |
|
| 38 |
\begin{figure}[t!] |
| 39 |
\centerline{\includegraphics[angle=-90,width=4.in]{figures/waveu}} |
| 40 |
\caption{Input Displacement at Source Point ($\alpha=0.7$, $t\hackscore{0}=3$, $U\hackscore{0}=1$).} |
| 41 |
\label{WAVE FIG 1.2} |
| 42 |
\end{figure} |
| 43 |
|
| 44 |
\begin{figure}[t!] |
| 45 |
\centerline{\includegraphics[angle=-90,width=4.in]{figures/wavea}} |
| 46 |
\caption{Input Acceleration at Source Point ($\alpha=0.7$, $t\hackscore{0}=3$, $U\hackscore{0}=1$).} |
| 47 |
\label{WAVE FIG 1.1} |
| 48 |
\end{figure} |
| 49 |
|
| 50 |
Here we are modelling a point source at the point $x\hackscore C$ in the $x\hackscore{0}$-direction |
| 51 |
which is a negative pulse of amplitude $U\hackscore{0}$ followed by the same |
| 52 |
positive pulse. In mathematical terms we use |
| 53 |
\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}} \ |
| 55 |
\end{eqnarray} |
| 56 |
for all $t\ge0$ where $\alpha$ is the width of the pulse and $t\hackscore{0}$ is the time when |
| 57 |
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 |
| 59 |
$x\hackscore C$. |
| 60 |
|
| 61 |
We use an explicit time integration scheme to calculate the displacement field $u$ at |
| 62 |
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$ |
| 63 |
which is defined by |
| 64 |
\begin{eqnarray} \label{WAVE dyn 2} |
| 65 |
u^{(n)}=2u^{(n-1)}-u^{(n-2)} + h^2 a^{(n)} \\ |
| 66 |
\end{eqnarray} |
| 67 |
for all $n=2,3,\ldots$. It is designed to solve a system of equations of the form |
| 68 |
\begin{eqnarray} \label{WAVE dyn 2b} |
| 69 |
u\hackscore{,tt}=G(u) |
| 70 |
\end{eqnarray} |
| 71 |
where one sets $a^{(n)}=G(u^{(n-1)})$. |
| 72 |
|
| 73 |
In our case $a^{(n)}$ is given by |
| 74 |
\begin{eqnarray}\label{WAVE dyn 3} |
| 75 |
\rho a^{(n)}\hackscore{i}=\sigma^{(n-1)}\hackscore{ij,j} |
| 76 |
\end{eqnarray} |
| 77 |
and boundary conditions |
| 78 |
\begin{eqnarray} \label{WAVE natural at n} |
| 79 |
\sigma^{(n-1)}\hackscore{ij}n\hackscore{j}=0 |
| 80 |
\end{eqnarray} |
| 81 |
derived from \eqn{WAVE natural} where |
| 82 |
\begin{eqnarray} \label{WAVE dyn 3a} |
| 83 |
\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}). |
| 84 |
\end{eqnarray} |
| 85 |
We also need to apply the constraint |
| 86 |
\begin{eqnarray} \label{WAVE dyn 3aa} |
| 87 |
a^{(n)}\hackscore{0}(x\hackscore C,t)= U\hackscore{0} |
| 88 |
\frac{\sqrt(2.)}{\alpha^2} (4\frac{(t-t\hackscore{0})^3}{\alpha^3}-6\frac{t-t\hackscore{0}}{\alpha})e^{\frac{1}{2}-\frac{(t-t\hackscore{0})^2}{\alpha^2}} |
| 89 |
\end{eqnarray} |
| 90 |
which is derived from equation~\ref{WAVE source} by calculating the second order time derivative, |
| 91 |
see Figure~\ref{WAVE FIG 1.1}. Now we have converted our problem for displacement, $u^{(n)}$, into a problem for |
| 92 |
acceleration, $a^{(n)}$, which now depends |
| 93 |
on the solution at the previous two time steps, $u^{(n-1)}$ and $u^{(n-2)}$. |
| 94 |
|
| 95 |
In each time step we have to solve this problem to get the acceleration $a^{(n)}$, and we will |
| 96 |
use the \LinearPDE class of the \linearPDEs to do so. The general form of the PDE defined through |
| 97 |
the \LinearPDE class is discussed in \Sec{SEC LinearPDE}. The form which is relevant here is |
| 98 |
\begin{equation}\label{WAVE dyn 100} |
| 99 |
D\hackscore{ij} a^{(n)}\hackscore{j} = - X\hackscore{ij,j}\; . |
| 100 |
\end{equation} |
| 101 |
The natural boundary condition |
| 102 |
\begin{equation}\label{WAVE dyn 101} |
| 103 |
n\hackscore{j}X\hackscore{ij} =0 |
| 104 |
\end{equation} |
| 105 |
is used. |
| 106 |
With $u=a^{(n)}$ we can identify the values to be assigned to $D$ and $X$: |
| 107 |
\begin{equation}\label{WAVE dyn 6} |
| 108 |
\begin{array}{ll} |
| 109 |
D\hackscore{ij}=\rho \delta\hackscore{ij}& |
| 110 |
X\hackscore{ij}=-\sigma^{(n-1)}\hackscore{ij} \\ |
| 111 |
\end{array} |
| 112 |
\end{equation} |
| 113 |
Moreover we need to define the location $r$ where the constraint~\ref{WAVE dyn 3aa} is applied. We will apply |
| 114 |
the constraint on a small sphere of radius $R$ around $x\hackscore C$ (we will use 3p.c. of the width of the domain): |
| 115 |
\begin{equation}\label{WAVE dyn 6 1} |
| 116 |
q\hackscore{i}(x) = |
| 117 |
\left\{ |
| 118 |
\begin{array}{rc} |
| 119 |
1 & \|x-x_c\|\le R \\ |
| 120 |
0 & \mbox{otherwise} |
| 121 |
\end{array} |
| 122 |
\right. |
| 123 |
\end{equation} |
| 124 |
The following script defines a the function \function{wavePropagation} which |
| 125 |
implements the Verlet scheme to solve our wave propagation problem. |
| 126 |
The argument \var{domain} which is a \Domain class object |
| 127 |
defines the domain of the problem. \var{h} and \var{tend} are the time step size |
| 128 |
and the end time of the simulation. \var{lam}, \var{mu} and |
| 129 |
\var{rho} are material properties. |
| 130 |
\begin{python} |
| 131 |
def wavePropagation(domain,h,tend,lam,mu,rho, x_c, src_radius, U0): |
| 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() |
| 136 |
# ... open new PDE ... |
| 137 |
mypde=LinearPDE(domain) |
| 138 |
mypde.getSolverOptions().setSolverMethod(mypde.getSolverOptions().LUMPING) |
| 139 |
kronecker=identity(mypde.getDim()) |
| 140 |
dunit=numpy.array([1.,0.,0.]) # defines direction of point source |
| 141 |
mypde.setValue(D=kronecker*rho, q=whereNegative(length(x-xc)-src_radius)*dunit) |
| 142 |
# ... set initial values .... |
| 143 |
n=0 |
| 144 |
# for first two time steps |
| 145 |
u=Vector(0.,Solution(domain)) |
| 146 |
u_last=Vector(0.,Solution(domain)) |
| 147 |
t=0 |
| 148 |
# define the location of the point source |
| 149 |
L=Locator(domain,xc) |
| 150 |
# find potential at point source |
| 151 |
u_pc=L.getValue(u) |
| 152 |
print "u at point charge=",u_pc |
| 153 |
# open file to save displacement at point source |
| 154 |
u_pc_data=FileWriter('./data/U_pc.out') |
| 155 |
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: |
| 158 |
t+=h |
| 159 |
# ... get current stress .... |
| 160 |
g=grad(u) |
| 161 |
stress=lam*trace(g)*kronecker+mu*(g+transpose(g)) |
| 162 |
# ... get new acceleration .... |
| 163 |
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) |
| 166 |
a=mypde.getSolution() |
| 167 |
# ... get new displacement ... |
| 168 |
u_new=2*u-u_last+h**2*a |
| 169 |
# ... shift displacements .... |
| 170 |
u_last=u |
| 171 |
u=u_new |
| 172 |
n+=1 |
| 173 |
print n,"-th time step t ",t |
| 174 |
u_pc=L.getValue(u) |
| 175 |
print "u at point charge=",u_pc |
| 176 |
# save displacements at point source to file for t > 0 |
| 177 |
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 |
| 181 |
if n==1 or n%10==0: saveVTK("./data/usoln.%i.vtu"%(n/10), \ |
| 182 |
acceleration=length(a)/9.81, |
| 183 |
displacement = length(u), \ |
| 184 |
tensor = stress, Ux = u[0] ) |
| 185 |
|
| 186 |
return ts, u_pc0, u_pc1, u_pc2 |
| 187 |
\end{python} |
| 188 |
Notice that |
| 189 |
all coefficients of the PDE which are independent of time $t$ are set outside the \code{while} |
| 190 |
loop. This is very efficient since it allows the \LinearPDE class to reuse information as much as possible |
| 191 |
when iterating over time. |
| 192 |
|
| 193 |
The statement |
| 194 |
\begin{python} |
| 195 |
mypde.getSolverOptions().setSolverMethod(mypde.getSolverOptions().LUMPING) |
| 196 |
\end{python} |
| 197 |
switches on the use of an aggressive approximation of the PDE operator as a diagonal matrix |
| 198 |
formed from the coefficient \var{D}. |
| 199 |
The approximation allows, at the cost of |
| 200 |
additional error, very fast |
| 201 |
solution of the PDE. When using lumping the presence of \var{A}, \var{B} or \var{C} will produce wrong results. |
| 202 |
|
| 203 |
There are a few new \escript functions in this example: |
| 204 |
\code{grad(u)} returns the gradient $u\hackscore{i,j}$ of $u$ (in fact \var{grad(g)[i,j]}=$=u\hackscore{i,j}$). |
| 205 |
There are restrictions on the argument of the \function{grad} function, for instance |
| 206 |
the statement \code{grad(grad(u))} will raise an exception. |
| 207 |
\code{trace(g)} returns the sum of the main diagonal elements \var{g[k,k]} of \var{g} |
| 208 |
and \code{transpose(g)} returns the matrix transpose of \var{g} (ie. $\var{transpose(g)[i,j]}=\var{g[j,i]}$). |
| 209 |
|
| 210 |
We initialize the values of \code{u} and \code{u_last} to be zero. It is important |
| 211 |
to initialize both with the \SolutionFS \FunctionSpace as they have to be seen as solutions of PDEs from previous time steps. In fact, the \function{grad} does not accept arguments with a certain \FunctionSpace, for more details see \Sec{SEC ESCRIPT DATA}. |
| 212 |
|
| 213 |
The \class{Locator} is designed to extract values at a given location (in this case $x\hackscore C$) from functions such as the displacement vector \code{u}. Typically the \class{Locator} is used in the following form: |
| 214 |
\begin{python} |
| 215 |
L=Locator(domain,xc) |
| 216 |
u=... |
| 217 |
u_pc=L.getValue(u) |
| 218 |
\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 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 |
access these data. |
| 223 |
|
| 224 |
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). |
| 226 |
The problem becomes so simple because we use the stress from the last time step rather then the stress which is |
| 227 |
actually present at the current time step. Schemes using this approach are called an explicit time integration |
| 228 |
schemes \index{explicit scheme} \index{time integration!explicit}. The |
| 229 |
backward Euler scheme we have used in \Chap{DIFFUSION CHAP} is |
| 230 |
an example of an implicit scheme |
| 231 |
\index{implicit scheme} \index{time integration!implicit}. In this case one uses the actual status of |
| 232 |
each variable at a particular time rather then values from previous time steps. This will lead to a problem which is |
| 233 |
more expensive to solve, in particular for non-linear problems. |
| 234 |
Although |
| 235 |
explicit time integration schemes are cheap to finalize a single time step, they need significantly smaller time |
| 236 |
steps then implicit schemes and can suffer from stability problems. Therefore they need a |
| 237 |
very careful selection of the time step size $h$. |
| 238 |
|
| 239 |
An easy, heuristic way of choosing an appropriate |
| 240 |
time step size is the Courant condition \index{Courant condition} \index{explicit scheme!Courant condition} |
| 241 |
which says that within a time step a information should not travel further than a cell used in the |
| 242 |
discretization scheme. In the case of the wave equation the velocity of a (p-) wave is given as |
| 243 |
$\sqrt{\frac{\lambda+2\mu}{\rho}}$ so one should choose $h$ from |
| 244 |
\begin{eqnarray}\label{WAVE dyn 66} |
| 245 |
h= \frac{1}{5} \sqrt{\frac{\rho}{\lambda+2\mu}} \Delta x |
| 246 |
\end{eqnarray} |
| 247 |
where $\Delta x$ is the cell diameter. The factor $\frac{1}{5}$ is a safety factor considering the heuristics of |
| 248 |
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 |
| 266 |
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}). |
| 268 |
The \code{ne} gives the number of elements in $x\hackscore{0}$ and $x\hackscore 1$ directions. |
| 269 |
The depth of the block is aligned with the $x\hackscore{2}$-direction. |
| 270 |
The depth (\code{l2}) of the block in |
| 271 |
the $x\hackscore{2}$-direction is chosen so that there are 10 elements and the magnitude of the |
| 272 |
of the depth is chosen such that the elements |
| 273 |
become cubic. We chose 10 for the number of elements in $x\hackscore{2}$-direction so that the |
| 274 |
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 |
| 275 |
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]$. |
| 276 |
\begin{python} |
| 277 |
from esys.finley import Brick |
| 278 |
ne=32 # number of cells in x_0 and x_1 directions |
| 279 |
width=10000. # length in x_0 and x_1 directions |
| 280 |
lam=3.462e9 |
| 281 |
mu=3.462e9 |
| 282 |
rho=1154. |
| 283 |
tend=60 |
| 284 |
U0=1. # amplitude of point source |
| 285 |
# spherical source at middle of bottom face |
| 286 |
xc=[width/2.,width/2.,0.] |
| 287 |
# define small radius around point xc |
| 288 |
src_radius = 0.03*width |
| 289 |
print "src_radius = ",src_radius |
| 290 |
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()) |
| 292 |
print "time step size = ",h |
| 293 |
ts, u_pc0, u_pc1, u_pc2 = \ |
| 294 |
wavePropagation(mydomain,h,tend,lam,mu,rho,xc, src_radius, U0) |
| 295 |
\end{python} |
| 296 |
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. |
| 298 |
|
| 299 |
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: |
| 301 |
\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 |
| 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. |
| 307 |
|
| 308 |
\begin{figure}[t!] |
| 309 |
\centerline{\includegraphics[width=4.in, angle=-90]{figures/WavePC}} |
| 310 |
\caption{Amplitude at Point source from the Simulation.} |
| 311 |
\label{WAVE FIG 1} |
| 312 |
\end{figure} |
| 313 |
|
| 314 |
It remains to show some possibilities to inspect the collected data \var{u_pc0}, \var{u_pc1} and \var{u_pc2}. |
| 315 |
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 |
file \file{./data/U_pc.out}: |
| 317 |
\begin{python} |
| 318 |
u_pc_data=FileWriter('./data/U_pc.out') |
| 319 |
for i in xrange(len(ts)) : |
| 320 |
u_pc_data.write("%f %f %f %f\n"%(ts[i],u_pc0[i],u_pc1[i],u_pc2[i])) |
| 321 |
u_pc_data.close() |
| 322 |
\end{python} |
| 323 |
The \code{U_pc.out} stores 4 columns of data: $t,u\hackscore x,u\hackscore y,u\hackscore z$ |
| 324 |
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 |
plotted easily using any plotting package. In \gnuplot the command: |
| 326 |
\begin{verbatim} |
| 327 |
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 |
| 328 |
lw 2, '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 |
|
Parent Directory
| 
Revision Log