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 
\section{Steadystate Heat Refraction} 
15 
\label{STEADYSTATE HEAT REFRACTION} 
16 

17 

18 
New concepts will include nonlinear boundaries and the ability to prescribe location specific variables. 
19 

20 
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 models. In conjunction with \gmsh we can develop finite element grids that conform to our domain's shape providing accurate modeling 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 \esc to define properties like material constants and source locations. 
21 

22 
\subsection{Creating the model with \pycad} 
23 
\sslist{heatrefraction_mesher001.py and cblib.py} 
24 

25 
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 upon the material properties and shape. 
26 

27 
\begin{figure}[h!] 
28 
\centerline{\includegraphics[width=4.in]{figures/anticlineheatrefraction}} 
29 
\caption{Heat refraction model with point and line labels.} 
30 
\label{fig:anticlinehrmodel} 
31 
\end{figure} 
32 

33 
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}. The dependencies come first and for this example we have: 
34 
\begin{python} 
35 
from esys.pycad import * #domain constructor 
36 
from esys.pycad.gmsh import Design #Finite Element meshing package 
37 
from esys.finley import MakeDomain #Converter for escript 
38 
import os #file path tool 
39 
import numpy as np #numerial python for arrays 
40 
from math import * # math package 
41 
#used to construct polygons for plotting from pycad 
42 
from cblib import getLoopCoords 
43 
\end{python} 
44 
There are two modes available in our example. When \verb modal=1 this indicates to the script that we wish to model an anticline. Otherwise when \verb modal=1 this will model 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. 
45 

46 
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: 
47 
\begin{python} 
48 
Point() #Create a point in space. 
49 
Line() #Creates a line from a number of points. 
50 
CurveLoop() #Creates a closed loop from a number of lines. 
51 
Surface() #Creates a surface based on a CurveLoop 
52 
\end{python} 
53 
We start by inputting the variables we need to construct the model. 
54 
\begin{python} 
55 
#Model Parameters 
56 
width=5000.0 #width of model 
57 
depth=6000.0 #depth of model 
58 
sspl=51 #number of discrete points in spline 
59 
dsp=width/(sspl1) #dx of spline steps for width 
60 
dep_sp=2500.0 #avg depth of spline 
61 
h_sp=1500.0 #heigh of spline 
62 
orit=1.0 #orientation of spline 1.0=>up 1.0=>down 
63 
\end{python} 
64 
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. 
65 
\begin{python} 
66 
# Overall Domain 
67 
p0=Point(0.0, 0.0, 0.0) 
68 
p1=Point(0.0, depth, 0.0) 
69 
p2=Point(width, depth, 0.0) 
70 
p3=Point(width, 0.0, 0.0) 
71 
\end{python} 
72 
Now lines are defined in an \textbf{anticlockwise} direction using our points. This forms a rectangle around our domain. 
73 
\begin{python} 
74 
l01=Line(p0, p1) 
75 
l12=Line(p1, p2) 
76 
l23=Line(p2, p3) 
77 
l30=Line(p3, p0) 
78 
\end{python} 
79 
These lines form the basis for our domain boundary, which is a closed loop. 
80 
\begin{python} 
81 
c=CurveLoop(l01, l12, l23, l30) 
82 
\end{python} 
83 
The curve defining our clinal structure is located in approximately the middle of the domain and has a sinusoidal 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. 
84 
\begin{python} 
85 
# Material Boundary 
86 
x=[ Point(i*dsp\ 
87 
,dep_sp+modal*orit*h_sp*cos(pi*i*dsp/dep_sp+pi))\ 
88 
for i in range(0,sspl)\ 
89 
] 
90 
mysp = Spline(*tuple(x)) #*tuple() forces x to become a tuple 
91 
\end{python} 
92 
The start and end points of the spline can be returned to help define the material boundaries. 
93 
\begin{python} 
94 
x1=Spline.getStartPoint(mysp) 
95 
x2=Spline.getEndPoint(mysp) 
96 
\end{python} 
97 
The 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. 
98 
\begin{python} 
99 
# TOP BLOCK 
100 
tbl1=Line(p0,x1) 
101 
tbl2=mysp 
102 
tbl3=Line(x2,p3) 
103 
tblockloop = CurveLoop(tbl1,tbl2,tbl3,l30) 
104 
tblock = PlaneSurface(tblockloop) 
105 
\end{python} 
106 
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. 
107 
\begin{python} 
108 
tpg = getLoopCoords(tblockloop) 
109 
np.savetxt(os.path.join(save_path,"toppg"),tpg,delimiter=" ") 
110 
\end{python} 
111 
We must repeat the above for every other subdomain. In this example there is 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. 
112 
\begin{python} 
113 
bbl4=mysp 
114 
\end{python} 
115 
This reverse spline option unfortunately does not work for the getLoopCoords command, however, the \modmpl polygon tool will accept clockwise oriented points so we can define a new curve. 
116 
\begin{python} 
117 
#clockwise check 
118 
bblockloop2=CurveLoop(mysp,Line(x2,p2),Line(p2,p1),Line(p1,x1)) 
119 
bpg = getLoopCoords(bblockloop2) 
120 
np.savetxt(os.path.join(save_path,"botpg"),bpg,delimiter=" ") 
121 
\end{python} 
122 
The last few steps in creating the model take the previously defined domain and subdomain points and generate a mesh that can be imported into \esc. 
123 
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 that need to be applied to the domain \verb element_size . It then becomes a simple task of adding the subdomains and flux boundaries to the design. Each element of our model can be given an identifier which makes it easier to define the subdomain properties in the solution script. This is done using the \verb PropertySet() function. The geometry and mesh are then saved so the \esc domain can be created. 
124 
\begin{python} 
125 
# Create a Design which can make the mesh 
126 
d=Design(dim=2, element_size=200) 
127 
# Add the subdomains and flux boundaries. 
128 
d.addItems(PropertySet("top",tblock),PropertySet("bottom",bblock),\ 
129 
PropertySet("linebottom",l12)) 
130 
# Create the geometry, mesh and \esc domain 
131 
d.setScriptFileName(os.path.join(save_path,"heatrefraction_mesh001.geo")) 
132 
d.setMeshFileName(os.path.join(save_path,"heatrefraction_mesh001.msh")) 
133 
domain=MakeDomain(d, integrationOrder=1, reducedIntegrationOrder=1,\ 
134 
optimizeLabeling=True) 
135 
# Create a file that can be read back in to python with mesh=ReadMesh(fileName) 
136 
domain.write(os.path.join(save_path,"heatrefraction_mesh001.fly")) 
137 
\end{python} 
138 
The creation of our domain and its mesh is complete. Now we must create a solution for our steady state problem. 
139 

140 
\subsection{Steadystate PDE solution} 
141 
\sslist{heatrefraction.py and cblib.py} 
142 
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. 
143 
\begin{equation} 
144 
\rho c\hackscore p \frac{\partial T}{\partial t}  \kappa \frac{\partial^{2} T}{\partial x^{2}} = q\hackscore H 
145 
\label{eqn:hd2} 
146 
\end{equation} 
147 
In the steady state PDE however, there is no time dependence. This means that \refEq{eqn:hd2} reduces to 
148 
\begin{equation} 
149 
 \kappa \frac{\partial^{2} T}{\partial x^{2}} = q\hackscore H 
150 
\label{eqn:sshd} 
151 
\end{equation} 
152 
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 
153 
\begin{equation} 
154 
u=r \textit{ where } q>0 
155 
\label{eqn:bdr} 
156 
\end{equation} 
157 
using this functionality forces the solution $u$ to equal the value or $r$ wherever the criterion supplied by $q>0$ is true. 
158 
The structure of our solution script is similar to those we have used already. We begin by importing the dependencies and defining the PDE and script variables. 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. 
159 
\begin{python} 
160 
##ESTABLISHING VARIABLES 
161 
qin=70*Milli*W/(m*m) #our heat source temperature is now zero 
162 
Ti=290.15*K # Kelvin #the starting temperature of our iron bar 
163 
width=5000.0*m 
164 
depth=6000.0*m 
165 
\end{python} 
166 
The mesh is now imported via the \verb ReadMesh() function and loaded into a domain variable \verb mymesh . This is followed by the structural polygons. 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. 
167 
\begin{python} 
168 
mymesh=ReadMesh(os.path.join(saved_path,"heatrefraction_mesh001.fly")) 
169 
tpg = np.loadtxt(os.path.join(saved_path,"toppg")) 
170 
tpgx = tpg[:,0] 
171 
tpgy = tpg[:,1] 
172 
bpg = np.loadtxt(os.path.join(saved_path,"botpg")) 
173 
bpgx = bpg[:,0] 
174 
bpgy = bpg[:,1] 
175 

176 
# set up kappa (thermal conductivity across domain) using tags 
177 
kappa=Scalar(0,Function(mymesh)) 
178 
kappa.setTaggedValue("top",2.0) 
179 
kappa.setTaggedValue("bottom",4.0) 
180 
\end{python} 
181 
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() . 
182 
\begin{python} 
183 
#... generate functionspace... 
184 
#... open PDE ... 
185 
mypde=LinearPDE(mymesh) 
186 
#define first coefficient 
187 
mypde.setValue(A=kappa*kronecker(mymesh)) 
188 

189 
# ... set initial temperature .... 
190 
x=mymesh.getX() 
191 

192 
qH=Scalar(0,FunctionOnBoundary(mymesh)) 
193 
qH.setTaggedValue("linebottom",qin) 
194 
mypde.setValue(q=whereZero(x[1]),r=Ti) 
195 
mypde.setValue(y=qH) 
196 

197 
# get steady state solution. 
198 
T=mypde.getSolution() 
199 
\end{python} 
200 
The problem is now solved and plotting is required to visualise the data. 
201 

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

232 

233 

234 

235 
\subsection{Quiver plots in \mpl} 
236 
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 temperature contours and thirdly, we will create some polygons from our meshing output files to show our material boundaries. 
237 

238 
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 mould 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 . 
239 
\begin{python} 
240 
# calculate gradient of solution for quiver plot 
241 
qu=kappa*grad(T) 
242 

243 
# rearrange mymesh to suit solution function space 
244 
oldspacecoords=mymesh.getX() 
245 
coords=Data(oldspacecoords, T.getFunctionSpace()) 
246 
\end{python} 
247 
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 data is then transformed to tuples for \modmpl . 
248 
\begin{python} 
249 
quivshape = [20,20] #quivers x and quivers y 
250 
# function to calculate quiver locations 
251 
qu,qulocs = toQuivLocs(quivshape,width,depth,qu) 
252 

253 
kappaT = kappa.toListOfTuples(scalarastuple=False) 
254 
coordsK = Data(oldspacecoords, kappa.getFunctionSpace()) 
255 
coordKX, coordKY = toXYTuple(coordsK) 
256 

257 
tempT = T.toListOfTuples(scalarastuple=False) 
258 
coordX, coordY = toXYTuple(coords) 
259 
\end{python} 
260 
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. Labelling is then added and finally the Quivers using the \verb quiver() function. More detail on each \modmpl function is available from the \modmpl website. 
261 
\begin{python} 
262 
xi = np.linspace(0.0,width,100) 
263 
yi = np.linspace(depth,0.0,100) 
264 
# grid the data. 
265 
zi = pl.matplotlib.mlab.griddata(coordX,coordY,tempT,xi,yi) 
266 
ziK = pl.matplotlib.mlab.griddata(coordKX,coordKY,kappaT,xi,yi) 
267 
# contour the gridded data, plotting dots at the randomly spaced data points. 
268 
pl.matplotlib.pyplot.autumn() 
269 

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

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

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

292 
\newpage 
293 
\subsection{Fault and Overburden Model} 
294 
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. 
295 
\begin{figure}[ht] 
296 
\centerline{\includegraphics[width=4.in]{figures/heatrefraction002contqu}} 
297 
\caption{Heat refraction model with three blocks and gradient quivers.} 
298 
\end{figure} 
299 
