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{2D model: granitic intrusion of sandstone country rock.} 
17 
\label{fig:twodhdmodel} 
18 
\end{figure} 
19 

20 
\sslist{twodheatdiff001.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. 
23 
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 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{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 didn't put much 
34 
attention to the second dimension. This will be changed now. 
35 
In our definition phase by creating a square domain in $x$ and $y$\footnote{In \esc the notation 
36 
$x\hackscore{0}$ and $x\hackscore{1}$ is used for $x$ and $y$, respectively.} that is $600$ meters along each side \reffig{fig:twodhdmodel}. The number of discrete spatial cells will be 100 in either direction. The radius of the intrusion will be $200$ meters with the centre located at the $300$ meter mark on the $x$axis. The domain variables are; 
37 
\begin{python} 
38 
mx = 600*m #meters  model length 
39 
my = 600*m #meters  model width 
40 
ndx = 150 #mesh steps in x direction 
41 
ndy = 150 #mesh steps in y direction 
42 
r = 200*m #meters  radius of intrusion 
43 
ic = [300*m, 0] #coordinates of the centre of intrusion (meters) 
44 
qH=0.*J/(sec*m**3) #our heat source temperature is zero 
45 
\end{python} 
46 
As before we use 
47 
\begin{python} 
48 
model = Rectangle(l0=mx,l1=my,n0=ndx, n1=ndy) 
49 
\end{python} 
50 
to generate the domain. 
51 

52 
There are two fundamental changes to the PDE has we have discussed PDEs in \refSec{Sec:1DHDv00}. Firstly, 
53 
as the material coefficients for granite and sandstone is different, we need to deal with 
54 
PDE coefficients which vary with there location in the domain. Secondly, we need 
55 
to deal with the second spatial dimension. We will look at these two modification at the same time. 
56 
In fact our temperature model \refEq{eqn:hd} we have used in \refSec{Sec:1DHDv00} applies for the 
57 
1D case with constant material parameter only. For the more general case we are interested 
58 
in this chapter the correct model equation is 
59 
\begin{equation} 
60 
\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 
61 
\label{eqn:hd2} 
62 
\end{equation} 
63 
Notice, that for the 1D case we have $\frac{\partial T}{\partial y}=0$ and 
64 
for the case of constant material parameters $\frac{\partial }{\partial x} \kappa = \kappa \frac{\partial }{\partial x}$ so this new equation coincides with simplified model equation for this case. It is more convenient 
65 
to write this equation using the $\nabla$ notation as we have already seen in \refEq{eqn:commonform nabla}; 
66 
\begin{equation}\label{eqn:Tform nabla} 
67 
\rho c\hackscore p \frac{\partial T}{\partial t} 
68 
\nabla \cdot \kappa \nabla T = q\hackscore H 
69 
\end{equation} 
70 
We can easily apply the backward Euler scheme as before to end up with 
71 
\begin{equation} 
72 
\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)} 
73 
\label{eqn:hdgenf2} 
74 
\end{equation} 
75 
which is very similar to \refEq{eqn:hdgenf} used to define the temperature update in the simple 1D case. 
76 
The difference is in the second order derivate term $\nabla \cdot \kappa \nabla T^{(n)}$. Under 
77 
the light of the more general case we need to revisit the \esc PDE form as 
78 
shown in \ref{eqn:commonform2D}. For the 2D case with variable PDE coefficients 
79 
the form needs to be read as 
80 
\begin{equation}\label{eqn:commonform2D 2} 
81 
\frac{\partial }{\partial x} A\hackscore{00}\frac{\partial u}{\partial x} 
82 
\frac{\partial }{\partial x} A\hackscore{01}\frac{\partial u}{\partial y} 
83 
\frac{\partial }{\partial y} A\hackscore{10}\frac{\partial u}{\partial x} 
84 
\frac{\partial }{\partial x} A\hackscore{11}\frac{\partial u}{\partial y} 
85 
+ Du = f 
86 
\end{equation} 
87 
So besides the settings $u=T^{(n)}$, $D = \frac{\rho c \hackscore{p}}{h}$ and 
88 
$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 
89 
\begin{equation}\label{eqn: kappa general} 
90 
A\hackscore{00}=A\hackscore{11}=\kappa; A\hackscore{01}=A\hackscore{10}=0 
91 
\end{equation} 
92 
The fundamental difference to the 1D case is that $A\hackscore{11}$ is not set to zero but $\kappa$ 
93 
which brings in the second dimension. Important to notice that the fact that the coefficients 
94 
of the PDE may depend on their location in the domain now 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. 
95 

96 
A very convenient way to define the matrix $A$ is required in \refEq{eqn: kappa general} is using the 
97 
Kronecker $\delta$ symbol\footnote{see \url{http://en.wikipedia.org/wiki/Kronecker_delta}}. The 
98 
\esc function \verbkronecker returns this matrix; 
99 
\begin{equation} 
100 
\verbkronecker(model) = \left[ 
101 
\begin{array}{cc} 
102 
1 & 0 \\ 
103 
0 & 1 \\ 
104 
\end{array} 
105 
\right] 
106 
\end{equation} 
107 
As the argument \verbmodel represents a two dimensional domain the matrix is returned as $2 \times 2$ matrix 
108 
(In case of a threedimensional domain a $3 \times 3$ matrix is returned). The call 
109 
\begin{python} 
110 
mypde.setValue(A=kappa*kronecker(model),D=rhocp/h) 
111 
\end{python} 
112 
sets the PDE coefficients according to \refEq{eqn: kappa general}. 
113 

114 
Before we turn the question how we set $\kappa$ we need to check the boundary conditions. As 
115 
pointed out in \refEq{NEUMAN 2} makes certain assumptions on the boundary conditions. In our case 
116 
this assumptions translates to; 
117 
\begin{equation} 
118 
n \cdot \kappa \nabla T^{(n)} = 
119 
n\hackscore{0} \cdot \kappa \frac{\partial T^{(n)}}{\partial x}  n\hackscore{1} \cdot \kappa \frac{\partial T^{(n)}}{\partial y} = 0 
120 
\label{eq:hom flux} 
121 
\end{equation} 
122 
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. 
123 
On the left and right face of the domain where we have $(n\hackscore{0},n\hackscore{1} ) = (\mp 1,0)$ 
124 
this means $\frac{\partial T^{(n)}}{\partial x}=0$ and on the top and bottom on the domain 
125 
where we have $(n\hackscore{0},n\hackscore{1} ) = (\pm 1,0)$ this is $\frac{\partial T^{(n)}}{\partial y}=0$. 
126 

127 
\section{Setting Variable PDE Coefficients} 
128 
Now we need to look into the problem how we define the material coefficients 
129 
$\kappa$ and $\rho c\hackscore p$ depending on there location in the domain. 
130 
We have used the technique we discuss here already when we set up the initial 
131 
temperature in the granite block example in \refSec{Sec:1DHDv00}. However, 
132 
the situation is more complicated here as we have to deal with a 
133 
curved interface between the two subdomain. 
134 

135 
Prior to setting up the PDE the interface between the two materials must be established. 
136 
The distance $s\ge 0$ between two points $[x,y]$ and $[x\hackscore{0},y\hackscore{0}]$ in Cartesian coordinates is defined as: 
137 
\begin{equation} 
138 
(xx\hackscore{0})^{2}+(yy\hackscore{0})^{2} = s^{2} 
139 
\end{equation} 
140 
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$. 
141 
All the points that fall within the given radius $r$ of our intrusion will have a corresponding 
142 
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$ 
143 
for all country rock points. 
144 
Defining these conditions within the script is quite simple and is done using the following command: 
145 
\begin{python} 
146 
bound = length(xic)r #where the boundary will be located 
147 
\end{python} 
148 
This definition of the boundary can now be used with \verbwhereNegative command again to help define the material constants and temperatures in our domain. 
149 
If \verbkappai and \verbkappac are the 
150 
thermal conductivities for the intrusion material granite and for the surrounding sandstone we set; 
151 
\begin{python} 
152 
x=Function(model).getX() 
153 
bound = length(xic)r 
154 
kappa = kappai * whereNegative(bound) + kappac * (1whereNegative(bound)) 
155 
mypde.setValue(A=kappa*kronecker(model)) 
156 
\end{python} 
157 
Notice that we are using the sample points of the \verbFunction function space as expected for the 
158 
PDE coefficient \verbA\footnote{For the experience user: use \texttt{x=mypde.getFunctionSpace("A").getX()}.} 
159 
The corresponding statements are used to set $\rho c\hackscore p$. 
160 

161 
Our PDE has now been properly established. The initial conditions for temperature are set out in a similar matter: 
162 
\begin{python} 
163 
#defining the initial temperatures. 
164 
x=Solution(model).getX() 
165 
bound = length(xic)r 
166 
T= Ti*whereNegative(bound)+Tc*(1whereNegative(bound)) 
167 
\end{python} 
168 
where \verbTi and \verbTc are the initial temperature 
169 
in the regions of the granite and surrounding sandstone, respectively. It is important to 
170 
notice that we have reset \verbx and \verbbound to refer to the appropriate 
171 
sample points of a PDE solution\footnote{For the experience user: use \texttt{x=mypde.getFunctionSpace("r").getX()}.}. 
172 

173 
\begin{figure}[h] 
174 
\centerline{\includegraphics[width=4.in]{figures/heatrefraction001.png}} 
175 
\centerline{\includegraphics[width=4.in]{figures/heatrefraction030.png}} 
176 
\centerline{\includegraphics[width=4.in]{figures/heatrefraction200.png}} 
177 
\caption{2D model: Total temperature distribution ($T$) at time step $1$, $20$ and $200$. Contour lines show temperature.} 
178 
\label{fig:twodhdans} 
179 
\end{figure} 
180 

181 
\section{Contouring \esc data using \modmpl.} 
182 
\label{Sec:2DHD plot} 
183 
To complete our transition from a 1D to a 2D model we also need to modify the 
184 
plotting procedure. As before we use the \modmpl to do the plotting 
185 
in this case a contour plot. For 2D contour plots \modmpl expects that the 
186 
data are regularly gridded. We have no control on the location and ordering of the sample points 
187 
used to represent the solution. Therefore it is necessary to interpolate our solution onto a regular grid. 
188 

189 
In \refSec{sec: plot T} we have already learned how to extract the $x$ coordinates of sample points as 
190 
\verbnumpy array to hand the values to \modmpl. This can easily be extended to extract both the 
191 
$x$ and the $y$ coordinates; 
192 
\begin{python} 
193 
import numpy as np 
194 
def toXYTuple(coords): 
195 
coords = np.array(coords.toListOfTuples()) #convert to Tuple 
196 
coordX = coords[:,0] #X components. 
197 
coordY = coords[:,1] #Y components. 
198 
return coordX,coordY 
199 
\end{python} 
200 
For convenience we have put this function into \file{clib.py} file so we can use this 
201 
function in other scripts more easily. 
202 

203 

204 
We now generate a regular $100 \times 100$ grid over the domain ($mx$ and $my$ 
205 
are the dimensions in $x$ and $y$ direction) which is done using the \modnumpy function \verblinspace . 
206 
\begin{python} 
207 
from clib import toXYTuple 
208 
# get sample points for temperature as for contouring 
209 
coordX, coordY = toXYTuple(T.getFunctionSpace().getX()) 
210 
# create regular grid 
211 
xi = np.linspace(0.0,mx,75) 
212 
yi = np.linspace(0.0,my,75) 
213 
\end{python} 
214 
The values \verb[xi[k], yi[l]] are the grid points. 
215 

216 
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 a potentially 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 
217 
\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. 
218 
\begin{python} 
219 
#grid the data. 
220 
zi = pl.matplotlib.mlab.griddata(coordX,coordY,tempT,xi,yi) 
221 
# contour the gridded data, plotting dots at the randomly spaced data points. 
222 
pl.matplotlib.pyplot.autumn() 
223 
pl.contourf(xi,yi,zi,10) 
224 
CS = pl.contour(xi,yi,zi,5,linewidths=0.5,colors='k') 
225 
pl.clabel(CS, inline=1, fontsize=8) 
226 
pl.axis([0,600,0,600]) 
227 
pl.title("Heat diffusion from an intrusion.") 
228 
pl.xlabel("Horizontal Displacement (m)") 
229 
pl.ylabel("Depth (m)") 
230 
pl.savefig(os.path.join(save_path,"heatrefraction%03d.png") %i) 
231 
pl.clf() 
232 
\end{python} 
233 
The function \verbpl.contour is used to add labeled contour lines to the plot. 
234 
The results for selected time steps are shown in \reffig{fig:twodhdans}. 
235 

236 

237 

238 
\section{Advanced Visualization using VTK} 
239 

240 
\sslist{twodheatdiffvtk.py} 
241 
An alternative approach to \modmpl for visualization is the usage of a package which base on 
242 
visualization tool kit (VTK) library\footnote{see \url{http://www.vtk.org/}}. There is a variety 
243 
of package available. Here we will use the package \mayavi\footnote{see \url{http://code.enthought.com/projects/mayavi/}} as an example. 
244 

245 
\mayavi is an interactive, GUI driven tool which is 
246 
really designed to visualize large three dimensional data sets where \modmpl 
247 
is not suitable. But it is very useful when it comes to two dimensional problems. 
248 
The decision which tool is best is finally the user's decision. The main 
249 
difference between using \mayavi (and other VTK based tools) 
250 
or \modmpl is the fact that actually visualization is detached from the 
251 
calculation by writing the results to external files 
252 
and import them into \mayavi. In 3D where the best camera position for rendering a scene is not obvious 
253 
before the results are available. Therefore the user may need to try 
254 
different position before the best is found. Moreover, in many cases in 3D the interactive 
255 
visualization is the only way to really understand the results (e.g. using stereographic projection). 
256 

257 
To write the temperatures at each time step to data files in the VTK file format one 
258 
needs to insert a \verbsaveVTK call into the code; 
259 
\begin{python} 
260 
while t<=tend: 
261 
i+=1 #counter 
262 
t+=h #current time 
263 
mypde.setValue(Y=qH+T*rhocp/h) 
264 
T=mypde.getSolution() 
265 
saveVTK(os.path.join(save_path,"data.%03d.vtu"%i, T=T) 
266 
\end{python} 
267 
The data files, eg. \file{data.001.vtu}, contains all necessary information to 
268 
visualize the temperature and can directly processed by \mayavi. Notice that there is no 
269 
regridding required. It is recommended to use the file extension \file{.vtu} for files 
270 
created by \verbsaveVTK. 
271 

272 
\begin{figure}[ht] 
273 
\centerline{\includegraphics[width=4.in]{figures/ScreeshotMayavi2n1}} 
274 
\caption{\mayavi start up Window.} 
275 
\label{fig:mayavi window} 
276 
\end{figure} 
277 

278 
\begin{figure}[ht] 
279 
\centerline{\includegraphics[width=4.in]{figures/ScreeshotMayavi2n2}} 
280 
\caption{\mayavi data control window.} 
281 
\label{fig:mayavi window2} 
282 
\end{figure} 
283 
After you have run the script you will find the 
284 
result files \file{data.*.vtu} in the result directory \file{data/twodheatdiff}. Run the 
285 
command 
286 
\begin{python} 
287 
>> mayavi2 d data.001.vtu m Surface & 
288 
\end{python} 
289 
from the result directory. \mayavi will start up a window similar to \reffig{fig:mayavi window}. 
290 
The right hand side shows the temperature at the first time step. To show 
291 
the results at other time steps you can click at the item \texttt{VTK XML file (data.001.vtu) (timeseries)} 
292 
at the top left hand side. This will bring up a new window similar to tye 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. 
293 
You will also notice the text \textbf{T} next to the item \texttt{Point scalars name}. The 
294 
name \textbf{T} corresponds to the keyword argument name \texttt{T} we have used 
295 
in the \verbsaveVTK call. In this menu item you can select other results 
296 
you may have written to the output file for visualization. 
297 

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