19 


20 
\sslist{example03a.py and cblib.py} 
\sslist{example03a.py and cblib.py} 
21 


22 
Building upon our success from the 1D models, it is now prudent to expand our domain by another dimension. 
Building upon our success from the 1D models, it is now prudent to expand our 
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 
domain by another dimension. 
24 
we model a crosssection as shown in \reffig{fig:twodhdmodel}. We will implement the same 
For this example we use a very simple magmatic intrusion as the basis for our 
25 

model. The simulation will be a single event where some molten granite has 
26 

formed a cylindrical dome at the base of some cold sandstone country rock. 
27 

Assuming that the cylinder is very long 
28 

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

the same 
30 
diffusion model as we have use for the granite blocks in \refSec{Sec:1DHDv00} 
diffusion model as we have use for the granite blocks in \refSec{Sec:1DHDv00} 
31 
but will add the second spatial dimension and show how to define 
but will add the second spatial dimension and show how to define 
32 
variables depending on the location of the domain. 
variables depending on the location of the domain. 
33 
We use \file{onedheatdiff001b.py} as the starting point for develop this model. 
We use \file{onedheatdiff001b.py} as the starting point for develop this model. 
34 


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

Intrusion} 
37 
\label{Sec:2DHD} 
\label{Sec:2DHD} 
38 


39 
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 
To expand upon our 1D problem, the domain must first be expanded. In fact, we 
40 
$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; 
have solved a two dimensional problem already but essentially ignored the second 
41 

dimension. In our definition phase we create a square domain in $x$ and 
42 

$y$\footnote{In \esc the notation 
43 

$x\hackscore{0}$ and $x\hackscore{1}$ is used for $x$ and $y$, respectively.} 
44 

that is $600$ meters along each side \reffig{fig:twodhdmodel}. Now we set the 
45 

number of discrete spatial cells to 150 in both direction and the radius of the 
46 

intrusion to $200$ meters with the centre located at the $300$ meter mark on the 
47 

$x$axis. Thus, the domain variables are; 
48 
\begin{python} 
\begin{python} 
49 
mx = 600*m #meters  model length 
mx = 600*m #meters  model length 
50 
my = 600*m #meters  model width 
my = 600*m #meters  model width 
60 
\end{python} 
\end{python} 
61 
to generate the domain. 
to generate the domain. 
62 


63 
There are two fundamental changes to the PDE that we have discussed in \refSec{Sec:1DHDv00}. Firstly, 
There are two fundamental changes to the PDE that we have discussed in 
64 
because the material coefficients for granite and sandstone are different, we need to deal with 
\refSec{Sec:1DHDv00}. Firstly, 
65 
PDE coefficients which vary with their location in the domain. Secondly, we need 
because the material coefficients for granite and sandstone are different, we 
66 
to deal with the second spatial dimension. We can investigate these two modification at the same time. 
need to deal with 
67 
In fact, the temperature model \refEq{eqn:hd} we utilised in \refSec{Sec:1DHDv00} applied for the 
PDE coefficients which vary with their location in the domain. Secondly, we 
68 
1D case with a constant material parameter only. For the more general case examined in this chapter, the correct model equation is 
need 
69 

to deal with the second spatial dimension. We can investigate these two 
70 

modification at the same time. 
71 

In fact, the temperature model \refEq{eqn:hd} we utilised in 
72 

\refSec{Sec:1DHDv00} applied for the 
73 

1D case with a constant material parameter only. For the more general case 
74 

examined in this chapter, the correct model equation is 
75 
\begin{equation} 
\begin{equation} 
76 
\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 
\rho c\hackscore p \frac{\partial T}{\partial t}  \frac{\partial }{\partial x} 
77 

\kappa \frac{\partial T}{\partial x}  \frac{\partial }{\partial y} \kappa 
78 

\frac{\partial T}{\partial y} = q\hackscore H 
79 
\label{eqn:hd2} 
\label{eqn:hd2} 
80 
\end{equation} 
\end{equation} 
81 
Notice, that for the 1D case we have $\frac{\partial T}{\partial y}=0$ and 
Notice, that for the 1D case we have $\frac{\partial T}{\partial y}=0$ and 
82 
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 
for the case of constant material parameters $\frac{\partial }{\partial x} 
83 
to write this equation using the $\nabla$ notation as we have already seen in \refEq{eqn:commonform nabla}; 
\kappa = \kappa \frac{\partial }{\partial x}$ thus this new equation coincides 
84 

with a simplified model equation for this case. It is more convenient 
85 

to write this equation using the $\nabla$ notation as we have already seen in 
86 

\refEq{eqn:commonform nabla}; 
87 
\begin{equation}\label{eqn:Tform nabla} 
\begin{equation}\label{eqn:Tform nabla} 
88 
\rho c\hackscore p \frac{\partial T}{\partial t} 
\rho c\hackscore p \frac{\partial T}{\partial t} 
89 
\nabla \cdot \kappa \nabla T = q\hackscore H 
\nabla \cdot \kappa \nabla T = q\hackscore H 
90 
\end{equation} 
\end{equation} 
91 
We can easily apply the backward Euler scheme as before to end up with 
We can easily apply the backward Euler scheme as before to end up with 
92 
\begin{equation} 
\begin{equation} 
93 
\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)} 
\frac{\rho c\hackscore p}{h} T^{(n)} \nabla \cdot \kappa \nabla T^{(n)} = 
94 

q\hackscore H + \frac{\rho c\hackscore p}{h} T^{(n1)} 
95 
\label{eqn:hdgenf2} 
\label{eqn:hdgenf2} 
96 
\end{equation} 
\end{equation} 
97 
which is very similar to \refEq{eqn:hdgenf} used to define the temperature solution in the simple 1D case. 
which is very similar to \refEq{eqn:hdgenf} used to define the temperature 
98 
The difference is in the second order derivate term $\nabla \cdot \kappa \nabla T^{(n)}$. Under 
update in the simple 1D case. 
99 

The difference is in the second order derivate term $\nabla \cdot \kappa \nabla 
100 

T^{(n)}$. Under 
101 
the light of the more general case we need to revisit the \esc PDE form as 
the light of the more general case we need to revisit the \esc PDE form as 
102 
shown in \ref{eqn:commonform2D}. For the 2D case with variable PDE coefficients 
shown in \ref{eqn:commonform2D}. For the 2D case with variable PDE coefficients 
103 
the form needs to be read as 
the form needs to be read as 
109 
+ Du = f 
+ Du = f 
110 
\end{equation} 
\end{equation} 
111 
So besides the settings $u=T^{(n)}$, $D = \frac{\rho c \hackscore{p}}{h}$ and 
So besides the settings $u=T^{(n)}$, $D = \frac{\rho c \hackscore{p}}{h}$ and 
112 
$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 
$f = q \hackscore{H} + \frac{\rho c\hackscore p}{h} T^{(n1)}$ as we have used 
113 

before (see \refEq{ESCRIPT SET}) we need to set 
114 
\begin{equation}\label{eqn: kappa general} 
\begin{equation}\label{eqn: kappa general} 
115 
A\hackscore{00}=A\hackscore{11}=\kappa; A\hackscore{01}=A\hackscore{10}=0 
A\hackscore{00}=A\hackscore{11}=\kappa; A\hackscore{01}=A\hackscore{10}=0 
116 
\end{equation} 
\end{equation} 
117 
The fundamental difference to the 1D case is that $A\hackscore{11}$ is not set to zero but $\kappa$, 
The fundamental difference to the 1D case is that $A\hackscore{11}$ is not set 
118 

to zero but $\kappa$, 
119 
which brings in the second dimension. Important to notice that the coefficients 
which brings in the second dimension. Important to notice that the coefficients 
120 
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. 
of the PDE may depend on their location in the domain and does not influence the 
121 

usage of the PDE form. This is very convenient as we can introduce spatial 
122 
A very convenient way to define the matrix $A$ in \refEq{eqn: kappa general} can be carried out using the 
dependence to the PDE coefficients without modification to the way we call the 
123 
Kronecker $\delta$ symbol\footnote{see \url{http://en.wikipedia.org/wiki/Kronecker_delta}}. The 
PDE solver. 
124 


125 

A very convenient way to define the matrix $A$ in \refEq{eqn: kappa general} can 
126 

be carried out using the 
127 

Kronecker $\delta$ symbol\footnote{see 
128 

\url{http://en.wikipedia.org/wiki/Kronecker_delta}}. The 
129 
\esc function \verbkronecker returns this matrix; 
\esc function \verbkronecker returns this matrix; 
130 
\begin{equation} 
\begin{equation} 
131 
\verbkronecker(model) = \left[ 
\verbkronecker(model) = \left[ 
135 
\end{array} 
\end{array} 
136 
\right] 
\right] 
137 
\end{equation} 
\end{equation} 
138 
As the argument \verbmodel represents a two dimensional domain the matrix is returned as $2 \times 2$ matrix 
As the argument \verbmodel represents a two dimensional domain the matrix is 
139 
(In case of a threedimensional domain a $3 \times 3$ matrix is returned). The call 
returned as $2 \times 2$ matrix 
140 

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

call 
142 
\begin{python} 
\begin{python} 
143 
mypde.setValue(A=kappa*kronecker(model),D=rhocp/h) 
mypde.setValue(A=kappa*kronecker(model),D=rhocp/h) 
144 
\end{python} 
\end{python} 
145 
sets the PDE coefficients according to \refEq{eqn: kappa general}. 
sets the PDE coefficients according to \refEq{eqn: kappa general}. 
146 


147 
We need to check the boundary conditions before we turn to the question: how we set $\kappa$. As 
We need to check the boundary conditions before we turn to the question: how we 
148 
pointed out in \refEq{NEUMAN 2} makes certain assumptions on the boundary conditions. In our case 
set $\kappa$. As 
149 

pointed out in \refEq{NEUMAN 2} makes certain assumptions on the boundary 
150 

conditions. In our case 
151 
this assumptions translates to; 
this assumptions translates to; 
152 
\begin{equation} 
\begin{equation} 
153 
n \cdot \kappa \nabla T^{(n)} = 
n \cdot \kappa \nabla T^{(n)} = 
154 
n\hackscore{0} \cdot \kappa \frac{\partial T^{(n)}}{\partial x}  n\hackscore{1} \cdot \kappa \frac{\partial T^{(n)}}{\partial y} = 0 
n\hackscore{0} \cdot \kappa \frac{\partial T^{(n)}}{\partial x}  
155 

n\hackscore{1} \cdot \kappa \frac{\partial T^{(n)}}{\partial y} = 0 
156 
\label{eq:hom flux} 
\label{eq:hom flux} 
157 
\end{equation} 
\end{equation} 
158 
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. 
which sets the normal component of the heat flux $ \kappa \cdot (\frac{\partial 
159 
On the left and right face of the domain where we have $(n\hackscore{0},n\hackscore{1} ) = (\mp 1,0)$ 
T^{(n)}}{\partial x}, \frac{\partial T^{(n)}}{\partial y})$ to zero. This means 
160 
this means $\frac{\partial T^{(n)}}{\partial x}=0$ and on the top and bottom on the domain 
that the regions is insulated which is what we want. 
161 
where we have $(n\hackscore{0},n\hackscore{1} ) = (\pm 1,0)$ this is $\frac{\partial T^{(n)}}{\partial y}=0$. 
On the left and right face of the domain where we have 
162 

$(n\hackscore{0},n\hackscore{1} ) = (\mp 1,0)$ 
163 

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

the domain 
165 

where we have $(n\hackscore{0},n\hackscore{1} ) = (\pm 1,0)$ this is 
166 

$\frac{\partial T^{(n)}}{\partial y}=0$. 
167 


168 
\section{Setting Variable PDE Coefficients} 
\section{Setting Variable PDE Coefficients} 
169 
Now we need to look into the problem of how we define the material coefficients 
Now we need to look into the problem of how we define the material coefficients 
170 
$\kappa$ and $\rho c\hackscore p$ depending on their location in the domain. 
$\kappa$ and $\rho c\hackscore p$ depending on their location in the domain. 
171 
We can make use of the technique used in the granite block example in \refSec{Sec:1DHDv00} 
We can make use of the technique used in the granite block example in 
172 

\refSec{Sec:1DHDv00} 
173 
to set up the initial temperature. However, 
to set up the initial temperature. However, 
174 
the situation is more complicated here as we have to deal with a 
the situation is more complicated here as we have to deal with a 
175 
curved interface between the two subdomain. 
curved interface between the two subdomain. 
176 


177 
Prior to setting up the PDE, the interface between the two materials must be established. 
Prior to setting up the PDE, the interface between the two materials must be 
178 
The distance $s\ge 0$ between two points $[x,y]$ and $[x\hackscore{0},y\hackscore{0}]$ in Cartesian coordinates is defined as: 
established. 
179 

The distance $s\ge 0$ between two points $[x,y]$ and 
180 

$[x\hackscore{0},y\hackscore{0}]$ in Cartesian coordinates is defined as: 
181 
\begin{equation} 
\begin{equation} 
182 
(xx\hackscore{0})^{2}+(yy\hackscore{0})^{2} = s^{2} 
(xx\hackscore{0})^{2}+(yy\hackscore{0})^{2} = s^{2} 
183 
\end{equation} 
\end{equation} 
184 
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$. 
If we define the point $[x\hackscore{0},y\hackscore{0}]$ as $ic$ which denotes 
185 
All the points that fall within the given radius $r$ of our intrusion will have a corresponding 
the centre of the semicircle of our intrusion, then for all the points $[x,y]$ 
186 
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$ 
in our model we can calculate a distance to $ic$. 
187 

All the points that fall within the given radius $r$ of our intrusion will have 
188 

a corresponding 
189 

value $s < r$ and all those belonging to the country rock will have a value $s > 
190 

r$. By subtracting $r$ from both of these conditions we find $sr < 0$ for all 
191 

intrusion points and $sr > 0$ 
192 
for all country rock points. 
for all country rock points. 
193 
Defining these conditions within the script is quite simple and is done using the following command: 
Defining these conditions within the script is quite simple and is done using 
194 

the following command: 
195 
\begin{python} 
\begin{python} 
196 
bound = length(xic)r #where the boundary will be located 
bound = length(xic)r #where the boundary will be located 
197 
\end{python} 
\end{python} 
198 
This definition of the boundary can now be used with \verbwhereNegative command again to help define the material constants and temperatures in our domain. 
This definition of the boundary can now be used with \verbwhereNegative 
199 

command again to help define the material constants and temperatures in our 
200 

domain. 
201 
If \verbkappai and \verbkappac are the 
If \verbkappai and \verbkappac are the 
202 
thermal conductivities for the intrusion material granite and for the surrounding sandstone, then we set; 
thermal conductivities for the intrusion material granite and for the 
203 

surrounding sandstone, then we set; 
204 
\begin{python} 
\begin{python} 
205 
x=Function(model).getX() 
x=Function(model).getX() 
206 
bound = length(xic)r 
bound = length(xic)r 
207 
kappa = kappai * whereNegative(bound) + kappac * (1whereNegative(bound)) 
kappa = kappai * whereNegative(bound) + kappac * (1whereNegative(bound)) 
208 
mypde.setValue(A=kappa*kronecker(model)) 
mypde.setValue(A=kappa*kronecker(model)) 
209 
\end{python} 
\end{python} 
210 
Notice that we are using the sample points of the \verbFunction function space as expected for the 
Notice that we are using the sample points of the \verbFunction function space 
211 
PDE coefficient \verbA\footnote{For the experienced user: use \texttt{x=mypde.getFunctionSpace("A").getX()}.} 
as expected for the 
212 

PDE coefficient \verbA\footnote{For the experienced user: use 
213 

\texttt{x=mypde.getFunctionSpace("A").getX()}.} 
214 
The corresponding statements are used to set $\rho c\hackscore p$. 
The corresponding statements are used to set $\rho c\hackscore p$. 
215 


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

temperature are set out in a similar manner: 
218 
\begin{python} 
\begin{python} 
219 
#defining the initial temperatures. 
#defining the initial temperatures. 
220 
x=Solution(model).getX() 
x=Solution(model).getX() 
222 
T= Ti*whereNegative(bound)+Tc*(1whereNegative(bound)) 
T= Ti*whereNegative(bound)+Tc*(1whereNegative(bound)) 
223 
\end{python} 
\end{python} 
224 
where \verbTi and \verbTc are the initial temperature 
where \verbTi and \verbTc are the initial temperature 
225 
in the regions of the granite and surrounding sandstone, respectively. It is important to 
in the regions of the granite and surrounding sandstone, respectively. It is 
226 

important to 
227 
notice that we reset \verbx and \verbbound to refer to the appropriate 
notice that we reset \verbx and \verbbound to refer to the appropriate 
228 
sample points of a PDE solution\footnote{For the experienced user: use \texttt{x=mypde.getFunctionSpace("r").getX()}.}. 
sample points of a PDE solution\footnote{For the experienced user: use 
229 

\texttt{x=mypde.getFunctionSpace("r").getX()}.}. 
230 


231 
\begin{figure}[ht] 
\begin{figure}[ht] 
232 
\centerline{\includegraphics[width=4.in]{figures/heatrefraction001.png}} 
\centerline{\includegraphics[width=4.in]{figures/heatrefraction001.png}} 
233 
\centerline{\includegraphics[width=4.in]{figures/heatrefraction030.png}} 
\centerline{\includegraphics[width=4.in]{figures/heatrefraction030.png}} 
234 
\centerline{\includegraphics[width=4.in]{figures/heatrefraction200.png}} 
\centerline{\includegraphics[width=4.in]{figures/heatrefraction200.png}} 
235 
\caption{Example 3a: 2D model: Total temperature distribution ($T$) at time step $1$, $20$ and $200$. Contour lines show temperature.} 
\caption{Example 3a: 2D model: Total temperature distribution ($T$) at time step 
236 

$1$, $20$ and $200$. Contour lines show temperature.} 
237 
\label{fig:twodhdans} 
\label{fig:twodhdans} 
238 
\end{figure} 
\end{figure} 
239 


242 
To complete our transition from a 1D to a 2D model we also need to modify the 
To complete our transition from a 1D to a 2D model we also need to modify the 
243 
plotting procedure. As before we use the \modmpl to do the plotting 
plotting procedure. As before we use the \modmpl to do the plotting 
244 
in this case a contour plot. For 2D contour plots \modmpl expects that the 
in this case a contour plot. For 2D contour plots \modmpl expects that the 
245 
data are regularly gridded. We have no control over the location and ordering of the sample points 
data are regularly gridded. We have no control over the location and ordering of 
246 
used to represent the solution. Therefore it is necessary to interpolate our solution onto a regular grid. 
the sample points 
247 

used to represent the solution. Therefore it is necessary to interpolate our 
248 
In \refSec{sec: plot T} we have already learned how to extract the $x$ coordinates of sample points as 
solution onto a regular grid. 
249 
\verbnumpy array to hand the values to \modmpl. This can easily be extended to extract both the 

250 

In \refSec{sec: plot T} we have already learned how to extract the $x$ 
251 

coordinates of sample points as 
252 

\verbnumpy array to hand the values to \modmpl. This can easily be extended to 
253 

extract both the 
254 
$x$ and the $y$ coordinates; 
$x$ and the $y$ coordinates; 
255 
\begin{python} 
\begin{python} 
256 
import numpy as np 
import numpy as np 
260 
coordY = coords[:,1] #Y components. 
coordY = coords[:,1] #Y components. 
261 
return coordX,coordY 
return coordX,coordY 
262 
\end{python} 
\end{python} 
263 
For convenience we have put this function into \file{clib.py} file so we can use this 
For convenience we have put this function into \file{clib.py} file so we can use 
264 

this 
265 
function in other scripts more easily. 
function in other scripts more easily. 
266 


267 


268 
We now generate a regular $100 \times 100$ grid over the domain ($mx$ and $my$ 
We now generate a regular $100 \times 100$ grid over the domain ($mx$ and $my$ 
269 
are the dimensions in the $x$ and $y$ directions) which is done using the \modnumpy function \verblinspace . 
are the dimensions in the $x$ and $y$ directions) which is done using the 
270 

\modnumpy function \verblinspace . 
271 
\begin{python} 
\begin{python} 
272 
from clib import toXYTuple 
from clib import toXYTuple 
273 
# get sample points for temperature as for contouring 
# get sample points for temperature as for contouring 
278 
\end{python} 
\end{python} 
279 
The values \verb[xi[k], yi[l]] are the grid points. 
The values \verb[xi[k], yi[l]] are the grid points. 
280 


281 
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 
The remainder of our contouring commands reside within a \verb while loop so 
282 
\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. 
that a new contour is generated for each time step. For each time step the 
283 

solution must be regridded for \modmpl using the \verb griddata function. This 
284 

function interprets irregularly located values \verb tempT at locations defined 
285 

by \verb coordX and \verb coordY as values at the new coordinates of a 
286 

rectangular grid defined by 
287 

\verb xi and \verb yi . The output is \verb zi . It is now possible to use the 
288 

\verb contourf function which generates colour filled contours. The colour 
289 

gradient of our plots is set via the command 
290 

\verb pl.matplotlib.pyplot.autumn() , other colours are listed on the \modmpl web page\footnote{see 
291 

\url{http://matplotlib.sourceforge.net/api/}}. Our results are then contoured, 
292 

visually adjusted using the \modmpl functions and then saved to file. \verb 
293 

pl.clf() clears the figure in readiness for the next time iteration. 
294 
\begin{python} 
\begin{python} 
295 
#grid the data. 
#grid the data. 
296 
zi = pl.matplotlib.mlab.griddata(coordX,coordY,tempT,xi,yi) 
zi = pl.matplotlib.mlab.griddata(coordX,coordY,tempT,xi,yi) 
306 
pl.savefig(os.path.join(save_path,"Tcontour%03d.png") %i) 
pl.savefig(os.path.join(save_path,"Tcontour%03d.png") %i) 
307 
pl.clf() 
pl.clf() 
308 
\end{python} 
\end{python} 
309 
The function \verbpl.contour is used to add labeled contour lines to the plot. 
The function \verbpl.contour is used to add labeled contour lines to the 
310 

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


313 


315 
\section{Advanced Visualization using VTK} 
\section{Advanced Visualization using VTK} 
316 


317 
\sslist{example03b.py} 
\sslist{example03b.py} 
318 
An alternative approach to \modmpl for visualization is the usage of a package which base on 
An alternative approach to \modmpl for visualization is the usage of a package 
319 
visualization tool kit (VTK) library\footnote{see \url{http://www.vtk.org/}}. There are variety 
which base on 
320 
of packages available. Here we use the package \mayavi\footnote{see \url{http://code.enthought.com/projects/mayavi/}} as an example. 
visualization tool kit (VTK) library\footnote{see \url{http://www.vtk.org/}}. 
321 

There are variety 
322 

of packages available. Here we use the package \mayavi\footnote{see 
323 

\url{http://code.enthought.com/projects/mayavi/}} as an example. 
324 


325 
\mayavi is an interactive, GUI driven tool which is 
\mayavi is an interactive, GUI driven tool which is 
326 
really designed to visualize large three dimensional data sets where \modmpl 
really designed to visualize large three dimensional data sets where \modmpl 
327 
is not suitable. But it is very useful when it comes to two dimensional problems. 
is not suitable. But it is very useful when it comes to two dimensional 
328 
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) 
problems. 
329 

The decision of which tool is the best can be subjective and the user should 
330 

determine which package they require and are most comfortable with. The main 
331 

difference between using \mayavi (and other VTK based tools) 
332 
or \modmpl is that the actual visualization is detached from the 
or \modmpl is that the actual visualization is detached from the 
333 
calculation by writing the results to external files 
calculation by writing the results to external files 
334 
and importing them into \mayavi. In 3D the best camera position for rendering a scene is not obvious 
and importing them into \mayavi. In 3D the best camera position for rendering a 
335 

scene is not obvious 
336 
before the results are available. Therefore the user may need to try 
before the results are available. Therefore the user may need to try 
337 
different position before the best is found. Moreover, in many cases a 3D interactive 
different position before the best is found. Moreover, in many cases a 3D 
338 
visualization is the only way to really understand the results (e.g. using stereographic projection). 
interactive 
339 

visualization is the only way to really understand the results (e.g. using 
340 

stereographic projection). 
341 


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

one 
344 
needs to insert a \verbsaveVTK call into the code; 
needs to insert a \verbsaveVTK call into the code; 
345 
\begin{python} 
\begin{python} 
346 
while t<=tend: 
while t<=tend: 
351 
saveVTK(os.path.join(save_path,"data.%03d.vtu"%i, T=T) 
saveVTK(os.path.join(save_path,"data.%03d.vtu"%i, T=T) 
352 
\end{python} 
\end{python} 
353 
The data files, eg. \file{data.001.vtu}, contains all necessary information to 
The data files, eg. \file{data.001.vtu}, contains all necessary information to 
354 
visualize the temperature and can directly processed by \mayavi. Notice that there is no 
visualize the temperature and can directly processed by \mayavi. Notice that 
355 
regridding required. It is recommended to use the file extension \file{.vtu} for files 
there is no 
356 

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

files 
358 
created by \verbsaveVTK. 
created by \verbsaveVTK. 
359 


360 
\begin{figure}[ht] 
\begin{figure}[ht] 
369 
\label{fig:mayavi window2} 
\label{fig:mayavi window2} 
370 
\end{figure} 
\end{figure} 
371 
After you run the script you will find the 
After you run the script you will find the 
372 
result files \file{data.*.vtu} in the result directory \file{data/example03}. Run the 
result files \file{data.*.vtu} in the result directory \file{data/example03}. 
373 

Run the 
374 
command 
command 
375 
\begin{python} 
\begin{python} 
376 
>> mayavi2 d data.001.vtu m Surface & 
>> mayavi2 d data.001.vtu m Surface & 
377 
\end{python} 
\end{python} 
378 
from the result directory. \mayavi will start up a window similar to \reffig{fig:mayavi window}. 
from the result directory. \mayavi will start up a window similar to 
379 

\reffig{fig:mayavi window}. 
380 
The right hand side shows the temperature at the first time step. To show 
The right hand side shows the temperature at the first time step. To show 
381 
the results at other time steps you can click at the item \texttt{VTK XML file (data.001.vtu) (timeseries)} 
the results at other time steps you can click at the item \texttt{VTK XML file 
382 
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. 
(data.001.vtu) (timeseries)} 
383 
You will also notice the text \textbf{T} next to the item \texttt{Point scalars name}. The 
at the top left hand side. This will bring up a new window similar to the window 
384 
name \textbf{T} corresponds to the keyword argument name \texttt{T} we have used 
shown in \reffig{fig:mayavi window2}. By clicking at the arrows in the top 
385 

right corner you can move forwards and backwards in time. 
386 

You will also notice the text \textbf{T} next to the item \texttt{Point scalars 
387 

name}. The 
388 

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

used 
390 
in the \verbsaveVTK call. In this menu item you can select other results 
in the \verbsaveVTK call. In this menu item you can select other results 
391 
you may have written to the output file for visualization. 
you may have written to the output file for visualization. 
392 



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


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


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


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

393 

\textbf{For the advanced user}: Using the \modmpl to visualize spatially 
394 

distributed data 
395 

is not MPI compatible. However, the \verbsaveVTK function can be used with 
396 

MPI. In fact, 
397 

the function collects the values of the sample points spread across processor 
398 

ranks into a single. 
399 

For more details on writing scripts for parallel computing please consult the 
400 

\emph{user's guide}. 