1 
ahallam 
2597 

2 


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
3 


% 
4 
jfenwick 
2881 
% Copyright (c) 20032010 by University of Queensland 
5 
ahallam 
2597 
% 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 
gross 
2870 
\section{Steadystate Heat Refraction} 
15 


\label{STEADYSTATE HEAT REFRACTION} 
16 
gross 
2931 

17 
ahallam 
2975 
In this chapter we demonstrate how to handle more complex geometries. 
18 
gross 
2931 

19 
artak 
2972 
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 generate finite element meshes 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 \esc to define properties like material constants and source locations. 
20 
gross 
2931 

21 
ahallam 
2977 
We proceed in this chapter by first looking at a very simple geometry. Whilst a simple rectangular domain is not very interesting the example is elaborated upon later by introducing an internal curved interface. 
22 
ahallam 
2597 

23 
gross 
2948 
\begin{figure}[ht] 
24 


\centerline{\includegraphics[width=4.in]{figures/pycadrec}} 
25 
gross 
2952 
\caption{Example 4: Rectangular Domain for \pycad.} 
26 
gross 
2948 
\label{fig:pycad rec} 
27 


\end{figure} 
28 
ahallam 
2634 

29 
artak 
2972 
\section{Example 4: Creating the domain with \pycad} 
30 
ahallam 
2977 
\label{example4} 
31 
gross 
2951 
\sslist{example04a.py} 
32 
ahallam 
2597 

33 
artak 
2978 
We modify the example in Chapter~\ref{CHAP HEAT 2a} in two ways: we look the steady state 
34 


case with slightly modified boundary conditions and use a more flexible tool 
35 


to generate to generate the geometry. Lets look at the geometry first. 
36 
ahallam 
2597 

37 
ahallam 
2975 
We want to define a rectangular domain of width $5 km$ and depth $6 km$ below the surface of the Earth. The domain is subject to a few conditions. The temperature is known at the surface and the basement has a known heat flux. Each side of the domain is insulated and the aim is to calculate the final temperature distribution. 
38 
ahallam 
2597 

39 
gross 
2948 
In \pycad there are a few primary constructors that build upon each other to define domains and boundaries; 
40 
artak 
2972 
the ones we use are: 
41 
ahallam 
2801 
\begin{python} 
42 
gross 
2948 
from esys.pycad import * 
43 
ahallam 
2597 
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 
gross 
2948 
PlaneSurface() #Creates a surface based on a CurveLoop 
47 
ahallam 
2801 
\end{python} 
48 
ahallam 
2975 
So to construct our domain as shown in \reffig{fig:pycad rec}, we first need to create 
49 
gross 
2948 
the corner points. From the corner points we build the four edges of the rectangle. The four edges 
50 
ahallam 
2975 
then form a closed loop which defines our domain as a surface. 
51 
ahallam 
2634 
We start by inputting the variables we need to construct the model. 
52 
ahallam 
2801 
\begin{python} 
53 
gross 
2948 
width=5000.0*m #width of model 
54 


depth=6000.0*m #depth of model 
55 
ahallam 
2801 
\end{python} 
56 
gross 
2948 
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. 
57 
ahallam 
2801 
\begin{python} 
58 
ahallam 
2597 
# Overall Domain 
59 
ahallam 
2634 
p0=Point(0.0, 0.0, 0.0) 
60 


p1=Point(0.0, depth, 0.0) 
61 


p2=Point(width, depth, 0.0) 
62 


p3=Point(width, 0.0, 0.0) 
63 
ahallam 
2801 
\end{python} 
64 
gross 
2948 
Now lines are defined using our points. This forms a rectangle around our domain; 
65 
ahallam 
2801 
\begin{python} 
66 
ahallam 
2597 
l01=Line(p0, p1) 
67 


l12=Line(p1, p2) 
68 


l23=Line(p2, p3) 
69 


l30=Line(p3, p0) 
70 
ahallam 
2801 
\end{python} 
71 
gross 
2948 
Notice that lines have a direction. These lines form the basis for our domain boundary, which is a closed loop. 
72 
ahallam 
2801 
\begin{python} 
73 
ahallam 
2597 
c=CurveLoop(l01, l12, l23, l30) 
74 
ahallam 
2801 
\end{python} 
75 
gross 
2948 
Be careful to define the curved loop in an \textbf{anticlockwise} manner otherwise the meshing algorithm may fail. 
76 


Finally we can define the domain as 
77 
ahallam 
2801 
\begin{python} 
78 
gross 
2948 
rec = PlaneSurface(c) 
79 


\end{python} 
80 


At this point the introduction of the curved loop seems to be unnecessary but this concept plays an important role 
81 


if holes are introduced. 
82 



83 


Now we are ready to handover the domain \verbrec to a mesher which subdivides the domain into triangles (or tetrahedron in 3D). In our case we use \gmsh. We create 
84 


an instance of the \verbDesign class which will handle the interface to \gmsh: 
85 


\begin{python} 
86 


from esys.pycad.gmsh import Design 
87 


d=Design(dim=2, element_size=200*m) 
88 


\end{python} 
89 
ahallam 
2977 
The argument \verbdim defines the spatial dimension of the domain\footnote{If \texttt{dim}=3 the rectangle would be interpreted as a surface in the three dimensional space.}. The second argument \verbelement_size defines element size which is the maximum length of a triangle edge in the mesh. The element size needs to be chosen with care in order to avoid very dense meshes. If the mesh is too dense computational time will be long but if the mesh is too sparse the model result will be poor. In our case with an element size of $200$m 
90 
artak 
2972 
and a domain length of $6000$m we will end up with about $\frac{6000m}{200m}=30$ triangles in each spatial direction. So we end up with about $30 \times 30 = 900$ triangles which is a size that can be handled easily. 
91 
gross 
2948 
We can easily add our domain \verbrec to the \verbDesign; 
92 


\begin{python} 
93 


d.addItem(rec) 
94 


\end{python} 
95 


We have the plan to set a heat flux on the bottom of the domain. One can use the masking technique to do this 
96 
artak 
2972 
but \pycad offers a more convenient technique called tagging. With this technique items in the domain are 
97 
gross 
2948 
named using the \verbPropertySet class. We can then later use this name to set values secifically for 
98 


those sample points located on the named items. Here we name the bottom face of the 
99 


domain where we will set the heat influx: 
100 


\begin{python} 
101 


ps=PropertySet("linebottom",l12)) 
102 


d.addItem(ps) 
103 


\end{python} 
104 


Now we are ready to hand over the \verbDesign to \FINLEY: 
105 


\begin{python} 
106 


from esys.finley import MakeDomain 
107 


domain=MakeDomain(d) 
108 


\end{python} 
109 


The \verbdomain object can now be used in the same way like the return object of the \verbRectangle 
110 


object we have used previously to generate a mesh. It is common practice to separate the 
111 


mesh generation from the PDE solution. The main reason for this is that mesh generation can be computationally very expensive in particular in 3D. So it is more efficient to generate the mesh once and write it to a file. The mesh 
112 


can then be read in every time a new simulation is run. \FINLEY supports this in the following 
113 


way~\footnote{An alternative are the \texttt{dump} and \texttt{load} functions. They using a binary format and tend to be much smaller.} 
114 


\begin{python} 
115 


# write domain to an text file 
116 
gross 
2951 
domain.write("example04.fly") 
117 
gross 
2948 
\end{python} 
118 
artak 
2972 
and then for reading in another script; 
119 
gross 
2948 
\begin{python} 
120 


# read domain from text file 
121 


from esys.finley import ReadMesh 
122 
gross 
2951 
domain =ReadMesh("example04.fly") 
123 
gross 
2948 
\end{python} 
124 



125 


\begin{figure}[ht] 
126 


\centerline{\includegraphics[width=4.in]{figures/simplemesh}} 
127 
gross 
2952 
\caption{Example 4a: Mesh over rectangular domain, see \reffig{fig:pycad rec}.} 
128 
gross 
2948 
\label{fig:pycad rec mesh} 
129 


\end{figure} 
130 



131 
ahallam 
2977 
Before we discuss how to solve the PDE for this 
132 


problem, it is useful to present two additional options of the \verbDesign class. 
133 
artak 
2972 
They allow the user accessing the script which is used by \gmsh to generate the mesh as well as 
134 
gross 
2948 
the mesh as it has been generated by \gmsh. This is done by setting specific names for these files: 
135 


\begin{python} 
136 
gross 
2951 
d.setScriptFileName("example04.geo") 
137 


d.setMeshFileName("example04.msh") 
138 
gross 
2948 
\end{python} 
139 


Usually the extension \texttt{geo} is used for the script file of the \gmsh geometry and 
140 


the extension \texttt{msh} for the mesh file. Normally these files are deleted after usage. 
141 
ahallam 
2977 
Accessing these files can be helpful to debug the generation of more complex geometries. The geometry and the mesh can be visualised from the command line using 
142 
gross 
2948 
\begin{verbatim} 
143 
gross 
2951 
gmsh example04.geo # show geometry 
144 


gmsh example04.msh # show mesh 
145 
gross 
2948 
\end{verbatim} 
146 


The mesh is shown in \reffig{fig:pycad rec mesh}. 
147 


\begin{figure}[ht] 
148 


\centerline{\includegraphics[width=4.in]{figures/simpleheat}} 
149 
gross 
2952 
\caption{Example 4b: Result of simple steady state heat problem.} 
150 
gross 
2948 
\label{fig:steady state heat} 
151 


\end{figure} 
152 



153 



154 


\section{The Steadystate Heat Equation} 
155 
gross 
2951 
\sslist{example04b.py, cblib} 
156 
ahallam 
2977 
A temperature equilibrium or steady state is reached when the temperature distribution in the model does not change with time. To calculate the 
157 


the steady state solution the time derivative term in \refEq{eqn:Tform nabla} needs to be set to zero; 
158 
gross 
2948 
\begin{equation}\label{eqn:Tform nabla steady} 
159 


\nabla \cdot \kappa \nabla T = q\hackscore H 
160 


\end{equation} 
161 
ahallam 
2977 
This PDE is easier to solve than the PDE in \refEq{eqn:hdgenf2}, as no time steps (iterations) are required. The \verbD term can be dropped from the equation; 
162 
gross 
2948 
\begin{python} 
163 


mypde=LinearPDE(domain) 
164 


mypde.setValue(A=kappa*kronecker(model), Y=qH) 
165 


\end{python} 
166 


The temperature at the top face of the domain is known as \verbTtop (=$20 C$). In \refSec{Sec:1DHDv0} we have 
167 


already discussed how this constraint is added to the PDE: 
168 


\begin{python} 
169 


x=Solution(domain).getX() 
170 


mypde.setValue(q=whereZero(x[1]sup(x[1])),r=Ttop) 
171 


\end{python} 
172 
artak 
2972 
Notice that we use the \verbsup function to calculate the maximum of $y$ coordinates of the relevant sample points. 
173 
gross 
2948 

174 


In all cases so far we have assumed that the domain is insulated which translates 
175 


into a zero normal flux $n \cdot \kappa \nabla T$, see \refEq{eq:hom flux}. In the modeling 
176 


setup of this chapter we want to set the normal heat flux at the bottom to \verbqin while we still 
177 


maintain insulation at the left and right face. Mathematically we can express this as the format 
178 


\begin{equation} 
179 


n \cdot \kappa \nabla T = q\hackscore{S} 
180 


\label{eq:inhom flux} 
181 


\end{equation} 
182 


where $q\hackscore{S}$ is a function of its location on the boundary. Its value becomes zero 
183 


for locations on the left or right face of the domain while it has the value \verbqin at the bottom face. 
184 


Notice that the value of $q\hackscore{S}$ at the top face is not relevant as we prescribe the temperature here. 
185 
artak 
2972 
We could define $q\hackscore{S}$ by using the masking techniques demonstrated earlier. The tagging mechanism provides an alternative and in many cases more convenient way of defining piecewise 
186 
ahallam 
2977 
constant functions such as $q\hackscore{S}$. Recall now that the bottom face was 
187 


denoted with the name \verblinebottom when we defined the domain. 
188 
gross 
2948 
We can use this now to create $q\hackscore{S}$; 
189 


\begin{python} 
190 


qS=Scalar(0,FunctionOnBoundary(domain)) 
191 


qS.setTaggedValue("linebottom",qin) 
192 


\end{python} 
193 


In the first line \verbqS is defined as a scalar value over the sample points on the boundary of the domain. It is 
194 


initialized to zero for all sample points. In the second statement the values for those sample points which 
195 


on the line marked by \verblinebottom are set to \verbqin. 
196 



197 
ahallam 
2977 
The Neuman boundary condition assumed by \esc has the form 
198 
gross 
2948 
\begin{equation}\label{NEUMAN 2b} 
199 


n\cdot A \cdot\nabla u = y 
200 


\end{equation} 
201 


In comparison to the version in \refEq{NEUMAN 2} we have used so far the right hand side is now 
202 
ahallam 
2977 
the new PDE coefficient $y$. As we have not specified $y$ in our previous examples, \esc has assumed 
203 
gross 
2948 
the value zero for $y$. A comparison of \refEq{NEUMAN 2b} and \refEq{eq:inhom flux} reveals that one need to 
204 


choose $y=q\hackscore{S}$; 
205 


\begin{python} 
206 


qS=Scalar(0,FunctionOnBoundary(domain)) 
207 


qS.setTaggedValue("linebottom",qin) 
208 


mypde.setValue(y=qS) 
209 


\end{python} 
210 


To plot the results we are using the \modmpl library as shown \refSec{Sec:2DHD plot}. For convenience 
211 


the interpolation of the temperature to a rectangular grid for contour plotting is made available 
212 


via the \verbtoRegGrid function in the \verbcblib module. Your result should look similar to 
213 


\reffig{fig:steady state heat}. 
214 



215 


\begin{figure}[ht] 
216 


\centerline{\includegraphics[width=4.in]{figures/anticlineheatrefraction}} 
217 
gross 
2952 
\caption{Example 5a: Heat refraction model with point and line labels.} 
218 
gross 
2948 
\label{fig:anticlinehrmodel} 
219 


\end{figure} 
220 



221 


\begin{figure}[ht] 
222 


\centerline{\includegraphics[width=4.in]{figures/heatrefraction}} 
223 
gross 
2952 
\caption{Example 5a: Temperature Distribution in the Heat Refraction Model.} 
224 
gross 
2948 
\label{fig:anticlinetemp} 
225 


\end{figure} 
226 



227 
gross 
2952 
\section{Example 5: A Heat Refraction model} 
228 
ahallam 
2977 
\label{example5} 
229 
gross 
2952 
\sslist{example05a.py and cblib.py} 
230 
ahallam 
2977 
Our heat refraction model will be a large anticlinal structure that is subject to a constant temperature at the surface and experiencing 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. The heat flow pattern demonstrates the dependence upon the material properties and the shape of the interface. 
231 
gross 
2948 

232 
ahallam 
2977 
The script of Example \ref{example4m} is modifed by subdividing the block into two parts. The curve 
233 
gross 
2948 
separating the two blocks is given as a spline, see \reffig{fig:anticlinehrmodel}. The data points 
234 
ahallam 
2977 
used to define the curve may be used from a database of measurements, but for simplicity it is assumed here that there coordinates are 
235 


known in an analytic form. 
236 
gross 
2948 

237 
ahallam 
2977 
There are two modes available in this example. When \verb modal=1 , this indicates to the script that the model should be an anticline. Otherwise, when \verb modal=1 , thel model is 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. 
238 
gross 
2948 

239 


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



241 
artak 
2978 
The curve defining our clinal structure is located approximately in 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. 
242 
gross 
2948 
\begin{python} 
243 
ahallam 
2634 
# Material Boundary 
244 


x=[ Point(i*dsp\ 
245 


,dep_sp+modal*orit*h_sp*cos(pi*i*dsp/dep_sp+pi))\ 
246 


for i in range(0,sspl)\ 
247 


] 
248 
ahallam 
2597 
mysp = Spline(*tuple(x)) #*tuple() forces x to become a tuple 
249 
ahallam 
2801 
\end{python} 
250 
ahallam 
2597 
The start and end points of the spline can be returned to help define the material boundaries. 
251 
ahallam 
2801 
\begin{python} 
252 
gross 
2948 
x1=mysp.getStartPoint() 
253 


x2=mysp.getEndPoint() 
254 
ahallam 
2801 
\end{python} 
255 
ahallam 
2977 
The top block or material above the clinal/spline boundary is defined in an \textbf{anticlockwise} manner by creating lines and then a closed loop. By meshing the subdomain we also need to assign it a planar surface. 
256 
ahallam 
2801 
\begin{python} 
257 
ahallam 
2597 
# TOP BLOCK 
258 


tbl1=Line(p0,x1) 
259 


tbl2=mysp 
260 


tbl3=Line(x2,p3) 
261 
gross 
2948 
l30=Line(p3, p0) 
262 
ahallam 
2597 
tblockloop = CurveLoop(tbl1,tbl2,tbl3,l30) 
263 


tblock = PlaneSurface(tblockloop) 
264 
ahallam 
2801 
\end{python} 
265 
ahallam 
2977 
This process is repeated for every other subdomain. In this example there is only one other, the bottom block. The process is similar to the top block but with a few differences. The spline points must be reversed by setting the spline as negative. 
266 
ahallam 
2801 
\begin{python} 
267 
ahallam 
2597 
bbl4=mysp 
268 
ahallam 
2801 
\end{python} 
269 


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. 
270 


\begin{python} 
271 
ahallam 
2597 
#clockwise check 
272 
gross 
2948 
bblockloop=CurveLoop(mysp,Line(x2,p2),Line(p2,p1),Line(p1,x1)) 
273 
ahallam 
2801 
\end{python} 
274 
ahallam 
2977 
The last few steps in creating the domain require that the previously defined domain and subdomain points are submitted to generate a mesh that can be imported into \esc. 
275 
ahallam 
2801 
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. 
276 


\begin{python} 
277 
ahallam 
2597 
# Create a Design which can make the mesh 
278 


d=Design(dim=2, element_size=200) 
279 
ahallam 
2801 
# Add the subdomains and flux boundaries. 
280 
ahallam 
2658 
d.addItems(PropertySet("top",tblock),PropertySet("bottom",bblock),\ 
281 


PropertySet("linebottom",l12)) 
282 
ahallam 
2801 
# Create the geometry, mesh and \esc domain 
283 
gross 
2952 
d.setScriptFileName(os.path.join(save_path,"example05.geo")) 
284 


d.setMeshFileName(os.path.join(save_path,"example05.msh")) 
285 
gross 
2948 
domain=MakeDomain(d,optimizeLabeling=True) 
286 
ahallam 
2801 
\end{python} 
287 
gross 
2948 
The creation of our domain and its mesh is complete. 
288 
ahallam 
2597 

289 
gross 
2948 
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. 
290 
ahallam 
2801 
\begin{python} 
291 
ahallam 
2634 
# set up kappa (thermal conductivity across domain) using tags 
292 
gross 
2948 
kappa=Scalar(0,Function(domain)) 
293 


kappa.setTaggedValue("top",2.0*W/m/K) 
294 


kappa.setTaggedValue("bottom",4.0*W/m/K) 
295 
ahallam 
2801 
\end{python} 
296 
gross 
2948 
No further changes are required to the PDE solution step, see \reffig{fig:anticlinetemp} for the result. 
297 
ahallam 
2634 

298 
gross 
2948 
\begin{figure} 
299 


\centering 
300 


\subfigure[Temperature Depth Profile]{\includegraphics[width=3in]{figures/heatrefractiontdp.png}\label{fig:tdp}} 
301 


\subfigure[Temperature Gradient Depth Profile]{\includegraphics[width=3in]{figures/heatrefractiontgdp.png}\label{fig:tgdp}} 
302 


\subfigure[Thermal Conductivity Profile]{\includegraphics[width=3in]{figures/heatrefractiontcdp.png}\label{fig:tcdp}} 
303 


\subfigure[Heat Flow Depth Profile]{\includegraphics[width=3in]{figures/heatrefractionhf.png}\label{fig:hf}} 
304 
gross 
2952 
\caption{Example 5b: Depth profiles down centre of model.} 
305 
gross 
2948 
\label{figs:dps} 
306 


\end{figure} 
307 
ahallam 
2634 

308 
gross 
2948 
\section{Line profiles of 2D data} 
309 
gross 
2952 
\sslist{example05b.py and cblib.py} 
310 
gross 
2948 
We want to investigate the profile of the data of the last example. 
311 
artak 
2972 
We are particularly interested in the depth profile of the heat flux which is 
312 
ahallam 
2977 
the second component of $\kappa \nabla T$. The script from the previous section is extended 
313 


to show how a vertical profile can be plotted. 
314 
ahallam 
2634 

315 
ahallam 
2977 
The first important piece of information, is that \esc assumes that $\kappa \nabla T$ is not smooth and 
316 


that these values are defined at numerical interpolation points. This assumption is reasonable as 
317 
gross 
2948 
the flux is the product of the piecewise constant function $\kappa$ and 
318 
ahallam 
2977 
the gradient of the temperature $T$ which has a discontinuity at the rock interface. 
319 


Before plotting this function we need to smooth the solution using the 
320 
artak 
2972 
\verbProjector() class; 
321 
ahallam 
2801 
\begin{python} 
322 
gross 
2948 
from esys.escript.pdetools import Projector 
323 


proj=Projector(domain) 
324 


qu=proj(kappa*grad(T)) 
325 
ahallam 
2801 
\end{python} 
326 
ahallam 
2977 
The \verbproj object provides a mechanism to distribute values given at the numerical interpolation points to the nodes 
327 


of the FEM mesh  the heat flux in this example. \verbqu has the same function space 
328 


as the temperature \verbT. The smoothed flux is interpolated 
329 


to a regular $200\times 200$ grid via; 
330 
ahallam 
2801 
\begin{python} 
331 
gross 
2948 
xiq,yiq,ziq = toRegGrid(qu[1],200,200) 
332 
ahallam 
2801 
\end{python} 
333 
artak 
2972 
using the \verbtoRegGrid function from the cookbook library which we are using for the contour plot. 
334 
gross 
2948 
At return \verbziq[j,i] is the value of vertical heat flux at point 
335 


(\verbxiq[i],\verbyiq[j]). We can easily create deep profiles now by 
336 


plotting slices \verbziq[:,i] over \verbyiq. The following script 
337 


creates a deep profile at $x_{0}=\frac{width}{2}$; 
338 
ahallam 
2801 
\begin{python} 
339 
gross 
2948 
cut=int(len(xiq)/2) 
340 


pl.plot(ziq[:,cut]*1000.,yiq) 
341 


pl.title("Vertical Heat Flow Depth Profile") 
342 


pl.xlabel("Heat Flow (mW/m^2)") 
343 
ahallam 
2681 
pl.ylabel("Depth (m)") 
344 
gross 
2952 
pl.savefig(os.path.join(save_path,"hf.png")) 
345 
ahallam 
2801 
\end{python} 
346 


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}. 
347 
gross 
2948 

348 


\begin{figure}[ht] 
349 


\centerline{\includegraphics[width=5.in]{figures/heatrefractionflux}} 
350 
gross 
2952 
\caption{Example 5c: Heat refraction model with gradient indicated by vectors.} 
351 
gross 
2948 
\label{fig:hr001qumodel} 
352 
ahallam 
2681 
\end{figure} 
353 



354 
gross 
2948 
\section{Arrow plots in \mpl} 
355 
gross 
2952 
\sslist{example05c.py and cblib.py} 
356 
ahallam 
2977 
The distribution of the flux $\kappa \nabla T$ is now visualised over the domain 
357 


and the results are plotted in \reffig{fig:hr001qumodel}. 
358 
gross 
2948 
The plot puts together three components. A contour plot of the temperature, 
359 


a colored representation of the two subdomains where colour represents the thermal conductivity 
360 
ahallam 
2977 
in the particular region and finally the arrows representing the local direction of the steepest gradient of the flux. 
361 
ahallam 
2681 

362 
gross 
2948 
Contours have already been discussed in \refSec{Sec:2DHD plot}. To show subdomains, 
363 


we need to go back to \pycad data to get the points used to describe the boundary of the 
364 


subdomains. We have created the \verbCurveLoop class object 
365 


\verbtblockloop to define the boundary of the upper subdomain. 
366 


We use the \verbgetPolygon() method of \verbCurveLoop to get 
367 


access to the \verbPoints used top define the boundary. The statement 
368 
ahallam 
2801 
\begin{python} 
369 
gross 
2948 
[ p.getCoordinates() for p in tblockloop.getPolygon() ]) 
370 
ahallam 
2801 
\end{python} 
371 
gross 
2948 
creates a list of the node coordinates of all the points in question. In order 
372 


to simplify the selection of the $x$ and $y$ coordinates the list is converted 
373 


into \modnumpy array. To add the area colored in brown to the plot we use; 
374 
ahallam 
2801 
\begin{python} 
375 
gross 
2948 
import pylab as pl 
376 


import numarray as np 
377 


tpg=np.array([p.getCoordinates() for p in tblockloop.getPolygon() ]) 
378 


pl.fill(tpg[:,0],tpg[:,1],'brown',label='2 W/m/k',zorder=1000) 
379 


\end{python} 
380 
ahallam 
2977 
The same code is applied to \verbbblockloop to create the red area for this subdomain to the block. 
381 
ahallam 
2634 

382 
gross 
2948 
To plot vectors representing the flux orientation we use the 
383 


\verbquiver function in \pylab. The function places vectors at locations in the domain. 
384 


For instance one can plot vectors at the locations of the sample points used by \esc 
385 


to represent the flux \verbkappa*grad(T). As a vector is plotted at each sample point one typically ends 
386 


up with two many vectors. So one needs to select a subset of points: 
387 


First we create a coarse grid of points on a rectangular mesh, e.g. $20 \times 20$ points. Here we choose a grid of points which are located at the center of a \verbnx $\times$ \verbny grid; 
388 


\begin{python} 
389 


dx = width/nx # x spacing 
390 


dy = depth/ny # y spacing 
391 


grid = [ ] # the grid points 
392 


for j in xrange(0,ny1): 
393 


for i in xrange(0,nx1): 
394 


grid.append([dx/2+dx*i,dy/2+dy*j]) 
395 
ahallam 
2801 
\end{python} 
396 
ahallam 
2977 
With the \verbLocator function \esc provides a mechanism to identify sample points that are closest 
397 


to the the grid points we have selected and to rettieve the data at these data points; 
398 
ahallam 
2801 
\begin{python} 
399 
gross 
2948 
from esys.escript.pdetools import Locator 
400 


flux=kappa*grad(T) 
401 


fluxLoc = Locator(flux.getFunctionSpace(),grid) 
402 


subflux= fluxLoc(flux) 
403 


\end{python} 
404 


\verbsubflux gives now a list of flux component at certain sample points. To get the 
405 


list of the sample point coordinates one can use the \verbgetX() method of the 
406 


\verbLocator; 
407 


\begin{python} 
408 


subfluxloc = fluxLoc.getX() 
409 


\end{python} 
410 


To simplify the selection of $x$ and $y$ components it is convenient 
411 


to transform \verbsubflux and \verbsubfluxloc to \numpy arrays 
412 


\verbxflux, \verbflux. 
413 


This function is implemented in the \verbsubsample 
414 


in the \file{clib.py} file so we can use it in other examples. One can easily use this function 
415 


to create a vector plot of the flux; 
416 


\begin{python} 
417 


from cblib import subsample 
418 


xflux, flux=subsample(kappa*grad(T), nx=20, ny=20) 
419 


pl.quiver(xflux[:,0],xflux[:,1],flux[:,0],flux[:,1], angles='xy',color="white") 
420 


\end{python} 
421 


We add title and labels; 
422 


\begin{python} 
423 
ahallam 
2634 
pl.title("Heat Refraction across a clinal structure.") 
424 


pl.xlabel("Horizontal Displacement (m)") 
425 


pl.ylabel("Depth (m)") 
426 


pl.title("Heat Refraction across a clinal structure \n with gradient quivers.") 
427 
gross 
2952 
pl.savefig(os.path.join(saved_path,"flux.png")) 
428 
ahallam 
2801 
\end{python} 
429 
gross 
2948 
to get the desired result, see \reffig{fig:hr001qumodel}. 
430 
ahallam 
2634 

431 
ahallam 
2681 
\begin{figure}[ht] 
432 
gross 
2948 
\centerline{\includegraphics[width=4.in]{figures/heatrefraction2flux}} 
433 
gross 
2952 
\caption{Example 6: Heat refraction model with three blocks and heat flux.} 
434 
gross 
2948 
\label{fig:hr002qumodel} 
435 
ahallam 
2658 
\end{figure} 
436 



437 
gross 
2952 
\section{Example 6:Fault and Overburden Model} 
438 


\sslist{example06.py and cblib.py} 
439 
gross 
2948 
A slightly more complicated model can be found in the examples \textit{heatrefraction2_solver.py} where three blocks are used within the model, see~\reffig{fig:hr002qumodel}. It is left to the reader to work through this example. 
440 



441 


