1 
ahallam 
2401 

2 


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
3 


% 
4 


% Copyright (c) 20032009 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 


\section{Two Dimensional Heat Diffusion for a basic Magmatic Intrusion} 
15 
ahallam 
2658 
\sslist{twodheatdiff001.py and cblib.py} 
16 
gross 
2878 
\label{Sec:2DHD} 
17 
ahallam 
2801 
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 hemispherical 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. 
18 
ahallam 
2401 

19 
ahallam 
2606 
\begin{figure}[h!] 
20 


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


\caption{2D model: granitic intrusion of sandstone country rock.} 
22 


\label{fig:twodhdmodel} 
23 


\end{figure} 
24 



25 
ahallam 
2801 
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 in either direction. The radius of the intrusion will be 200 meters with the centre located at the 300 meter mark on the xaxis. The domain variables are; 
26 


\begin{python} 
27 
ahallam 
2645 
mx = 600*m #meters  model length 
28 


my = 600*m #meters  model width 
29 


ndx = 100 #mesh steps in x direction 
30 


ndy = 100 #mesh steps in y direction 
31 


r = 200*m #meters  radius of intrusion 
32 


ic = [300, 0] #centre of intrusion (meters) 
33 


q=0.*Celsius #our heat source temperature is now zero 
34 
ahallam 
2801 
\end{python} 
35 


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 in Section \ref{sec:key}. Prior to setting up the PDE the boundary between the two materials must be established. The distance $s$ between two points in Cartesian coordinates is defined as: 
36 
ahallam 
2401 
\begin{equation} 
37 


(x_{1}x_{0})^{2}+(y_{1}y_{0})^{2} = s^{2} 
38 


\end{equation} 
39 
ahallam 
2801 
If we define the point $[x_{0},y_{0}]$ as $c$ 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 $c$. All the 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: 
40 


\begin{python} 
41 
ahallam 
2401 
bound = length(xic)r #where the boundary will be located 
42 
ahallam 
2801 
\end{python} 
43 


This definition of the boundary can now be used with the \verb wherePositive() and \verb whereNegative() commands 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: 
44 


\begin{python} 
45 
ahallam 
2401 
A = (kappai)*whereNegative(bound)+(kappac)*wherePositive(bound) 
46 


D = (rhocpi/h)*whereNegative(bound)+(rhocpc/h)*wherePositive(bound) 
47 



48 


mypde.setValue(A=A*kronecker(model),D=D,d=eta,y=eta*Tc) 
49 
ahallam 
2801 
\end{python} 
50 
ahallam 
2401 
Our PDE has now been properly established. The initial conditions for temperature are set out in a similar matter: 
51 
ahallam 
2801 
\begin{python} 
52 
ahallam 
2658 
#defining the initial temperatures. 
53 


T= Ti*whereNegative(bound)+Tc*wherePositive(bound) 
54 
ahallam 
2801 
\end{python} 
55 
ahallam 
2645 
The iteration process now begins as before, but using our new conditions for \verb D as defined above. 
56 
ahallam 
2606 

57 


\subsection{Contouring escript data} 
58 
ahallam 
2801 
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 \verb numpy array known as a \verb tuple . Ths multidimensional array is then broken down into individual $x$ and $y$ arrays. This is a one step process using the function \verb toXYTuple . We also need to generate our regular grid which is done using the \modnumpy function \verb linspace . 
59 


\begin{python} 
60 
ahallam 
2645 
# rearrage mymesh to suit solution function space for contouring 
61 
ahallam 
2606 
oldspacecoords=model.getX() 
62 


coords=Data(oldspacecoords, T.getFunctionSpace()) 
63 
ahallam 
2645 
coordX, coordY = toXYTuple(coords) 
64 
ahallam 
2606 
# create regular grid 
65 
ahallam 
2645 
xi = np.linspace(0.0,mx,100) 
66 


yi = np.linspace(0.0,my,100) 
67 
ahallam 
2801 
\end{python} 
68 


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 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 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. 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. 
69 
ahallam 
2606 

70 
ahallam 
2801 
\begin{python} 
71 
ahallam 
2645 
#grid the data. 
72 


zi = pl.matplotlib.mlab.griddata(coordX,coordY,tempT,xi,yi) 
73 


# contour the gridded data, plotting dots at the randomly spaced data points. 
74 


pl.matplotlib.pyplot.autumn() 
75 


pl.contourf(xi,yi,zi,10) 
76 


CS = pl.contour(xi,yi,zi,5,linewidths=0.5,colors='k') 
77 


pl.clabel(CS, inline=1, fontsize=8) 
78 


pl.axis([0,600,0,600]) 
79 


pl.title("Heat diffusion from an intrusion.") 
80 


pl.xlabel("Horizontal Displacement (m)") 
81 


pl.ylabel("Depth (m)") 
82 


pl.savefig(os.path.join(save_path,"heatrefraction%03d.png") %i) 
83 


pl.clf() 
84 
ahallam 
2801 
\end{python} 
85 
ahallam 
2645 

86 


\begin{figure}[h!] 
87 


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


\caption{2D model: Total temperature distribution ($T$) at time $t=50$.} 
89 
ahallam 
2681 
\label{fig:twodhdans} 
90 
gross 
2878 
\end{figure} 
91 



92 


\subsection{Advanced Visualization using VTK} 
93 



94 


\subsubsection{Parallel scripts (MPI)} 
95 


In some of the example files for this cookbook the plotting commands are a little different. 
96 


For example, 
97 


\begin{python} 
98 


if getMPIRankWorld() == 0: 
99 


pl.savefig(os.path.join(save_path+"/totT","ttrodpyplot%03d.png")%i) 
100 


pl.clf() 
101 


\end{python} 
102 



103 


The additional \verb if statement is not necessary for normal desktop use. 
104 


It becomes important for scripts run on parallel computers. 
105 


Its purpose is to ensure that only one copy of the file is written. 
106 


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