1 

2 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
3 
% Copyright (c) 20032012 by University of Queensland 
4 
% http://www.uq.edu.au 
5 
% 
6 
% Primary Business: Queensland, Australia 
7 
% Licensed under the Open Software License version 3.0 
8 
% http://www.opensource.org/licenses/osl3.0.php 
9 
% 
10 
% Development until 2012 by Earth Systems Science Computational Center (ESSCC) 
11 
% Development since 2012 by School of Earth Sciences 
12 
% 
13 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
14 

15 
\section{Example 5: A Heat Refraction Model} 
16 
\label{example5} 
17 
\sslist{example05a.py and cblib.py} 
18 
Our heat refraction model will be a large anticlinal structure that is subject 
19 
to a constant temperature at the surface and experiencing a steady heat flux at 
20 
its base. Our aim is to show that the temperature flux across the surface is 
21 
not linear from bottom to top, but is in fact warped by the structure of the 
22 
model. The heat flow pattern demonstrates the dependence upon the material 
23 
properties and the shape of the interface. 
24 

25 
The script of \refSec{example4} is modified by subdividing the block into two 
26 
parts. The curve separating the two blocks is given as a spline, see 
27 
\reffig{fig:anticlinehrmodel}. The data points 
28 
used to define the curve may be imported from a database of measurements 
29 
(\textit{e.g.} borehole depth data), but for simplicity it is assumed here that 
30 
the coordinates are known in an analytic form. 
31 

32 
\begin{figure}[ht] 
33 
\centerline{\includegraphics[width=4.in]{figures/anticlineheatrefraction}} 
34 
\caption{Example 5a: Heat refraction model with point and line labels} 
35 
\label{fig:anticlinehrmodel} 
36 
\end{figure} 
37 

38 
There are two modes available in this example. When \verbmodal=1, this 
39 
indicates to the script that the model should be an anticline. Otherwise, when 
40 
\verbmodal=1, the model is a syncline. The modal operator simply changes the 
41 
orientation of the boundary function so that it is either upwards or downwards 
42 
curving. A \verbsave_path has also been defined so that we can easily separate 
43 
our data from other examples and our scripts. 
44 

45 
It is now possible to start defining our domain and boundaries. 
46 

47 
The curve defining our clinal structure is located approximately in the middle 
48 
of the domain and has a sinusoidal shape. We define the curve by generating 
49 
points at discrete intervals; $51$ in this case, and then create a smooth curve 
50 
through the points using the \verbSpline() function. 
51 
\begin{python} 
52 
# Material Boundary 
53 
x=[ Point(i*dsp\ 
54 
,dep_sp+modal*orit*h_sp*cos(pi*i*dsp/dep_sp+pi))\ 
55 
for i in range(0,sspl)\ 
56 
] 
57 
mysp = Spline(*tuple(x)) #*tuple() forces x to become a tuple 
58 
\end{python} 
59 
The start and end points of the spline can be returned to help define the 
60 
material boundaries. 
61 
\begin{python} 
62 
x1=mysp.getStartPoint() 
63 
x2=mysp.getEndPoint() 
64 
\end{python} 
65 
The top block or material above the clinal/spline boundary is defined in an 
66 
\textbf{anticlockwise} manner by creating lines and then a closed loop. By 
67 
meshing the subdomain we also need to assign it a planar surface. 
68 
\begin{python} 
69 
# TOP BLOCK 
70 
tbl1=Line(p0,x1) 
71 
tbl2=mysp 
72 
tbl3=Line(x2,p3) 
73 
l30=Line(p3, p0) 
74 
tblockloop = CurveLoop(tbl1,tbl2,tbl3,l30) 
75 
tblock = PlaneSurface(tblockloop) 
76 
\end{python} 
77 
This process is repeated for every other subdomain. In this example there is 
78 
only one other, the bottom block. The process is similar to the top block but 
79 
with a few differences. The spline points must be reversed by setting the spline 
80 
as negative. 
81 
\begin{python} 
82 
bbl4=mysp 
83 
\end{python} 
84 
This reverse spline option unfortunately does not work for the 
85 
\verbgetLoopCoords command, however, the \modmpl polygon tool will accept 
86 
clockwise oriented points so we can define a new curve. 
87 
\begin{python} 
88 
#clockwise check 
89 
bblockloop=CurveLoop(mysp,Line(x2,p2),Line(p2,p1),Line(p1,x1)) 
90 
\end{python} 
91 
The last few steps in creating the domain require that the previously defined 
92 
domain and subdomain points are submitted to generate a mesh that can be 
93 
imported into \esc. 
94 
To initialise the mesh it first needs some design parameters. In this case we 
95 
have 2 dimensions \verbdim and a specified number of finite elements that need 
96 
to be applied to the domain \verbelement_size. It then becomes a simple task 
97 
of adding the subdomains and flux boundaries to the design. Each element of our 
98 
model can be given an identifier which makes it easier to define the subdomain 
99 
properties in the solution script. This is done using the 
100 
\verbPropertySet() function. The geometry and mesh are then saved so the 
101 
\esc domain can be created. 
102 
\begin{python} 
103 
# Create a Design which can make the mesh 
104 
d=Design(dim=2, element_size=200) 
105 
# Add the subdomains and flux boundaries. 
106 
d.addItems(PropertySet("top",tblock),PropertySet("bottom",bblock),\ 
107 
PropertySet("linebottom",l12)) 
108 
# Create the geometry, mesh and escript domain 
109 
d.setScriptFileName(os.path.join(save_path,"example05.geo")) 
110 
d.setMeshFileName(os.path.join(save_path,"example05.msh")) 
111 
domain=MakeDomain(d,optimizeLabeling=True) 
112 
\end{python} 
113 
The creation of our domain and its mesh is now complete. 
114 

115 
With the mesh imported it is now possible to use our tagging property to set up 
116 
our PDE coefficients. In this case $\kappa$ is set via the 
117 
\verbsetTaggedValue() function which takes two arguments, the name of the 
118 
tagged points and the value to assign to them. 
119 
\begin{python} 
120 
# set up kappa (thermal conductivity across domain) using tags 
121 
kappa=Scalar(0,Function(domain)) 
122 
kappa.setTaggedValue("top",2.0*W/m/K) 
123 
kappa.setTaggedValue("bottom",4.0*W/m/K) 
124 
\end{python} 
125 
No further changes are required to the PDE solution step, see 
126 
\reffig{fig:anticlinetemp} for the result. 
127 

128 
\begin{figure}[ht] 
129 
\centerline{\includegraphics[width=4.in]{figures/heatrefraction}} 
130 
\caption{Example 5a: Temperature Distribution in the Heat Refraction Model} 
131 
\label{fig:anticlinetemp} 
132 
\end{figure} 
133 
\clearpage 
134 

135 
\section{Line Profiles of 2D Data} 
136 
\sslist{example05b.py and cblib.py} 
137 
We want to investigate the profile of the data of the last example. 
138 
Of particular interest is the depth profile of the heat flux which is the 
139 
second component of $\kappa \nabla T$. The script from the previous section 
140 
is extended to show how a vertical profile can be plotted. 
141 

142 
The first important piece of information, is that \esc assumes that $\kappa 
143 
\nabla T$ is not smooth and that the point values of this solution are defined 
144 
at numerical interpolation points. This assumption is reasonable as 
145 
the flux is the product of the piecewise constant function $\kappa$ and 
146 
the gradient of the temperature $T$ which has a discontinuity at the rock 
147 
interface. 
148 
Before plotting this function we need to smooth the solution using the 
149 
\verbProjector() class; 
150 
\begin{python} 
151 
from esys.escript.pdetools import Projector 
152 
proj=Projector(domain) 
153 
qu=proj(kappa*grad(T)) 
154 
\end{python} 
155 
The \verbproj object provides a mechanism to distribute values given at the 
156 
numerical interpolation points to the nodes of the FEM mesh  the heat flux in 
157 
this example. \verbqu has the same function space as the temperature 
158 
\verbT. The smoothed flux is interpolated to a regular $200\times 200$ grid 
159 
via; 
160 
\begin{python} 
161 
xiq,yiq,ziq = toRegGrid(qu[1],200,200) 
162 
\end{python} 
163 
using the \verbtoRegGrid function from the cookbook library which we are 
164 
using for the contour plot. 
165 
At return \verbziq[j,i] is the value of vertical heat flux at point 
166 
(\verbxiq[i],\verbyiq[j]). We can easily create deep profiles now by 
167 
plotting slices \verbziq[:,i] over \verbyiq. The following script 
168 
creates a deep profile at $x_{0}=\frac{width}{2}$; 
169 
\begin{python} 
170 
cut=int(len(xiq)/2) 
171 
pl.plot(ziq[:,cut]*1000.,yiq) 
172 
pl.title("Vertical Heat Flow Depth Profile") 
173 
pl.xlabel("Heat Flow (mW/m^2)") 
174 
pl.ylabel("Depth (m)") 
175 
pl.savefig(os.path.join(save_path,"hf.png")) 
176 
\end{python} 
177 
This process can be repeated for various variations of the solution. 
178 
\reffig{figs:dps} shows temperature, temperature gradient, thermal conductivity 
179 
and heat flow. 
180 

181 
\begin{figure}[htp] 
182 
\centering 
183 
\subfigure[Temperature Depth 
184 
Profile]{\includegraphics[width=3in]{figures/heatrefractiontdp.png}\label{ 
185 
fig:tdp}} 
186 
\subfigure[Temperature Gradient Depth 
187 
Profile]{\includegraphics[width=3in]{figures/heatrefractiontgdp.png}\label{ 
188 
fig:tgdp}} 
189 
\subfigure[Thermal Conductivity 
190 
Profile]{\includegraphics[width=3in]{figures/heatrefractiontcdp.png}\label{ 
191 
fig:tcdp}} 
192 
\subfigure[Heat Flow Depth 
193 
Profile]{\includegraphics[width=3in]{figures/heatrefractionhf.png}\label{fig:hf} 
194 
} 
195 
\caption{Example 5b: Depth profiles down centre of model} 
196 
\label{figs:dps} 
197 
\end{figure} 
198 
\clearpage 
199 

200 
\section{Arrow Plots in \mpl} 
201 
\sslist{example05c.py and cblib.py} 
202 
The distribution of the flux $\kappa \nabla T$ is now visualised over the 
203 
domain and the results are plotted in \reffig{fig:hr001qumodel}. 
204 
The plot puts together three components. A contour plot of the temperature, 
205 
a coloured representation of the two subdomains where colour represents the 
206 
thermal conductivity in the particular region, and finally the arrows 
207 
representing the local direction of the steepest gradient of the flux. 
208 

209 
Contours have already been discussed in \refSec{Sec:2DHD plot}. To show 
210 
subdomains, we need to go back to \pycad data to get the points used to 
211 
describe the boundary of the subdomains. We have created the \verbCurveLoop 
212 
class object \verbtblockloop to define the boundary of the upper subdomain. 
213 
We use the \verbgetPolygon() method of \verbCurveLoop to get 
214 
access to the \verbPoints used to define the boundary. The statement 
215 
\begin{python} 
216 
[ p.getCoordinates() for p in tblockloop.getPolygon() ] 
217 
\end{python} 
218 
creates a list of the node coordinates of all the points in question. In order 
219 
to simplify the selection of the $x$ and $y$ coordinates the list is converted 
220 
into \modnumpy array. To add the area coloured in brown to the plot we use; 
221 
\begin{python} 
222 
import pylab as pl 
223 
import numarray as np 
224 
tpg=np.array([p.getCoordinates() for p in tblockloop.getPolygon() ]) 
225 
pl.fill(tpg[:,0],tpg[:,1],'brown',label='2 W/m/k',zorder=1000) 
226 
\end{python} 
227 
The same code is applied to \verbbblockloop to create the red area for this 
228 
subdomain. 
229 

230 
To plot vectors representing the flux orientation we use the 
231 
\verbquiver function in \pylab. The function places vectors at locations in 
232 
the domain. 
233 
For instance one can plot vectors at the locations of the sample points used by 
234 
\esc to represent the flux \verbkappa*grad(T). As a vector is plotted at 
235 
each sample point one typically ends up with too many vectors. So one needs to 
236 
select a subset of points as follows: 
237 

238 
First we create a coarse grid of points on a rectangular mesh, e.g. $20 \times 
239 
20$ points. Here we choose a grid of points which are located at the centre of 
240 
a \verbnx $\times$ \verbny grid; 
241 
\begin{python} 
242 
dx = width/nx # x spacing 
243 
dy = depth/ny # y spacing 
244 
grid = [ ] # the grid points 
245 
for j in xrange(0,ny1): 
246 
for i in xrange(0,nx1): 
247 
grid.append([dx/2+dx*i,dy/2+dy*j]) 
248 
\end{python} 
249 
With the \verbLocator function \esc provides a mechanism to identify sample 
250 
points that are closest to the grid points we have selected and to retrieve the 
251 
data at these points; 
252 
\begin{python} 
253 
from esys.escript.pdetools import Locator 
254 
flux=kappa*grad(T) 
255 
fluxLoc = Locator(flux.getFunctionSpace(),grid) 
256 
subflux= fluxLoc(flux) 
257 
\end{python} 
258 
\verbsubflux now contains a list of flux components at certain sample points. 
259 
To get the list of the sample point coordinates one can use the \verbgetX() 
260 
method of the \verbLocator; 
261 
\begin{python} 
262 
subfluxloc = fluxLoc.getX() 
263 
\end{python} 
264 
To simplify the selection of $x$ and $y$ components it is convenient to 
265 
transform \verbsubflux and \verbsubfluxloc to \numpy arrays 
266 
\verbxflux, \verbflux. 
267 
This function is implemented in the \verbsubsample function within the 
268 
\file{clib.py} file so we can use it in other examples. One can easily use this 
269 
function to create a vector plot of the flux; 
270 
\begin{python} 
271 
from cblib import subsample 
272 
xflux, flux=subsample(kappa*grad(T), nx=20, ny=20) 
273 
pl.quiver(xflux[:,0],xflux[:,1],flux[:,0],flux[:,1], angles='xy',color="white") 
274 
\end{python} 
275 
Finally, we add a title and labels; 
276 
\begin{python} 
277 
pl.title("Heat Refraction across a clinal structure.") 
278 
pl.xlabel("Horizontal Displacement (m)") 
279 
pl.ylabel("Depth (m)") 
280 
pl.title("Heat Refraction across a clinal structure \n with gradient quivers.") 
281 
pl.savefig(os.path.join(saved_path,"flux.png")) 
282 
\end{python} 
283 
to get the desired result, see \reffig{fig:hr001qumodel}. 
284 

285 
\begin{figure}[ht] 
286 
\centerline{\includegraphics[width=5.in]{figures/heatrefractionflux}} 
287 
\caption{Example 5c: Heat refraction model with gradient indicated by vectors} 
288 
\label{fig:hr001qumodel} 
289 
\end{figure} 
290 
\clearpage 
291 

292 
\section{Example 6:Fault and Overburden Model} 
293 
\sslist{example06.py and cblib.py} 
294 
A slightly more complicated model can be found in the example file 
295 
\textit{heatrefraction2_solver.py} where three blocks are used within the 
296 
model, see~\reffig{fig:hr002qumodel}. It is left to the reader to work through 
297 
this example. 
298 

299 
\begin{figure}[ht] 
300 
\centerline{\includegraphics[width=4.in]{figures/heatrefraction2flux}} 
301 
\caption{Example 6: Heat refraction model with three blocks and heat flux} 
302 
\label{fig:hr002qumodel} 
303 
\end{figure} 
304 
\clearpage 