1 


2 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
3 
% 
% 
4 
% Copyright (c) 20032009 by University of Queensland 
% Copyright (c) 20032010 by University of Queensland 
5 
% Earth Systems Science Computational Center (ESSCC) 
% Earth Systems Science Computational Center (ESSCC) 
6 
% http://www.uq.edu.au/esscc 
% http://www.uq.edu.au/esscc 
7 
% 
% 
11 
% 
% 
12 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
13 


14 
\section{Two Dimensional Heat Diffusion for a basic Magmatic Intrusion} 
\begin{figure}[t] 

\sslist{twodheatdiff001.py and cblib.py} 


%\label{Sec:2DHD} 


Building upon our success from the 1D models it is now prudent to expand our domain by another dimension. For this example we will be using a very simple magmatic intrusion as the basis for our model. The simulation will be a single event where some molten granite has formed a hemisphericle dome at the base of some cold sandstone country rock. A hemisphere is symmetric so taking a crosssection through its centre will effectively model a 3D problem in 2D. New concepts will include nonlinear boundaries and the ability to prescribe location specific variables. 





\begin{figure}[h!] 

15 
\centerline{\includegraphics[width=4.in]{figures/twodheatdiff}} 
\centerline{\includegraphics[width=4.in]{figures/twodheatdiff}} 
16 
\caption{2D model: granitic intrusion of sandstone country rock.} 
\caption{Example 3: 2D model: granitic intrusion of sandstone country rock.} 
17 
\label{fig:twodhdmodel} 
\label{fig:twodhdmodel} 
18 
\end{figure} 
\end{figure} 
19 


20 
To expand upon our 1D problem, the domain must first be expanded. This will be done in our definition phase by creating a square domain in $x$ and $y$ that is 600 meters along each side \reffig{fig:twodhdmodel}. The number of discrete spatial cells will be 100. The radius of the intrusion will be 200 meters And the location of the centre of the intrusion will be at the 300 meter mark on the xaxis. The domain variables are; 
\sslist{example03a.py and cblib.py} 
21 
\begin{verbatim} 

22 

Building upon our success from the 1D models, it is now prudent to expand our domain by another dimension. 
23 

For this example we use a very simple magmatic intrusion as the basis for our model. The simulation will be a single event where some molten granite has formed a cylindrical dome at the base of some cold sandstone country rock. Assuming that the cylinder is very long 
24 

we model a crosssection as shown in \reffig{fig:twodhdmodel}. We will implement the same 
25 

diffusion model as we have use for the granite blocks in \refSec{Sec:1DHDv00} 
26 

but will add the second spatial dimension and show how to define 
27 

variables depending on the location of the domain. 
28 

We use \file{onedheatdiff001b.py} as the starting point for develop this model. 
29 


30 

\section{Example 3: Two Dimensional Heat Diffusion for a basic Magmatic Intrusion} 
31 

\label{Sec:2DHD} 
32 


33 

To expand upon our 1D problem, the domain must first be expanded. In fact, we have solved a two dimensional problem already but essentially ignored the second dimension. In our definition phase we create a square domain in $x$ and $y$\footnote{In \esc the notation 
34 

$x\hackscore{0}$ and $x\hackscore{1}$ is used for $x$ and $y$, respectively.} that is $600$ meters along each side \reffig{fig:twodhdmodel}. Now we set the number of discrete spatial cells to 150 in both direction and the radius of the intrusion to $200$ meters with the centre located at the $300$ meter mark on the $x$axis. Thus, the domain variables are; 
35 

\begin{python} 
36 
mx = 600*m #meters  model length 
mx = 600*m #meters  model length 
37 
my = 600*m #meters  model width 
my = 600*m #meters  model width 
38 
ndx = 100 #mesh steps in x direction 
ndx = 150 #mesh steps in x direction 
39 
ndy = 100 #mesh steps in y direction 
ndy = 150 #mesh steps in y direction 
40 
r = 200*m #meters  radius of intrusion 
r = 200*m #meters  radius of intrusion 
41 
ic = [300, 0] #centre of intrusion (meters) 
ic = [300*m, 0] #coordinates of the centre of intrusion (meters) 
42 
q=0.*Celsius #our heat source temperature is now zero 
qH=0.*J/(sec*m**3) #our heat source temperature is zero 
43 
\end{verbatim} 
\end{python} 
44 
The next step is to define our variables for each material in the model in a manner similar to the previous tutorial. Note that each material has its own unique set of values. The time steps and set up for the domain remain as before. Prior to setting up the PDE the boundary between the two materials must be established. The distance $s$ between two points in car Cartesian coordinates is defined as: 
As before we use 
45 

\begin{python} 
46 

model = Rectangle(l0=mx,l1=my,n0=ndx, n1=ndy) 
47 

\end{python} 
48 

to generate the domain. 
49 


50 

There are two fundamental changes to the PDE that we have discussed in \refSec{Sec:1DHDv00}. Firstly, 
51 

because the material coefficients for granite and sandstone are different, we need to deal with 
52 

PDE coefficients which vary with their location in the domain. Secondly, we need 
53 

to deal with the second spatial dimension. We can investigate these two modification at the same time. 
54 

In fact, the temperature model \refEq{eqn:hd} we utilised in \refSec{Sec:1DHDv00} applied for the 
55 

1D case with a constant material parameter only. For the more general case examined in this chapter, the correct model equation is 
56 

\begin{equation} 
57 

\rho c\hackscore p \frac{\partial T}{\partial t}  \frac{\partial }{\partial x} \kappa \frac{\partial T}{\partial x}  \frac{\partial }{\partial y} \kappa \frac{\partial T}{\partial y} = q\hackscore H 
58 

\label{eqn:hd2} 
59 

\end{equation} 
60 

Notice, that for the 1D case we have $\frac{\partial T}{\partial y}=0$ and 
61 

for the case of constant material parameters $\frac{\partial }{\partial x} \kappa = \kappa \frac{\partial }{\partial x}$ thus this new equation coincides with a simplified model equation for this case. It is more convenient 
62 

to write this equation using the $\nabla$ notation as we have already seen in \refEq{eqn:commonform nabla}; 
63 

\begin{equation}\label{eqn:Tform nabla} 
64 

\rho c\hackscore p \frac{\partial T}{\partial t} 
65 

\nabla \cdot \kappa \nabla T = q\hackscore H 
66 

\end{equation} 
67 

We can easily apply the backward Euler scheme as before to end up with 
68 

\begin{equation} 
69 

\frac{\rho c\hackscore p}{h} T^{(n)} \nabla \cdot \kappa \nabla T^{(n)} = q\hackscore H + \frac{\rho c\hackscore p}{h} T^{(n1)} 
70 

\label{eqn:hdgenf2} 
71 

\end{equation} 
72 

which is very similar to \refEq{eqn:hdgenf} used to define the temperature solution in the simple 1D case. 
73 

The difference is in the second order derivate term $\nabla \cdot \kappa \nabla T^{(n)}$. Under 
74 

the light of the more general case we need to revisit the \esc PDE form as 
75 

shown in \ref{eqn:commonform2D}. For the 2D case with variable PDE coefficients 
76 

the form needs to be read as 
77 

\begin{equation}\label{eqn:commonform2D 2} 
78 

\frac{\partial }{\partial x} A\hackscore{00}\frac{\partial u}{\partial x} 
79 

\frac{\partial }{\partial x} A\hackscore{01}\frac{\partial u}{\partial y} 
80 

\frac{\partial }{\partial y} A\hackscore{10}\frac{\partial u}{\partial x} 
81 

\frac{\partial }{\partial x} A\hackscore{11}\frac{\partial u}{\partial y} 
82 

+ Du = f 
83 

\end{equation} 
84 

So besides the settings $u=T^{(n)}$, $D = \frac{\rho c \hackscore{p}}{h}$ and 
85 

$f = q \hackscore{H} + \frac{\rho c\hackscore p}{h} T^{(n1)}$ as we have used before (see \refEq{ESCRIPT SET}) we need to set 
86 

\begin{equation}\label{eqn: kappa general} 
87 

A\hackscore{00}=A\hackscore{11}=\kappa; A\hackscore{01}=A\hackscore{10}=0 
88 

\end{equation} 
89 

The fundamental difference to the 1D case is that $A\hackscore{11}$ is not set to zero but $\kappa$, 
90 

which brings in the second dimension. Important to notice that the coefficients 
91 

of the PDE may depend on their location in the domain and does not influence the usage of the PDE form. This is very convenient as we can introduce spatial dependence to the PDE coefficients without modification to the way we call the PDE solver. 
92 


93 

A very convenient way to define the matrix $A$ in \refEq{eqn: kappa general} can be carried out using the 
94 

Kronecker $\delta$ symbol\footnote{see \url{http://en.wikipedia.org/wiki/Kronecker_delta}}. The 
95 

\esc function \verbkronecker returns this matrix; 
96 

\begin{equation} 
97 

\verbkronecker(model) = \left[ 
98 

\begin{array}{cc} 
99 

1 & 0 \\ 
100 

0 & 1 \\ 
101 

\end{array} 
102 

\right] 
103 

\end{equation} 
104 

As the argument \verbmodel represents a two dimensional domain the matrix is returned as $2 \times 2$ matrix 
105 

(In case of a threedimensional domain a $3 \times 3$ matrix is returned). The call 
106 

\begin{python} 
107 

mypde.setValue(A=kappa*kronecker(model),D=rhocp/h) 
108 

\end{python} 
109 

sets the PDE coefficients according to \refEq{eqn: kappa general}. 
110 


111 

We need to check the boundary conditions before we turn to the question: how we set $\kappa$. As 
112 

pointed out in \refEq{NEUMAN 2} makes certain assumptions on the boundary conditions. In our case 
113 

this assumptions translates to; 
114 

\begin{equation} 
115 

n \cdot \kappa \nabla T^{(n)} = 
116 

n\hackscore{0} \cdot \kappa \frac{\partial T^{(n)}}{\partial x}  n\hackscore{1} \cdot \kappa \frac{\partial T^{(n)}}{\partial y} = 0 
117 

\label{eq:hom flux} 
118 

\end{equation} 
119 

which sets the normal component of the heat flux $ \kappa \cdot (\frac{\partial T^{(n)}}{\partial x}, \frac{\partial T^{(n)}}{\partial y})$ to zero. This means that the regions is insulated which is what we want. 
120 

On the left and right face of the domain where we have $(n\hackscore{0},n\hackscore{1} ) = (\mp 1,0)$ 
121 

this means $\frac{\partial T^{(n)}}{\partial x}=0$ and on the top and bottom on the domain 
122 

where we have $(n\hackscore{0},n\hackscore{1} ) = (\pm 1,0)$ this is $\frac{\partial T^{(n)}}{\partial y}=0$. 
123 


124 

\section{Setting Variable PDE Coefficients} 
125 

Now we need to look into the problem of how we define the material coefficients 
126 

$\kappa$ and $\rho c\hackscore p$ depending on their location in the domain. 
127 

We can make use of the technique used in the granite block example in \refSec{Sec:1DHDv00} 
128 

to set up the initial temperature. However, 
129 

the situation is more complicated here as we have to deal with a 
130 

curved interface between the two subdomain. 
131 


132 

Prior to setting up the PDE, the interface between the two materials must be established. 
133 

The distance $s\ge 0$ between two points $[x,y]$ and $[x\hackscore{0},y\hackscore{0}]$ in Cartesian coordinates is defined as: 
134 
\begin{equation} 
\begin{equation} 
135 
(x_{1}x_{0})^{2}+(y_{1}y_{0})^{2} = s^{2} 
(xx\hackscore{0})^{2}+(yy\hackscore{0})^{2} = s^{2} 
136 
\end{equation} 
\end{equation} 
137 
If $[x_{0},y_{0}]$ is the point $c$ the centre of the semicircle that defines our intrusion then for all the points $[x,y]$ in our solution space we can define a distance to $c$. Hence, and points that fall within the radius $r$ of our intrusion will have a corresponding value $s < r$ and all those belonging to the country rock will have a value $s > r$. By subtracting $r$ from both of these conditions we find $sr < 0$ for all intrusion points and $sr > 0$ for all country rock points. Defining these conditions within the script is quite simple and is done using the following command: 
If we define the point $[x\hackscore{0},y\hackscore{0}]$ as $ic$ which denotes the centre of the semicircle of our intrusion, then for all the points $[x,y]$ in our model we can calculate a distance to $ic$. 
138 
\begin{verbatim} 
All the points that fall within the given radius $r$ of our intrusion will have a corresponding 
139 

value $s < r$ and all those belonging to the country rock will have a value $s > r$. By subtracting $r$ from both of these conditions we find $sr < 0$ for all intrusion points and $sr > 0$ 
140 

for all country rock points. 
141 

Defining these conditions within the script is quite simple and is done using the following command: 
142 

\begin{python} 
143 
bound = length(xic)r #where the boundary will be located 
bound = length(xic)r #where the boundary will be located 
144 
\end{verbatim} 
\end{python} 
145 
This definition of the boundary can now be used with the \verb wherePositive() and \verb whereNegative() commands from before to help define the material constants and temperatures in our domain. By examining the general form we solved in the earlier tutorials, it is obvious that both \verb A and \verb D depend on the predefined variables. To set these variables accordingly and complete our PDE we use: 
This definition of the boundary can now be used with \verbwhereNegative command again to help define the material constants and temperatures in our domain. 
146 
\begin{verbatim} 
If \verbkappai and \verbkappac are the 
147 
A = (kappai)*whereNegative(bound)+(kappac)*wherePositive(bound) 
thermal conductivities for the intrusion material granite and for the surrounding sandstone, then we set; 
148 
D = (rhocpi/h)*whereNegative(bound)+(rhocpc/h)*wherePositive(bound) 
\begin{python} 
149 

x=Function(model).getX() 
150 
mypde.setValue(A=A*kronecker(model),D=D,d=eta,y=eta*Tc) 
bound = length(xic)r 
151 
\end{verbatim} 
kappa = kappai * whereNegative(bound) + kappac * (1whereNegative(bound)) 
152 
Our PDE has now been properly established. The initial conditions for temperature are set out in a similar matter: 
mypde.setValue(A=kappa*kronecker(model)) 
153 
\begin{verbatim} 
\end{python} 
154 

Notice that we are using the sample points of the \verbFunction function space as expected for the 
155 

PDE coefficient \verbA\footnote{For the experienced user: use \texttt{x=mypde.getFunctionSpace("A").getX()}.} 
156 

The corresponding statements are used to set $\rho c\hackscore p$. 
157 


158 

Our PDE has now been properly established. The initial conditions for temperature are set out in a similar manner: 
159 

\begin{python} 
160 
#defining the initial temperatures. 
#defining the initial temperatures. 
161 
T= Ti*whereNegative(bound)+Tc*wherePositive(bound) 
x=Solution(model).getX() 
162 
\end{verbatim} 
bound = length(xic)r 
163 
The iteration process now begins as before, but using our new conditions for \verb D as defined above. 
T= Ti*whereNegative(bound)+Tc*(1whereNegative(bound)) 
164 

\end{python} 
165 
\subsection{Contouring escript data} 
where \verbTi and \verbTc are the initial temperature 
166 
It is possible to contour our solution using \modmpl . Unfortunately the \modmpl contouring function only accepts regularly gridded data. As our solution is not regularly gridded, it is necessary to interpolate our solution onto a regular grid. First we extract the model coordinates using \verb getX these are then transformed to a numpy array and into individual $x$ and $y$ arrays. We also need to generate our regular grid which is done using the \modnumpy function \verb linspace . 
in the regions of the granite and surrounding sandstone, respectively. It is important to 
167 
\begin{verbatim} 
notice that we reset \verbx and \verbbound to refer to the appropriate 
168 
# rearrage mymesh to suit solution function space for contouring 
sample points of a PDE solution\footnote{For the experienced user: use \texttt{x=mypde.getFunctionSpace("r").getX()}.}. 
169 
oldspacecoords=model.getX() 

170 
coords=Data(oldspacecoords, T.getFunctionSpace()) 
\begin{figure}[ht] 
171 
coordX, coordY = toXYTuple(coords) 
\centerline{\includegraphics[width=4.in]{figures/heatrefraction001.png}} 
172 

\centerline{\includegraphics[width=4.in]{figures/heatrefraction030.png}} 
173 

\centerline{\includegraphics[width=4.in]{figures/heatrefraction200.png}} 
174 

\caption{Example 3a: 2D model: Total temperature distribution ($T$) at time step $1$, $20$ and $200$. Contour lines show temperature.} 
175 

\label{fig:twodhdans} 
176 

\end{figure} 
177 


178 

\section{Contouring \esc data using \modmpl.} 
179 

\label{Sec:2DHD plot} 
180 

To complete our transition from a 1D to a 2D model we also need to modify the 
181 

plotting procedure. As before we use the \modmpl to do the plotting 
182 

in this case a contour plot. For 2D contour plots \modmpl expects that the 
183 

data are regularly gridded. We have no control over the location and ordering of the sample points 
184 

used to represent the solution. Therefore it is necessary to interpolate our solution onto a regular grid. 
185 


186 

In \refSec{sec: plot T} we have already learned how to extract the $x$ coordinates of sample points as 
187 

\verbnumpy array to hand the values to \modmpl. This can easily be extended to extract both the 
188 

$x$ and the $y$ coordinates; 
189 

\begin{python} 
190 

import numpy as np 
191 

def toXYTuple(coords): 
192 

coords = np.array(coords.toListOfTuples()) #convert to Tuple 
193 

coordX = coords[:,0] #X components. 
194 

coordY = coords[:,1] #Y components. 
195 

return coordX,coordY 
196 

\end{python} 
197 

For convenience we have put this function into \file{clib.py} file so we can use this 
198 

function in other scripts more easily. 
199 


200 


201 

We now generate a regular $100 \times 100$ grid over the domain ($mx$ and $my$ 
202 

are the dimensions in the $x$ and $y$ directions) which is done using the \modnumpy function \verblinspace . 
203 

\begin{python} 
204 

from clib import toXYTuple 
205 

# get sample points for temperature as for contouring 
206 

coordX, coordY = toXYTuple(T.getFunctionSpace().getX()) 
207 
# create regular grid 
# create regular grid 
208 
xi = np.linspace(0.0,mx,100) 
xi = np.linspace(0.0,mx,75) 
209 
yi = np.linspace(0.0,my,100) 
yi = np.linspace(0.0,my,75) 
210 
\end{verbatim} 
\end{python} 
211 
The remainder of our contouring commands reside within the \verb while loop so that a new contour is generated for each time step. For each time step the solution much be regridded for \modmpl using the \verb griddata function. This function interprets an irregular grid and solution from \verb tempT , \verb xi and \verb yi this is transformed to the new coordinates defined by \verb coordX and \verb coordY with an output \verb zi . It is now possible to use the \verb contourf function which generates a colour filled contour. The colour gradient of our plots is set via the command \verb pl.matplotlib.pyplot.autumn() , other colours are listed on the \modmpl web page. Our results are then contoured, visually adjusted using the \modmpl functions and then saved to file. \verb pl.clf() clears the figure in readiness for the next time iteration. 
The values \verb[xi[k], yi[l]] are the grid points. 
212 


213 
\begin{verbatim} 
The remainder of our contouring commands reside within a \verb while loop so that a new contour is generated for each time step. For each time step the solution must be regridded for \modmpl using the \verb griddata function. This function interprets irregularly located values \verb tempT at locations defined by \verb coordX and \verb coordY as values at the new coordinates of a rectangular grid defined by 
214 

\verb xi and \verb yi . The output is \verb zi . It is now possible to use the \verb contourf function which generates colour filled contours. The colour gradient of our plots is set via the command \verb pl.matplotlib.pyplot.autumn() , other colours are listed on the \modmpl web page\footnote{see \url{http://matplotlib.sourceforge.net/api/}}. Our results are then contoured, visually adjusted using the \modmpl functions and then saved to file. \verb pl.clf() clears the figure in readiness for the next time iteration. 
215 

\begin{python} 
216 
#grid the data. 
#grid the data. 
217 
zi = pl.matplotlib.mlab.griddata(coordX,coordY,tempT,xi,yi) 
zi = pl.matplotlib.mlab.griddata(coordX,coordY,tempT,xi,yi) 
218 
# contour the gridded data, plotting dots at the randomly spaced data points. 
# contour the gridded data, plotting dots at the randomly spaced data points. 
224 
pl.title("Heat diffusion from an intrusion.") 
pl.title("Heat diffusion from an intrusion.") 
225 
pl.xlabel("Horizontal Displacement (m)") 
pl.xlabel("Horizontal Displacement (m)") 
226 
pl.ylabel("Depth (m)") 
pl.ylabel("Depth (m)") 
227 
pl.savefig(os.path.join(save_path,"heatrefraction%03d.png") %i) 
pl.savefig(os.path.join(save_path,"Tcontour%03d.png") %i) 
228 
pl.clf() 
pl.clf() 
229 
\end{verbatim} 
\end{python} 
230 

The function \verbpl.contour is used to add labeled contour lines to the plot. 
231 

The results for selected time steps are shown in \reffig{fig:twodhdans}. 
232 



\begin{figure}[h!] 


\centerline{\includegraphics[width=4.in]{figures/heatrefraction050}} 


\caption{2D model: Total temperature distribution ($T$) at time $t=50$.} 


\label{fig:twodhdmodel} 


\end{figure} 

233 


234 


235 

\section{Advanced Visualization using VTK} 
236 


237 

\sslist{example03b.py} 
238 

An alternative approach to \modmpl for visualization is the usage of a package which base on 
239 

visualization tool kit (VTK) library\footnote{see \url{http://www.vtk.org/}}. There are variety 
240 

of packages available. Here we use the package \mayavi\footnote{see \url{http://code.enthought.com/projects/mayavi/}} as an example. 
241 


242 

\mayavi is an interactive, GUI driven tool which is 
243 

really designed to visualize large three dimensional data sets where \modmpl 
244 

is not suitable. But it is very useful when it comes to two dimensional problems. 
245 

The decision of which tool is the best can be subjective and the user should determine which package they require and are most comfortable with. The main difference between using \mayavi (and other VTK based tools) 
246 

or \modmpl is that the actual visualization is detached from the 
247 

calculation by writing the results to external files 
248 

and importing them into \mayavi. In 3D the best camera position for rendering a scene is not obvious 
249 

before the results are available. Therefore the user may need to try 
250 

different position before the best is found. Moreover, in many cases a 3D interactive 
251 

visualization is the only way to really understand the results (e.g. using stereographic projection). 
252 


253 

To write the temperatures at each time step to data files in the VTK file format one 
254 

needs to insert a \verbsaveVTK call into the code; 
255 

\begin{python} 
256 

while t<=tend: 
257 

i+=1 #counter 
258 

t+=h #current time 
259 

mypde.setValue(Y=qH+T*rhocp/h) 
260 

T=mypde.getSolution() 
261 

saveVTK(os.path.join(save_path,"data.%03d.vtu"%i, T=T) 
262 

\end{python} 
263 

The data files, eg. \file{data.001.vtu}, contains all necessary information to 
264 

visualize the temperature and can directly processed by \mayavi. Notice that there is no 
265 

regridding required. It is recommended to use the file extension \file{.vtu} for files 
266 

created by \verbsaveVTK. 
267 


268 

\begin{figure}[ht] 
269 

\centerline{\includegraphics[width=4.in]{figures/ScreeshotMayavi2n1}} 
270 

\caption{Example 3b: \mayavi start up Window.} 
271 

\label{fig:mayavi window} 
272 

\end{figure} 
273 


274 

\begin{figure}[ht] 
275 

\centerline{\includegraphics[width=4.in]{figures/ScreeshotMayavi2n2}} 
276 

\caption{Example 3b: \mayavi data control window.} 
277 

\label{fig:mayavi window2} 
278 

\end{figure} 
279 

After you run the script you will find the 
280 

result files \file{data.*.vtu} in the result directory \file{data/example03}. Run the 
281 

command 
282 

\begin{python} 
283 

>> mayavi2 d data.001.vtu m Surface & 
284 

\end{python} 
285 

from the result directory. \mayavi will start up a window similar to \reffig{fig:mayavi window}. 
286 

The right hand side shows the temperature at the first time step. To show 
287 

the results at other time steps you can click at the item \texttt{VTK XML file (data.001.vtu) (timeseries)} 
288 

at the top left hand side. This will bring up a new window similar to the window shown in \reffig{fig:mayavi window2}. By clicking at the arrows in the top right corner you can move forwards and backwards in time. 
289 

You will also notice the text \textbf{T} next to the item \texttt{Point scalars name}. The 
290 

name \textbf{T} corresponds to the keyword argument name \texttt{T} we have used 
291 

in the \verbsaveVTK call. In this menu item you can select other results 
292 

you may have written to the output file for visualization. 
293 


294 

\textbf{For the advanced user}: Using the \modmpl to visualize spatially distributed data 
295 

is not MPI compatible. However, the \verbsaveVTK function can be used with MPI. In fact, 
296 

the function collects the values of the sample points spread across processor ranks into a single. 
297 

For more details on writing scripts for parallel computing please consult the \emph{user's guide}. 