1 

2 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
3 
% 
4 
% Copyright (c) 20032010 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/osl3.0.php 
11 
% 
12 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
13 

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

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

22 
Building upon our success from the 1D models, it is now prudent to expand our 
23 
domain by another dimension. 
24 
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 used for the granite blocks in \refSec{Sec:1DHDv00} 
31 
but will add the second spatial dimension and show how to define 
32 
variables depending on the location of the domain. 
33 
We use \file{onedheatdiff001b.py} as the starting point to develop this model. 
34 

35 
\section{Example 3: Two Dimensional Heat Diffusion for a basic Magmatic 
36 
Intrusion} 
37 
\label{Sec:2DHD} 
38 

39 
To expand upon our 1D problem, the domain must first be expanded. In fact, we 
40 
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 directions and the radius of 
46 
the intrusion to $200$ meters with the centre located at the $300$ meter mark 
47 
on the $x$axis. Thus, the domain variables are; 
48 
\begin{python} 
49 
mx = 600*m #meters  model length 
50 
my = 600*m #meters  model width 
51 
ndx = 150 #mesh steps in x direction 
52 
ndy = 150 #mesh steps in y direction 
53 
r = 200*m #meters  radius of intrusion 
54 
ic = [300*m, 0] #coordinates of the centre of intrusion (meters) 
55 
qH=0.*J/(sec*m**3) #our heat source temperature is zero 
56 
\end{python} 
57 
As before we use 
58 
\begin{python} 
59 
model = Rectangle(l0=mx,l1=my,n0=ndx, n1=ndy) 
60 
\end{python} 
61 
to generate the domain. 
62 

63 
There are two fundamental changes to the PDE that we have discussed in 
64 
\refSec{Sec:1DHDv00}. Firstly, 
65 
because the material coefficients for granite and sandstone are different, we 
66 
need to deal with 
67 
PDE coefficients which vary with their location in the domain. Secondly, we 
68 
need to deal with the second spatial dimension. We can investigate these two 
69 
modifications at the same time. 
70 
In fact, the temperature model \refEq{eqn:hd} we utilised in 
71 
\refSec{Sec:1DHDv00} applied for the 
72 
1D case with a constant material parameter only. For the more general case 
73 
examined in this chapter, the correct model equation is 
74 
\begin{equation} 
75 
\rho c\hackscore p \frac{\partial T}{\partial t}  \frac{\partial }{\partial x} 
76 
\kappa \frac{\partial T}{\partial x}  \frac{\partial }{\partial y} \kappa 
77 
\frac{\partial T}{\partial y} = q\hackscore H 
78 
\label{eqn:hd2} 
79 
\end{equation} 
80 
Notice that for the 1D case we have $\frac{\partial T}{\partial y}=0$ and 
81 
for the case of constant material parameters $\frac{\partial }{\partial x} 
82 
\kappa = \kappa \frac{\partial }{\partial x}$ thus this new equation coincides 
83 
with a simplified model equation for this case. It is more convenient 
84 
to write this equation using the $\nabla$ notation as we have already seen in 
85 
\refEq{eqn:commonform nabla}; 
86 
\begin{equation}\label{eqn:Tform nabla} 
87 
\rho c\hackscore p \frac{\partial T}{\partial t} 
88 
\nabla \cdot \kappa \nabla T = q\hackscore H 
89 
\end{equation} 
90 
We can easily apply the backward Euler scheme as before to end up with 
91 
\begin{equation} 
92 
\frac{\rho c\hackscore p}{h} T^{(n)} \nabla \cdot \kappa \nabla T^{(n)} = 
93 
q\hackscore H + \frac{\rho c\hackscore p}{h} T^{(n1)} 
94 
\label{eqn:hdgenf2} 
95 
\end{equation} 
96 
which is very similar to \refEq{eqn:hdgenf} used to define the temperature 
97 
update in the simple 1D case. 
98 
The difference is in the second order derivative term 
99 
$\nabla \cdot \kappa \nabla T^{(n)}$. Under the light of the more general case 
100 
we need to revisit the \esc PDE form as shown in \ref{eqn:commonform2D}. 
101 
For the 2D case with variable PDE coefficients the form needs to be read as 
102 
\begin{equation}\label{eqn:commonform2D 2} 
103 
\frac{\partial }{\partial x} A\hackscore{00}\frac{\partial u}{\partial x} 
104 
\frac{\partial }{\partial x} A\hackscore{01}\frac{\partial u}{\partial y} 
105 
\frac{\partial }{\partial y} A\hackscore{10}\frac{\partial u}{\partial x} 
106 
\frac{\partial }{\partial x} A\hackscore{11}\frac{\partial u}{\partial y} 
107 
+ Du = f 
108 
\end{equation} 
109 
So besides the settings $u=T^{(n)}$, $D = \frac{\rho c \hackscore{p}}{h}$ and 
110 
$f = q \hackscore{H} + \frac{\rho c\hackscore p}{h} T^{(n1)}$ as we have used 
111 
before (see \refEq{ESCRIPT SET}) we need to set 
112 
\begin{equation}\label{eqn: kappa general} 
113 
A\hackscore{00}=A\hackscore{11}=\kappa; A\hackscore{01}=A\hackscore{10}=0 
114 
\end{equation} 
115 
The fundamental difference to the 1D case is that $A\hackscore{11}$ is not set 
116 
to zero but $\kappa$, 
117 
which brings in the second dimension. It is important to note that the 
118 
coefficients of the PDE may depend on their location in the domain which does 
119 
not influence the usage of the PDE form. This is very convenient as we can 
120 
introduce spatial dependence to the PDE coefficients without modification to 
121 
the way we call the PDE solver. 
122 

123 
A very convenient way to define the matrix $A$ in \refEq{eqn: kappa general} 
124 
can be carried out using the Kronecker $\delta$ symbol\footnote{see 
125 
\url{http://en.wikipedia.org/wiki/Kronecker_delta}}. The 
126 
\esc function \verbkronecker returns this matrix; 
127 
\begin{equation} 
128 
\verbkronecker(model) = \left[ 
129 
\begin{array}{cc} 
130 
1 & 0 \\ 
131 
0 & 1 \\ 
132 
\end{array} 
133 
\right] 
134 
\end{equation} 
135 
As the argument \verbmodel represents a two dimensional domain the matrix is 
136 
returned as a $2 \times 2$ matrix 
137 
(in the case of a threedimensional domain a $3 \times 3$ matrix is returned). 
138 
The call 
139 
\begin{python} 
140 
mypde.setValue(A=kappa*kronecker(model),D=rhocp/h) 
141 
\end{python} 
142 
sets the PDE coefficients according to \refEq{eqn: kappa general}. 
143 

144 
We need to check the boundary conditions before we turn to the question: how do 
145 
we set $\kappa$. As pointed out \refEq{NEUMAN 2} makes certain assumptions on 
146 
the boundary conditions. In our case these assumptions translate to; 
147 
\begin{equation} 
148 
n \cdot \kappa \nabla T^{(n)} = 
149 
n\hackscore{0} \cdot \kappa \frac{\partial T^{(n)}}{\partial x}  
150 
n\hackscore{1} \cdot \kappa \frac{\partial T^{(n)}}{\partial y} = 0 
151 
\label{eq:hom flux} 
152 
\end{equation} 
153 
which sets the normal component of the heat flux $ \kappa \cdot (\frac{\partial 
154 
T^{(n)}}{\partial x}, \frac{\partial T^{(n)}}{\partial y})$ to zero. This means 
155 
that the region is insulated which is what we want. 
156 
On the left and right face of the domain where we have 
157 
$(n\hackscore{0},n\hackscore{1} ) = (\mp 1,0)$ 
158 
this means $\frac{\partial T^{(n)}}{\partial x}=0$ and on the top and bottom on 
159 
the domain 
160 
where we have $(n\hackscore{0},n\hackscore{1} ) = (\pm 1,0)$ this is 
161 
$\frac{\partial T^{(n)}}{\partial y}=0$. 
162 

163 
\section{Setting variable PDE Coefficients} 
164 
Now we need to look into the problem of how we define the material coefficients 
165 
$\kappa$ and $\rho c\hackscore p$ depending on their location in the domain. 
166 
We can make use of the technique used in the granite block example in 
167 
\refSec{Sec:1DHDv00} 
168 
to set up the initial temperature. However, 
169 
the situation is more complicated here as we have to deal with a 
170 
curved interface between the two subdomains. 
171 

172 
Prior to setting up the PDE, the interface between the two materials must be 
173 
established. 
174 
The distance $s\ge 0$ between two points $[x,y]$ and 
175 
$[x\hackscore{0},y\hackscore{0}]$ in Cartesian coordinates is defined as: 
176 
\begin{equation} 
177 
(xx\hackscore{0})^{2}+(yy\hackscore{0})^{2} = s^{2} 
178 
\end{equation} 
179 
If we define the point $[x\hackscore{0},y\hackscore{0}]$ as $ic$ which denotes 
180 
the centre of the semicircle of our intrusion, then for all the points $[x,y]$ 
181 
in our model we can calculate a distance to $ic$. 
182 
All the points that fall within the given radius $r$ of our intrusion will have 
183 
a corresponding 
184 
value $s < r$ and all those belonging to the country rock will have a value $s > 
185 
r$. By subtracting $r$ from both of these conditions we find $sr < 0$ for all 
186 
intrusion points and $sr > 0$ for all country rock points. 
187 
Defining these conditions within the script is quite simple and is done using 
188 
the following command: 
189 
\begin{python} 
190 
bound = length(xic)r #where the boundary will be located 
191 
\end{python} 
192 
This definition of the boundary can now be used with the \verbwhereNegative 
193 
command again to help define the material constants and temperatures in our 
194 
domain. 
195 
If \verbkappai and \verbkappac are the 
196 
thermal conductivities for the intrusion material granite and for the 
197 
surrounding sandstone, then we set; 
198 
\begin{python} 
199 
x=Function(model).getX() 
200 
bound = length(xic)r 
201 
kappa = kappai * whereNegative(bound) + kappac * (1whereNegative(bound)) 
202 
mypde.setValue(A=kappa*kronecker(model)) 
203 
\end{python} 
204 
Notice that we are using the sample points of the \verbFunction function space 
205 
as expected for the PDE coefficient \verbA\footnote{For the experienced user: use 
206 
\texttt{x=mypde.getFunctionSpace("A").getX()}.}. 
207 
The corresponding statements are used to set $\rho c\hackscore p$. 
208 

209 
Our PDE has now been properly established. The initial conditions for 
210 
temperature are set out in a similar manner: 
211 
\begin{python} 
212 
#defining the initial temperatures. 
213 
x=Solution(model).getX() 
214 
bound = length(xic)r 
215 
T= Ti*whereNegative(bound)+Tc*(1whereNegative(bound)) 
216 
\end{python} 
217 
where \verbTi and \verbTc are the initial temperature in the regions of the 
218 
granite and surrounding sandstone, respectively. It is important to 
219 
notice that we reset \verbx and \verbbound to refer to the appropriate 
220 
sample points of a PDE solution\footnote{For the experienced user: use 
221 
\texttt{x=mypde.getFunctionSpace("r").getX()}.}. 
222 

223 
\begin{figure}[ht] 
224 
\centerline{\includegraphics[width=4.in]{figures/heatrefraction001.png}} 
225 
\centerline{\includegraphics[width=4.in]{figures/heatrefraction030.png}} 
226 
\centerline{\includegraphics[width=4.in]{figures/heatrefraction200.png}} 
227 
\caption{Example 3a: 2D model: Total temperature distribution ($T$) at time step 
228 
$1$, $20$ and $200$. Contour lines show temperature.} 
229 
\label{fig:twodhdans} 
230 
\end{figure} 
231 

232 
\section{Contouring \esc Data using \modmpl.} 
233 
\label{Sec:2DHD plot} 
234 
To complete our transition from a 1D to a 2D model we also need to modify the 
235 
plotting procedure. As before we use \modmpl to do the plotting in this case a 
236 
contour plot. For 2D contour plots \modmpl expects that the data are regularly 
237 
gridded. We have no control over the location and ordering of the sample points 
238 
used to represent the solution. Therefore it is necessary to interpolate our 
239 
solution onto a regular grid. 
240 

241 
In \refSec{sec: plot T} we have already learned how to extract the $x$ 
242 
coordinates of sample points as 
243 
\verbnumpy array to hand the values to \modmpl. This can easily be extended 
244 
to extract both the $x$ and the $y$ coordinates; 
245 
\begin{python} 
246 
import numpy as np 
247 
def toXYTuple(coords): 
248 
coords = np.array(coords.toListOfTuples()) #convert to Tuple 
249 
coordX = coords[:,0] #X components. 
250 
coordY = coords[:,1] #Y components. 
251 
return coordX,coordY 
252 
\end{python} 
253 
For convenience we have put this function into \file{clib.py} so we can use 
254 
this function in other scripts more easily. 
255 

256 
We now generate a regular $100 \times 100$ grid over the domain ($mx$ and $my$ 
257 
are the dimensions in the $x$ and $y$ directions) which is done using the 
258 
\modnumpy function \verblinspace. 
259 
\begin{python} 
260 
from clib import toXYTuple 
261 
# get sample points for temperature as for contouring 
262 
coordX, coordY = toXYTuple(T.getFunctionSpace().getX()) 
263 
# create regular grid 
264 
xi = np.linspace(0.0,mx,75) 
265 
yi = np.linspace(0.0,my,75) 
266 
\end{python} 
267 
The values \verb[xi[k], yi[l]] are the grid points. 
268 

269 
The remainder of our contouring commands resides within a \verbwhile loop so 
270 
that a new contour is generated for each time step. For each time step the 
271 
solution must be regridded for \modmpl using the \verbgriddata function. This 
272 
function interprets irregularly located values \verbtempT at locations defined 
273 
by \verbcoordX and \verbcoordY as values at the new coordinates of a 
274 
rectangular grid defined by 
275 
\verbxi and \verbyi. The output is \verbzi. It is now possible to use the 
276 
\verbcontourf function which generates colour filled contours. The colour 
277 
gradient of our plots is set via the command 
278 
\verbpl.matplotlib.pyplot.autumn(), other colours are listed on the \modmpl web page\footnote{see 
279 
\url{http://matplotlib.sourceforge.net/api/}}. Our results are then contoured, 
280 
visually adjusted using the \modmpl functions and then saved to a file. 
281 
\verbpl.clf() clears the figure in readiness for the next time iteration. 
282 
\begin{python} 
283 
#grid the data. 
284 
zi = pl.matplotlib.mlab.griddata(coordX,coordY,tempT,xi,yi) 
285 
# contour the gridded data, plotting dots at the randomly spaced data points. 
286 
pl.matplotlib.pyplot.autumn() 
287 
pl.contourf(xi,yi,zi,10) 
288 
CS = pl.contour(xi,yi,zi,5,linewidths=0.5,colors='k') 
289 
pl.clabel(CS, inline=1, fontsize=8) 
290 
pl.axis([0,600,0,600]) 
291 
pl.title("Heat diffusion from an intrusion.") 
292 
pl.xlabel("Horizontal Displacement (m)") 
293 
pl.ylabel("Depth (m)") 
294 
pl.savefig(os.path.join(save_path,"Tcontour%03d.png") %i) 
295 
pl.clf() 
296 
\end{python} 
297 
The function \verbpl.contour is used to add labelled contour lines to the 
298 
plot. 
299 
The results for selected time steps are shown in \reffig{fig:twodhdans}. 
300 

301 

302 
\section{Advanced Visualisation using VTK} 
303 

304 
\sslist{example03b.py} 
305 
An alternative approach to \modmpl for visualisation is the usage of a package 
306 
which is based on the Visualization Toolkit (VTK) library\footnote{see \url{http://www.vtk.org/}}. 
307 
There is a variety of packages available. Here we use the package \mayavi\footnote{see 
308 
\url{http://code.enthought.com/projects/mayavi/}} as an example. 
309 

310 
\mayavi is an interactive, GUI driven tool which is 
311 
really designed to visualise large three dimensional data sets where \modmpl 
312 
is not suitable. But it is also very useful when it comes to two dimensional 
313 
problems. 
314 
The decision of which tool is the best can be subjective and users should 
315 
determine which package they require and are most comfortable with. The main 
316 
difference between using \mayavi (or other VTK based tools) and \modmpl is that 
317 
the actual visualisation is detached from the calculation by writing the 
318 
results to external files and importing them into \mayavi. In 3D the best 
319 
camera position for rendering a scene is not obvious before the results are 
320 
available. Therefore the user may need to try different settings before the 
321 
best is found. Moreover, in many cases a 3D interactive visualisation is the 
322 
only way to really understand the results (e.g. using stereographic projection). 
323 

324 
To write the temperatures at each time step to data files in the VTK file format 
325 
one needs to insert a \verbsaveVTK call into the code; 
326 
\begin{python} 
327 
while t<=tend: 
328 
i+=1 #counter 
329 
t+=h #current time 
330 
mypde.setValue(Y=qH+T*rhocp/h) 
331 
T=mypde.getSolution() 
332 
saveVTK(os.path.join(save_path,"data.%03d.vtu"%i, T=T) 
333 
\end{python} 
334 
The data files, e.g. \file{data.001.vtu}, contain all necessary information to 
335 
visualise the temperature and can be directly processed by \mayavi. Note that 
336 
there is no regridding required. It is recommended to use the file extension 
337 
\file{.vtu} for files created by \verbsaveVTK. 
338 

339 
\begin{figure}[ht] 
340 
\centerline{\includegraphics[width=4.in]{figures/ScreeshotMayavi2n1}} 
341 
\caption{Example 3b: \mayavi start up window} 
342 
\label{fig:mayavi window} 
343 
\end{figure} 
344 

345 
\begin{figure}[ht] 
346 
\centerline{\includegraphics[width=4.in]{figures/ScreeshotMayavi2n2}} 
347 
\caption{Example 3b: \mayavi data control window} 
348 
\label{fig:mayavi window2} 
349 
\end{figure} 
350 
After you run the script you will find the 
351 
result files \file{data.*.vtu} in the result directory \file{data/example03}. 
352 
Run the command 
353 
\begin{python} 
354 
>> mayavi2 d data.001.vtu m Surface & 
355 
\end{python} 
356 
from the result directory. \mayavi will start up a window similar to 
357 
\reffig{fig:mayavi window}. 
358 
The right hand side shows the temperature at the first time step. To show 
359 
the results at other time steps you can click at the item \texttt{VTK XML file 
360 
(data.001.vtu) (timeseries)} 
361 
at the top left hand side. This will bring up a new window similar to the window 
362 
shown in \reffig{fig:mayavi window2}. By clicking at the arrows in the top 
363 
right corner you can move forwards and backwards in time. 
364 
You will also notice the text \textbf{T} next to the item \texttt{Point scalars 
365 
name}. The 
366 
name \textbf{T} corresponds to the keyword argument name \texttt{T} we have 
367 
used in the \verbsaveVTK call. In this menu item you can select other results 
368 
you may have written to the output file for visualisation. 
369 

370 
\textbf{For the advanced user}: Using the \modmpl to visualise spatially 
371 
distributed data is not MPI compatible. However, the \verbsaveVTK function 
372 
can be used with MPI. In fact, the function collects the values of the sample 
373 
points spread across processor ranks into a single file. 
374 
For more details on writing scripts for parallel computing please consult the 
375 
\emph{user's guide}. 
376 
