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 

19 
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 infact warped by the structure of the model and is heavily dependant on the material properties and shape. 
20 

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

27 
We will start by defining our domain and material boundaries using \pycad. This example is located in the \esc directory under \exf and a labeled diagram is available in \reffig{fig:anticlinehrmodel}. As always we start with dependencies. For this example we have: 
28 
\begin{verbatim} 
29 
from esys.pycad import * #domain constructor 
30 
from esys.pycad.gmsh import Design #Finite Element meshing package 
31 
from esys.finley import MakeDomain #Converter for escript 
32 
import os #file path tool 
33 
import numpy as np #numerial python for arrays 
34 
from math import * # math package 
35 
#used to construct polygons for plotting from pycad 
36 
from cblib import getLoopCoords 
37 
\end{verbatim} 
38 
There are two modes available in our example. When \verb modal=1 we indicate to the script that we wish to model an anticline. Othewise when \verb modal=1 we wish to modal a syncline. The modal operator simply changes the phase 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. 
39 

40 
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 one we will use are: 
41 
\begin{verbatim} 
42 
Point() #Create a point in space. 
43 
Line() #Creates a line from a number of points. 
44 
CurveLoop() #Creates a closed loop from a number of lines. 
45 
Surface() #Creates a surface based on a CurveLoop 
46 
\end{verbatim} 
47 
We start by inputting the variables we need to construct the model. 
48 
\begin{verbatim} 
49 
#Model Parameters 
50 
width=5000.0 #width of model 
51 
depth=6000.0 #depth of model 
52 
sspl=51 #number of discrete points in spline 
53 
dsp=width/(sspl1) #dx of spline steps for width 
54 
dep_sp=2500.0 #avg depth of spline 
55 
h_sp=1500.0 #heigh of spline 
56 
orit=1.0 #orientation of spline 1.0=>up 1.0=>down 
57 
\end{verbatim} 
58 
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. 
59 
\begin{verbatim} 
60 
# Overall Domain 
61 
p0=Point(0.0, 0.0, 0.0) 
62 
p1=Point(0.0, depth, 0.0) 
63 
p2=Point(width, depth, 0.0) 
64 
p3=Point(width, 0.0, 0.0) 
65 
\end{verbatim} 
66 
Now lines are defined in an \textbf{anticlockwise} direction using our points. This forms a rectangle around our domain. 
67 
\begin{verbatim} 
68 
l01=Line(p0, p1) 
69 
l12=Line(p1, p2) 
70 
l23=Line(p2, p3) 
71 
l30=Line(p3, p0) 
72 
\end{verbatim} 
73 
These lines form the basis for our domain boundary, which is a closed loop. 
74 
\begin{verbatim} 
75 
c=CurveLoop(l01, l12, l23, l30) 
76 
\end{verbatim} 
77 
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. 
78 
\begin{verbatim} 
79 
# Material Boundary 
80 
x=[ Point(i*dsp\ 
81 
,dep_sp+modal*orit*h_sp*cos(pi*i*dsp/dep_sp+pi))\ 
82 
for i in range(0,sspl)\ 
83 
] 
84 
mysp = Spline(*tuple(x)) #*tuple() forces x to become a tuple 
85 
\end{verbatim} 
86 
The start and end points of the spline can be returned to help define the material boundaries. 
87 
\begin{verbatim} 
88 
x1=Spline.getStartPoint(mysp) 
89 
x2=Spline.getEndPoint(mysp) 
90 
\end{verbatim} 
91 
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. 
92 
\begin{verbatim} 
93 
# TOP BLOCK 
94 
tbl1=Line(p0,x1) 
95 
tbl2=mysp 
96 
tbl3=Line(x2,p3) 
97 
tblockloop = CurveLoop(tbl1,tbl2,tbl3,l30) 
98 
tblock = PlaneSurface(tblockloop) 
99 
\end{verbatim} 
100 
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 retreive its point coordinates with the function \verb getLoopCoords() and then output it with \modnumpy for our solution script. 
101 
\begin{verbatim} 
102 
tpg = getLoopCoords(tblockloop) 
103 
np.savetxt(os.path.join(save_path,"toppg"),tpg,delimiter=" ") 
104 
\end{verbatim} 
105 
We must repeat the above for our every other subdomain. We have only one other, the bottom block. The process is the fairly similar to the top block but with a few differences. The spline points must be reversed by setting the spline as negative. 
106 
\begin{verbatim} 
107 
bbl4=mysp 
108 
\end{verbatim} 
109 
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. 
110 
\begin{verbatim} 
111 
#clockwise check 
112 
bblockloop2=CurveLoop(mysp,Line(x2,p2),Line(p2,p1),Line(p1,x1)) 
113 
bpg = getLoopCoords(bblockloop2) 
114 
np.savetxt(os.path.join(save_path,"botpg"),bpg,delimiter=" ") 
115 
\end{verbatim} 
116 
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. 
117 
To initial 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 simlpe 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 scirpt. We then save the geometry, mesh and create our \esc domain. 
118 
\begin{verbatim} 
119 
# Create a Design which can make the mesh 
120 
d=Design(dim=2, element_size=200) 
121 
# Add the subdomains and flux boundaries. 
122 
d.addItems(PropertySet("top",tblock),PropertySet("bottom",bblock),PropertySet("linebottom",l12)) 
123 
# Create the geometry, mesh and Escript domain 
124 
d.setScriptFileName(os.path.join(save_path,"heatrefraction_mesh001.geo")) 
125 
d.setMeshFileName(os.path.join(save_path,"heatrefraction_mesh001.msh")) 
126 
domain=MakeDomain(d, integrationOrder=1, reducedIntegrationOrder=1, optimizeLabeling=True) 
127 
# Create a file that can be read back in to python with mesh=ReadMesh(fileName) 
128 
domain.write(os.path.join(save_path,"heatrefraction_mesh001.fly")) 
129 
\end{verbatim} 
130 
The creation of our domain and its mesh is complete. Now we must create a solution for our steady state problem. 
131 

132 
\subsection{Steadystate PDE solution} 
133 
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. 
134 
\begin{equation} 
135 
\rho c\hackscore p \frac{\partial T}{\partial t}  \kappa \frac{\partial^{2} T}{\partial x^{2}} = q\hackscore H 
136 
\label{eqn:hd2} 
137 
\end{equation} 
138 
In the steady state PDE however, there is no time dependance. This means that \refEq{eqn:hd2} reduces to 
139 
\begin{equation} 
140 
 \kappa \frac{\partial^{2} T}{\partial x^{2}} = q\hackscore H 
141 
\label{eqn:sshd} 
142 
\end{equation} 
143 
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 
144 
\begin{equation} 
145 
u=r \textit{ where } q>0 
146 
\label{eqn:bdr} 
147 
\end{equation} 
148 
using this functionality forces the solution $u$ to equal the value or $r$ wherever the criterion supplied by $q>0$ is true. 
149 
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. 
150 
\begin{verbatim} 
151 
##ESTABLISHING VARIABLES 
152 
qin=70*Milli*W/(m*m) #our heat source temperature is now zero 
153 
Ti=290.15*K # Kelvin #the starting temperature of our iron bar 
154 
width=5000.0*m 
155 
depth=6000.0*m 
156 
\end{verbatim} 
157 
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. 
158 
\begin{verbatim} 
159 
mymesh=ReadMesh(os.path.join(saved_path,"heatrefraction_mesh001.fly")) 
160 
tpg = np.loadtxt(os.path.join(saved_path,"toppg")) 
161 
tpgx = tpg[:,0] 
162 
tpgy = tpg[:,1] 
163 
bpg = np.loadtxt(os.path.join(saved_path,"botpg")) 
164 
bpgx = bpg[:,0] 
165 
bpgy = bpg[:,1] 
166 

167 
# set up kappa (thermal conductivity across domain) using tags 
168 
kappa=Scalar(0,Function(mymesh)) 
169 
kappa.setTaggedValue("top",2.0) 
170 
kappa.setTaggedValue("bottom",4.0) 
171 
\end{verbatim} 
172 
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() . 
173 
\begin{verbatim} 
174 
#... generate functionspace... 
175 
#... open PDE ... 
176 
mypde=LinearPDE(mymesh) 
177 
#define first coefficient 
178 
mypde.setValue(A=kappa*kronecker(mymesh)) 
179 

180 
# ... set initial temperature .... 
181 
x=mymesh.getX() 
182 

183 
qH=Scalar(0,FunctionOnBoundary(mymesh)) 
184 
qH.setTaggedValue("linebottom",qin) 
185 
mypde.setValue(q=whereZero(x[1]),r=Ti) 
186 
mypde.setValue(y=qH)#,r=17*Celsius) 
187 

188 
# get steady state solution and export to vtk. 
189 
T=mypde.getSolution() 
190 
\end{verbatim} 
191 

192 
\subsection{Quiver plots in \mpl} 
193 
\begin{verbatim} 
194 
# rearrage mymesh to suit solution function space 
195 
oldspacecoords=mymesh.getX() 
196 
coords=Data(oldspacecoords, T.getFunctionSpace()) 
197 

198 
# calculate gradient of solution for quiver plot 
199 
qu=kappa*grad(T) 
200 

201 
# create quiver locations 
202 
quivshape = [20,20] #quivers x and quivers y 
203 
numquiv = quivshape[0]*quivshape[1] # total number of quivers 
204 
maxx = 5000. # maximum x displacement 
205 
dx = maxx/quivshape[0]+1. # quiver x spacing 
206 
maxy = 6000. # maximum y displacement 
207 
dy = maxy/quivshape[1]+1. # quiver y spacing 
208 
qulocs=np.zeros([numquiv,2],float) # memory for quiver locations 
209 
# fill qulocs 
210 
for i in range(0,quivshape[0]1): 
211 
for j in range(0,quivshape[1]1): 
212 
qulocs[i*quivshape[0]+j,:] = [dx+dx*i,dy+dy*j] 
213 
# retreive values for quivers direction and shape from qu 
214 
quL = Locator(qu.getFunctionSpace(),qulocs.tolist()) 
215 
qu = quL(qu) #array of dx,dy for quivers 
216 
qu = np.array(qu) #turn into a numpy array 
217 
qulocs = quL.getX() #returns actual locations from data 
218 
qulocs = np.array(qulocs) #turn into a numpy array 
219 

220 
kappaT = kappa.toListOfTuples(scalarastuple=False) 
221 
coordsK = Data(oldspacecoords, kappa.getFunctionSpace()) 
222 
coordsK = np.array(coordsK.toListOfTuples()) 
223 
coordKX = coordsK[:,0] 
224 
coordKY = coordsK[:,1] 
225 

226 
coords = np.array(coords.toListOfTuples()) 
227 
tempT = T.toListOfTuples(scalarastuple=False) 
228 

229 
coordX = coords[:,0] 
230 
coordY = coords[:,1] 
231 

232 
xi = np.linspace(0.0,5000.0,100) 
233 
yi = np.linspace(6000.0,0.0,100) 
234 
# grid the data. 
235 
zi = pl.matplotlib.mlab.griddata(coordX,coordY,tempT,xi,yi) 
236 
ziK = pl.matplotlib.mlab.griddata(coordKX,coordKY,kappaT,xi,yi) 
237 
# contour the gridded data, plotting dots at the randomly spaced data points. 
238 

239 
pl.matplotlib.pyplot.autumn() 
240 
CKL = pl.fill(tpgx,tpgy,'brown',bpgx,bpgy,'red',zorder=1000) 
241 
#~ CK = pl.contourf(xi,yi,ziK,2) 
242 
CS = pl.contour(xi,yi,zi,5,linewidths=0.5,colors='k') 
243 
pl.clabel(CS, inline=1, fontsize=8) 
244 
pl.title("Heat Refraction across a clinal structure.") 
245 
pl.xlabel("Horizontal Displacement (m)") 
246 
pl.ylabel("Depth (m)") 
247 
#~ CB = pl.colorbar(CS, shrink=0.8, extend='both') 
248 
pl.savefig(os.path.join(saved_path,"heatrefraction001_cont.png")) 
249 

250 
QUIV=pl.quiver(qulocs[:,0],qulocs[:,1],qu[:,0],qu[:,1],angles='xy',color="white") 
251 
pl.title("Heat Refraction across a clinal structure \n with gradient quivers.") 
252 
pl.savefig(os.path.join(saved_path,"heatrefraction001_contqu.png")) 
253 
\end{verbatim} 
254 

255 

256 
\begin{figure} 
257 
\centerline{\includegraphics[width=4.in]{figures/heatrefraction001contqu}} 
258 
\caption{Heat refraction model with gradient indicated by quivers.} 
259 
\label{fig:hr001qumodel} 
260 
\end{figure} 
261 
