1 

2 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
3 
% Copyright (c) 20032018 by The University of Queensland 
4 
% http://www.uq.edu.au 
5 
% 
6 
% Primary Business: Queensland, Australia 
7 
% Licensed under the Apache License, version 2.0 
8 
% http://www.apache.org/licenses/LICENSE2.0 
9 
% 
10 
% Development until 2012 by Earth Systems Science Computational Center (ESSCC) 
11 
% Development 20122013 by School of Earth Sciences 
12 
% Development from 2014 by Centre for Geoscience Computing (GeoComp) 
13 
% 
14 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
15 

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

22 
\sslist{example03a.py and cblib.py} 
23 

24 
Building upon our success from the 1D models, it is now prudent to expand our 
25 
domain by another dimension. 
26 
For this example we use a very simple magmatic intrusion as the basis for our 
27 
model. The simulation will be a single event where some molten granite has 
28 
formed a cylindrical dome at the base of some cold sandstone country rock. 
29 
Assuming that the cylinder is very long 
30 
we model a crosssection as shown in \reffig{fig:twodhdmodel}. We will implement 
31 
the same 
32 
diffusion model as we have used for the granite blocks in \refSec{Sec:1DHDv00} 
33 
but will add the second spatial dimension and show how to define 
34 
variables depending on the location of the domain. 
35 
We use \file{onedheatdiff001b.py} as the starting point to develop this model. 
36 

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

41 
To expand upon our 1D problem, the domain must first be expanded. In fact, we 
42 
have solved a two dimensional problem already but essentially ignored the second 
43 
dimension. In our definition phase we create a square domain in $x$ and 
44 
$y$\footnote{In \esc the notation 
45 
$x_{0}$ and $x_{1}$ is used for $x$ and $y$, respectively.} 
46 
that is $600$ meters along each side \reffig{fig:twodhdmodel}. Now we set the 
47 
number of discrete spatial cells to $150$ in both directions and the radius of 
48 
the intrusion to $200$ meters with the centre located at the $300$ meter mark 
49 
on the $x$axis. Thus, the domain variables are; 
50 
\begin{python} 
51 
mx = 600*m #meters  model length 
52 
my = 600*m #meters  model width 
53 
ndx = 150 #mesh steps in x direction 
54 
ndy = 150 #mesh steps in y direction 
55 
r = 200*m #meters  radius of intrusion 
56 
ic = [300*m, 0] #coordinates of the centre of intrusion (meters) 
57 
qH=0.*J/(sec*m**3) #our heat source temperature is zero 
58 
\end{python} 
59 
As before we use 
60 
\begin{python} 
61 
model = Rectangle(l0=mx,l1=my,n0=ndx, n1=ndy) 
62 
\end{python} 
63 
to generate the domain. 
64 

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

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

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

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

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

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

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

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

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

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

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

303 

304 
\section{Advanced Visualisation using VTK} 
305 

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

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

326 
To write the temperatures at each time step to data files in the VTK file format 
327 
one needs to import \verbsaveVTK from the \weipa module and call it: 
328 
\begin{python} 
329 
from esys.weipa import saveVTK 
330 
while t<=tend: 
331 
i+=1 #counter 
332 
t+=h #current time 
333 
mypde.setValue(Y=qH+T*rhocp/h) 
334 
T=mypde.getSolution() 
335 
saveVTK(os.path.join(save_path,"data.%03d.vtu"%i, T=T) 
336 
\end{python} 
337 
The data files, e.g. \file{data.001.vtu}, contain all necessary information to 
338 
visualise the temperature and can be directly processed by \mayavi. Note that 
339 
there is no regridding required. The file extension \file{.vtu} is 
340 
automatically added if not supplied to \verbsaveVTK. 
341 

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

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

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