1 

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{Steadystate heat refraction} 
15 
Steadystate heat refraction will give us an opportunity to investigate some of the richer features that the \esc package has to offer. One of these is \pycad . The advantage of using \pycad is that it offers an easy method for developing and manipulating complex domains. In conjunction with \gmsh we can develop finite element grids that conform to our domain's shape providing accurate modelling of interfaces and boundaries. Another useful function of \pycad is that we can tag specific areas of our domain with labels as we construct them. These labels can then be used in escript to define properties like material constants and source locations. 
16 

17 
\subsection{Creating the model with \pycad} 
18 
\sslist{heatrefraction_mesher001.py and cblib.py} 
19 

20 
Our first heat refraction model will be a large anticlinal structure that is experiencing a constant temperature at the surface and a steady heat flux at it's base. Our aim is to show that the temperature flux across the surface is not linear from bottom to top but is in fact warped by the structure of the model and is heavily dependant on the material properties and shape. 
21 

22 
\begin{figure}[h!] 
23 
\centerline{\includegraphics[width=4.in]{figures/anticlineheatrefraction}} 
24 
\caption{Heat refraction model with point and line labels.} 
25 
\label{fig:anticlinehrmodel} 
26 
\end{figure} 
27 

28 
We will start by defining our domain and material boundaries using the \pycad module. This example is located in the \esc directory under \exf and a labelled diagram is available in \reffig{fig:anticlinehrmodel}. As always we start with dependencies. For this example we have: 
29 
\begin{verbatim} 
30 
from esys.pycad import * #domain constructor 
31 
from esys.pycad.gmsh import Design #Finite Element meshing package 
32 
from esys.finley import MakeDomain #Converter for escript 
33 
import os #file path tool 
34 
import numpy as np #numerial python for arrays 
35 
from math import * # math package 
36 
#used to construct polygons for plotting from pycad 
37 
from cblib import getLoopCoords 
38 
\end{verbatim} 
39 
There are two modes available in our example. When \verb modal=1 we indicate to the script that we wish to model an anticline. Otherwise when \verb modal=1 we wish to modal a syncline. The modal operator simply changes the orientation of the boundary function so that it is either upwards or downwards curving. A \verb save_path has also been defined so that we can easily separate our data from other examples and our scripts. 
40 

41 
It is now possible to start defining our domain and boundaries. There are a few primary constructors in \pycad that build upon each other to define domains and boundaries; the ones we will use are: 
42 
\begin{verbatim} 
43 
Point() #Create a point in space. 
44 
Line() #Creates a line from a number of points. 
45 
CurveLoop() #Creates a closed loop from a number of lines. 
46 
Surface() #Creates a surface based on a CurveLoop 
47 
\end{verbatim} 
48 
We start by inputting the variables we need to construct the model. 
49 
\begin{verbatim} 
50 
#Model Parameters 
51 
width=5000.0 #width of model 
52 
depth=6000.0 #depth of model 
53 
sspl=51 #number of discrete points in spline 
54 
dsp=width/(sspl1) #dx of spline steps for width 
55 
dep_sp=2500.0 #avg depth of spline 
56 
h_sp=1500.0 #heigh of spline 
57 
orit=1.0 #orientation of spline 1.0=>up 1.0=>down 
58 
\end{verbatim} 
59 
The variables are then used to construct the four corners of our domain, which from the origin has the dimensions of 5000 meters width and 6000 meters depth. This is done with the \verb Point() function which accepts x, y and z coordinates. Our domain is in two dimensions so z should always be zero. Be careful to define all your constructs in an \textbf{anticlockwise} manner otherwise the meshing algorithm may fail. 
60 
\begin{verbatim} 
61 
# Overall Domain 
62 
p0=Point(0.0, 0.0, 0.0) 
63 
p1=Point(0.0, depth, 0.0) 
64 
p2=Point(width, depth, 0.0) 
65 
p3=Point(width, 0.0, 0.0) 
66 
\end{verbatim} 
67 
Now lines are defined in an \textbf{anticlockwise} direction using our points. This forms a rectangle around our domain. 
68 
\begin{verbatim} 
69 
l01=Line(p0, p1) 
70 
l12=Line(p1, p2) 
71 
l23=Line(p2, p3) 
72 
l30=Line(p3, p0) 
73 
\end{verbatim} 
74 
These lines form the basis for our domain boundary, which is a closed loop. 
75 
\begin{verbatim} 
76 
c=CurveLoop(l01, l12, l23, l30) 
77 
\end{verbatim} 
78 
The curve defining our clinal structure is located in approximately the middle of the domain and has a cosinusoidal shape. We define the curve by generating points at discrete intervals; 51 in this case, and then create a smooth curve through the points using the \verb Spline() function. 
79 
\begin{verbatim} 
80 
# Material Boundary 
81 
x=[ Point(i*dsp\ 
82 
,dep_sp+modal*orit*h_sp*cos(pi*i*dsp/dep_sp+pi))\ 
83 
for i in range(0,sspl)\ 
84 
] 
85 
mysp = Spline(*tuple(x)) #*tuple() forces x to become a tuple 
86 
\end{verbatim} 
87 
The start and end points of the spline can be returned to help define the material boundaries. 
88 
\begin{verbatim} 
89 
x1=Spline.getStartPoint(mysp) 
90 
x2=Spline.getEndPoint(mysp) 
91 
\end{verbatim} 
92 
Our top block or material above the clinal/spline boundary is defined in a \textbf{anticlockwise} manner by creating lines and then a closed loop. As we will be meshing the subdomain we also need to assign it a planar surface. 
93 
\begin{verbatim} 
94 
# TOP BLOCK 
95 
tbl1=Line(p0,x1) 
96 
tbl2=mysp 
97 
tbl3=Line(x2,p3) 
98 
tblockloop = CurveLoop(tbl1,tbl2,tbl3,l30) 
99 
tblock = PlaneSurface(tblockloop) 
100 
\end{verbatim} 
101 
It is also necessary to create and export a polygon for \esc so that we can plot the boundaries of our subdomains. First we take the loop from our block and retrieve its point coordinates with the function \verb getLoopCoords() and then output it with \modnumpy for our solution script. 
102 
\begin{verbatim} 
103 
tpg = getLoopCoords(tblockloop) 
104 
np.savetxt(os.path.join(save_path,"toppg"),tpg,delimiter=" ") 
105 
\end{verbatim} 
106 
We must repeat the above for our every other subdomain. We have only one other, the bottom block. The process is fairly similar to the top block but with a few differences. The spline points must be reversed by setting the spline as negative. 
107 
\begin{verbatim} 
108 
bbl4=mysp 
109 
\end{verbatim} 
110 
This reverse spline option unfortunately does not work for the getLoopCoords command, however, the polygon tool will accept clockwise oriented points so we can define a new curve. 
111 
\begin{verbatim} 
112 
#clockwise check 
113 
bblockloop2=CurveLoop(mysp,Line(x2,p2),Line(p2,p1),Line(p1,x1)) 
114 
bpg = getLoopCoords(bblockloop2) 
115 
np.savetxt(os.path.join(save_path,"botpg"),bpg,delimiter=" ") 
116 
\end{verbatim} 
117 
The last few steps in creating our model take the previously defined domain and subdomain points and generate a mesh that can be imported into \esc. 
118 
To initialise the mesh it first needs some design parameters. In this case we have 2 dimensions \verb dim and a specified number of finite elements we wish to apply to the domain \verb element_size . It then becomes a simple task of adding the subdomains and flux boundaries to the design. If we wish to give the elements names as in the example we can use the \verb PropertySet() function. This will allow us to easily define the subdomain properties in the solution script. We then save the geometry, mesh and create our \esc domain. 
119 
\begin{verbatim} 
120 
# Create a Design which can make the mesh 
121 
d=Design(dim=2, element_size=200) 
122 
# Add the subdomains and flux boundaries. 
123 
d.addItems(PropertySet("top",tblock),PropertySet("bottom",bblock),\ 
124 
PropertySet("linebottom",l12)) 
125 
# Create the geometry, mesh and Escript domain 
126 
d.setScriptFileName(os.path.join(save_path,"heatrefraction_mesh001.geo")) 
127 
d.setMeshFileName(os.path.join(save_path,"heatrefraction_mesh001.msh")) 
128 
domain=MakeDomain(d, integrationOrder=1, reducedIntegrationOrder=1,\ 
129 
optimizeLabeling=True) 
130 
# Create a file that can be read back in to python with mesh=ReadMesh(fileName) 
131 
domain.write(os.path.join(save_path,"heatrefraction_mesh001.fly")) 
132 
\end{verbatim} 
133 
The creation of our domain and its mesh is complete. Now we must create a solution for our steady state problem. 
134 

135 
\subsection{Steadystate PDE solution} 
136 
\sslist{heatrefraction.py and cblib.py} 
137 
While a steadystate solution will not give an indication as to the quantitative heat flow in a model, it is useful because it allows us to visualise the direction of heat transport and the equilibrium state of the model. Also, a steadystate model only requires one iteration to find a solution. We know from \refCh{CHAP HEAT DIFF} that the full form of the heat diffusion PDE can be expressed as follows. 
138 
\begin{equation} 
139 
\rho c\hackscore p \frac{\partial T}{\partial t}  \kappa \frac{\partial^{2} T}{\partial x^{2}} = q\hackscore H 
140 
\label{eqn:hd2} 
141 
\end{equation} 
142 
In the steady state PDE however, there is no time dependence. This means that \refEq{eqn:hd2} reduces to 
143 
\begin{equation} 
144 
 \kappa \frac{\partial^{2} T}{\partial x^{2}} = q\hackscore H 
145 
\label{eqn:sshd} 
146 
\end{equation} 
147 
where we see our only variables are $\kappa$ the thermal conductivity and the heat source term $q\hackscore H$. Our boundary conditions are similar to earlier problems with some exceptions. Our heat source $q\hackscore H$ is now located at the base of the model and is implemented using our \modpycad identity \textit{linebottom} and we have a constant temperature along the surface of the model where $z=0$. The \modLPDE module's general form as a provision for these cases of the form 
148 
\begin{equation} 
149 
u=r \textit{ where } q>0 
150 
\label{eqn:bdr} 
151 
\end{equation} 
152 
using this functionality forces the solution $u$ to equal the value or $r$ wherever the criterion supplied by $q>0$ is true. 
153 
The structure of our solution script is similar to those we have used already. Firstly, the dependencies must be imported. Then the variables must be defined. In this case there are four major variables. \verb qin is the temperature of our source, \verb Ti is the temperature on the surface at $z=0$ and our width and depth of the model. 
154 
\begin{verbatim} 
155 
##ESTABLISHING VARIABLES 
156 
qin=70*Milli*W/(m*m) #our heat source temperature is now zero 
157 
Ti=290.15*K # Kelvin #the starting temperature of our iron bar 
158 
width=5000.0*m 
159 
depth=6000.0*m 
160 
\end{verbatim} 
161 
The mesh is now imported via the \verb ReadMesh() function and loaded into a domain variable \verb mymesh . The structural polygons are next. With the mesh imported it is now possible to use our tagging property to set up our PDE coefficients. In this case $\kappa$ is set via the \verb setTaggedValue() function which takes two arguments, the name of the tagged points and the value to assign to them. 
162 
\begin{verbatim} 
163 
mymesh=ReadMesh(os.path.join(saved_path,"heatrefraction_mesh001.fly")) 
164 
tpg = np.loadtxt(os.path.join(saved_path,"toppg")) 
165 
tpgx = tpg[:,0] 
166 
tpgy = tpg[:,1] 
167 
bpg = np.loadtxt(os.path.join(saved_path,"botpg")) 
168 
bpgx = bpg[:,0] 
169 
bpgy = bpg[:,1] 
170 

171 
# set up kappa (thermal conductivity across domain) using tags 
172 
kappa=Scalar(0,Function(mymesh)) 
173 
kappa.setTaggedValue("top",2.0) 
174 
kappa.setTaggedValue("bottom",4.0) 
175 
\end{verbatim} 
176 
Setting up the PDE is a little more complicated. The linear PDE \verb mypde is again assigned the domain \verb mymesh and \verb A is set. \verb qH is set on the boundary by creating \verb qH as a function space of type \verb FunctionOnBoundary which is in turn based on \verb mymesh . We can then set the tagged value \textit{linebottom} to the condition \verb qin . We also need to use our condition from \refEq{eqn:bdr}. By extracting \verb x we can set the domain at zero to 1 and everywhere else 0 via the \verb whereZero() function. This is put into \verb q and \verb r is set to \verb Ti our temperature at the surface. Note \verb y=qH . As there is no time derivative in the equation we do not have any iterative steps and there is one solution to the problem. This is found in the step \verb T=mypde.getSolution() . 
177 
\begin{verbatim} 
178 
#... generate functionspace... 
179 
#... open PDE ... 
180 
mypde=LinearPDE(mymesh) 
181 
#define first coefficient 
182 
mypde.setValue(A=kappa*kronecker(mymesh)) 
183 

184 
# ... set initial temperature .... 
185 
x=mymesh.getX() 
186 

187 
qH=Scalar(0,FunctionOnBoundary(mymesh)) 
188 
qH.setTaggedValue("linebottom",qin) 
189 
mypde.setValue(q=whereZero(x[1]),r=Ti) 
190 
mypde.setValue(y=qH) 
191 

192 
# get steady state solution. 
193 
T=mypde.getSolution() 
194 
\end{verbatim} 
195 
The problem is now solved and plotting is required to visualise the data. 
196 

197 
\subsection{Line profiles of 2D data} 
198 
It is possible to take slices or profile sections of 2d data using \mpl. To acheive this there are three major steps, firstly the solution must be smoothed across the domain using the \verb projector() function. 
199 
\begin{verbatim} 
200 
proj=Projector(mymesh) 
201 
smthT=proj(T) 
202 
\end{verbatim} 
203 
The data must then be moved to a regular grid. 
204 
\begin{verbatim} 
205 
xi,yi,zi = toRegGrid(smthT,mymesh,200,200,width,depth) 
206 
\end{verbatim} 
207 
The profile can then be plotted by taking a slice of the data output \verb zi . 
208 
\begin{verbatim} 
209 
cut=int(len(xi)/2) 
210 
pl.clf() 
211 
pl.plot(zi[:,cut],yi) 
212 
pl.title("Heat Refraction Temperature Depth Profile") 
213 
pl.xlabel("Temperature (K)") 
214 
pl.ylabel("Depth (m)") 
215 
\end{verbatim} 
216 
This process can be repeated for various variations of the solultion. In this case we have temperature, temperature gradient, thermal conductivity and heat flow \reffig{figs:dps}. 
217 
\begin{figure} 
218 
\centering 
219 
\subfigure[Temperature Depth Profile]{\includegraphics[width=3in]{figures/heatrefraction001tdp.png}\label{fig:tdp}} 
220 
\subfigure[Temperature Gradient Depth Profile]{\includegraphics[width=3in]{figures/heatrefraction001tgdp.png}\label{fig:tgdp}} 
221 
\subfigure[Thermal Conductivity Profile]{\includegraphics[width=3in]{figures/heatrefraction001tcdp.png}\label{fig:tcdp}} 
222 
\subfigure[Heat Flow Depth Profile]{\includegraphics[width=3in]{figures/heatrefraction001hf.png}\label{fig:hf}} 
223 
\caption{Depth profiles down centre of model.} 
224 
\label{figs:dps} 
225 
\end{figure} 
226 

227 

228 

229 

230 
\subsection{Quiver plots in \mpl} 
231 
To visualise this data we are going to do three things. Firstly we will generate the gradient of the solution and create some quivers or arrows that describe the gradient. Then we need to generate some temperatur contours and thirdly, we will create some polygons from our meshing output files to show our material boundaries. 
232 

233 
The first step is very easily provided by \esc, where we can use the \verb grad() function to generate our gradient. It is then necessary to mold the solution to the domain node points. We do this by first extracting the nodes to \verb oldspacecoords and transforming them via interpolation to suit our solution which is in the function space form of \verb T . 
234 
\begin{verbatim} 
235 
# calculate gradient of solution for quiver plot 
236 
qu=kappa*grad(T) 
237 

238 
# rearrage mymesh to suit solution function space 
239 
oldspacecoords=mymesh.getX() 
240 
coords=Data(oldspacecoords, T.getFunctionSpace()) 
241 
\end{verbatim} 
242 
Now the number of quivers we want is specified as $20 \times 20 = 400$ and their locations calculated based on our coordinates in \verb qu . All our date is then transformed to tuples for \modmpl . 
243 
\begin{verbatim} 
244 
quivshape = [20,20] #quivers x and quivers y 
245 
# function to calculate quiver locations 
246 
qu,qulocs = toQuivLocs(quivshape,width,depth,qu) 
247 

248 
kappaT = kappa.toListOfTuples(scalarastuple=False) 
249 
coordsK = Data(oldspacecoords, kappa.getFunctionSpace()) 
250 
coordKX, coordKY = toXYTuple(coordsK) 
251 

252 
tempT = T.toListOfTuples(scalarastuple=False) 
253 
coordX, coordY = toXYTuple(coords) 
254 
\end{verbatim} 
255 
In a similar fashion to the contouring, our data needs to be transferred to a regular grid to contour the temperature. Our polygons are handled by the \verb fill() function making sure they are ordered to the lowest level of the plot beneath all other data. Labeling is then added and finally the Quivers using the \verb quiver() function. More detail on each \modmpl function is available from the \modmpl website. 
256 
\begin{verbatim} 
257 
xi = np.linspace(0.0,width,100) 
258 
yi = np.linspace(depth,0.0,100) 
259 
# grid the data. 
260 
zi = pl.matplotlib.mlab.griddata(coordX,coordY,tempT,xi,yi) 
261 
ziK = pl.matplotlib.mlab.griddata(coordKX,coordKY,kappaT,xi,yi) 
262 
# contour the gridded data, plotting dots at the randomly spaced data points. 
263 
pl.matplotlib.pyplot.autumn() 
264 

265 
CKL = pl.fill(tpgx,tpgy,'brown',bpgx,bpgy,'red',zorder=1000) 
266 
#~ CK = pl.contourf(xi,yi,ziK,2) 
267 
CS = pl.contour(xi,yi,zi,5,linewidths=0.5,colors='k') 
268 

269 
pl.clabel(CS, inline=1, fontsize=8) 
270 
pl.title("Heat Refraction across a clinal structure.") 
271 
pl.xlabel("Horizontal Displacement (m)") 
272 
pl.ylabel("Depth (m)") 
273 
#~ CB = pl.colorbar(CS, shrink=0.8, extend='both') 
274 
pl.savefig(os.path.join(saved_path,"heatrefraction001_cont.png")) 
275 

276 
QUIV=pl.quiver(qulocs[:,0],qulocs[:,1],qu[:,0],qu[:,1],angles='xy',color="white") 
277 
pl.title("Heat Refraction across a clinal structure \n with gradient quivers.") 
278 
pl.savefig(os.path.join(saved_path,"heatrefraction001_contqu.png")) 
279 
\end{verbatim} 
280 
The data has now been plotted. 
281 
\begin{figure}[ht] 
282 
\centerline{\includegraphics[width=4.in]{figures/heatrefraction001contqu}} 
283 
\caption{Heat refraction model with gradient indicated by quivers.} 
284 
\label{fig:hr001qumodel} 
285 
\end{figure} 
286 

287 
\newpage 
288 
\subsection{Fault and Overburden Model} 
289 
A slightly more complicated model can be found in the examples \textit{heatrefraction_mesher002.py} and \textit{heatrefractoin002.py} where three blocks are used within the model. 
290 
\begin{figure}[ht] 
291 
\centerline{\includegraphics[width=4.in]{figures/heatrefraction002contqu}} 
292 
\caption{Heat refraction model with three blocks and gradient quivers.} 
293 
\end{figure} 
294 
