1 

2 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
3 
% Copyright (c) 20032018 by The University of Queensland 
4 
% http://www.uq.edu.au 
5 
% 
6 
% Primary Business: Queensland, Australia 
7 
% Licensed under the Apache License, version 2.0 
8 
% http://www.apache.org/licenses/LICENSE2.0 
9 
% 
10 
% Development until 2012 by Earth Systems Science Computational Center (ESSCC) 
11 
% Development 20122013 by School of Earth Sciences 
12 
% Development from 2014 by Centre for Geoscience Computing (GeoComp) 
13 
% 
14 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
15 

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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