1 
ahallam 
2597 

2 
jfenwick 
3989 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
3 
jfenwick 
3911 
% Copyright (c) 20032012 by University of Queensland 
4 
jfenwick 
3989 
% http://www.uq.edu.au 
5 
ahallam 
2597 
% 
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 
jfenwick 
3989 
% Development until 2012 by Earth Systems Science Computational Center (ESSCC) 
11 


% Development since 2012 by School of Earth Sciences 
12 


% 
13 


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
14 
ahallam 
2597 

15 
gross 
2870 
\section{Steadystate Heat Refraction} 
16 


\label{STEADYSTATE HEAT REFRACTION} 
17 
gross 
2931 

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

20 
ahallam 
2979 
Steadystate heat refraction will give us an opportunity to investigate some of 
21 


the richer features that the \esc package has to offer. One of these is \pycad . 
22 


The advantage of using \pycad is that it offers an easy method for developing 
23 


and manipulating complex domains. In conjunction with \gmsh we can generate 
24 


finite element meshes that conform to our domain's shape providing accurate 
25 


modelling of interfaces and boundaries. Another useful function of \pycad is 
26 


that we can tag specific areas of our domain with labels as we construct them. 
27 


These labels can then be used in \esc to define properties like material 
28 


constants and source locations. 
29 
gross 
2931 

30 
ahallam 
2979 
We proceed in this chapter by first looking at a very simple geometry. Whilst a 
31 


simple rectangular domain is not very interesting the example is elaborated upon 
32 


later by introducing an internal curved interface. 
33 
ahallam 
2597 

34 
caltinay 
2982 
\section{Example 4: Creating the Domain with \pycad} 
35 
ahallam 
2977 
\label{example4} 
36 
gross 
2951 
\sslist{example04a.py} 
37 
caltinay 
2982 
We modify the example in Chapter~\ref{CHAP HEAT 2a} in two ways: we look at the 
38 


steady state case with slightly modified boundary conditions and use a more 
39 


flexible tool to generate the geometry. Let us look at the geometry first. 
40 
ahallam 
2597 

41 
ahallam 
2979 
We want to define a rectangular domain of width $5 km$ and depth $6 km$ below 
42 


the surface of the Earth. The domain is subject to a few conditions. The 
43 


temperature is known at the surface and the basement has a known heat flux. Each 
44 


side of the domain is insulated and the aim is to calculate the final 
45 


temperature distribution. 
46 
ahallam 
2597 

47 
ahallam 
2979 
In \pycad there are a few primary constructors that build upon each other to 
48 
caltinay 
2982 
define domains and boundaries. The ones we use are: 
49 
ahallam 
2801 
\begin{python} 
50 
gross 
2948 
from esys.pycad import * 
51 
ahallam 
2597 
Point() #Create a point in space. 
52 


Line() #Creates a line from a number of points. 
53 


CurveLoop() #Creates a closed loop from a number of lines. 
54 
gross 
2948 
PlaneSurface() #Creates a surface based on a CurveLoop 
55 
ahallam 
2801 
\end{python} 
56 
ahallam 
2979 
So to construct our domain as shown in \reffig{fig:pycad rec}, we first need to 
57 
caltinay 
2982 
create the corner points. From the corner points we build the four edges of the 
58 


rectangle. The four edges then form a closed loop which defines our domain as a 
59 


surface. We start by inputting the variables we need to construct the model. 
60 
ahallam 
2801 
\begin{python} 
61 
gross 
2948 
width=5000.0*m #width of model 
62 


depth=6000.0*m #depth of model 
63 
ahallam 
2801 
\end{python} 
64 
ahallam 
3153 

65 


\begin{figure}[ht] 
66 


\centerline{\includegraphics[width=4.in]{figures/pycadrec}} 
67 


\caption{Example 4: Rectangular Domain for \pycad} 
68 


\label{fig:pycad rec} 
69 


\end{figure} 
70 



71 
ahallam 
2979 
The variables are then used to construct the four corners of our domain, which 
72 
caltinay 
2982 
from the origin has the dimensions of $5000$ meters width and $6000$ meters 
73 


depth. This is done with the \verbPoint() function which accepts x, y and z 
74 


coordinates. Our domain is in two dimensions so z should always be zero. 
75 
ahallam 
2801 
\begin{python} 
76 
ahallam 
2597 
# Overall Domain 
77 
ahallam 
2634 
p0=Point(0.0, 0.0, 0.0) 
78 


p1=Point(0.0, depth, 0.0) 
79 


p2=Point(width, depth, 0.0) 
80 


p3=Point(width, 0.0, 0.0) 
81 
ahallam 
2801 
\end{python} 
82 
ahallam 
2979 
Now lines are defined using our points. This forms a rectangle around our 
83 
caltinay 
2982 
domain: 
84 
ahallam 
2801 
\begin{python} 
85 
ahallam 
2597 
l01=Line(p0, p1) 
86 


l12=Line(p1, p2) 
87 


l23=Line(p2, p3) 
88 


l30=Line(p3, p0) 
89 
ahallam 
2801 
\end{python} 
90 
caltinay 
2982 
Note that lines have a direction. These lines form the basis for our domain 
91 
ahallam 
2979 
boundary, which is a closed loop. 
92 
ahallam 
2801 
\begin{python} 
93 
ahallam 
2597 
c=CurveLoop(l01, l12, l23, l30) 
94 
ahallam 
2801 
\end{python} 
95 
ahallam 
2979 
Be careful to define the curved loop in an \textbf{anticlockwise} manner 
96 


otherwise the meshing algorithm may fail. 
97 
gross 
2948 
Finally we can define the domain as 
98 
ahallam 
2801 
\begin{python} 
99 
gross 
2948 
rec = PlaneSurface(c) 
100 


\end{python} 
101 
ahallam 
2979 
At this point the introduction of the curved loop seems to be unnecessary but 
102 
caltinay 
2982 
this concept plays an important role if holes are introduced. 
103 
gross 
2948 

104 
caltinay 
2982 
Now we are ready to hand over the domain \verbrec to a mesher which 
105 


subdivides the domain into triangles (or tetrahedra in 3D). In our case we use 
106 


\gmsh. We create an instance of the \verbDesign class which will handle the 
107 


interface to \gmsh: 
108 
gross 
2948 
\begin{python} 
109 


from esys.pycad.gmsh import Design 
110 


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


\end{python} 
112 
ahallam 
2979 
The argument \verbdim defines the spatial dimension of the domain\footnote{If 
113 


\texttt{dim}=3 the rectangle would be interpreted as a surface in the three 
114 
caltinay 
2982 
dimensional space.}. The second argument \verbelement_size defines the element 
115 
ahallam 
2979 
size which is the maximum length of a triangle edge in the mesh. The element 
116 


size needs to be chosen with care in order to avoid very dense meshes. If the 
117 


mesh is too dense, the computational time will be long but if the mesh is too 
118 


sparse, the modelled result will be poor. In our case with an element size of 
119 
caltinay 
2982 
$200$m and a domain length of $6000$m we will end up with about $\frac{6000m}{200m}=30$ 
120 
ahallam 
2979 
triangles in each spatial direction. So we end up with about $30 \times 30 = 
121 


900$ triangles which is a size that can be handled easily. 
122 
caltinay 
2982 
The domain \verbrec can simply be added to the \verbDesign; 
123 
gross 
2948 
\begin{python} 
124 


d.addItem(rec) 
125 


\end{python} 
126 
caltinay 
2982 
We have the plan to set a heat flux on the bottom of the domain. One can use 
127 


the masking technique to do this but \pycad offers a more convenient technique 
128 


called tagging. With this technique items in the domain are named using the 
129 


\verbPropertySet class. We can then later use this name to set values 
130 


specifically for those sample points located on the named items. Here we name 
131 
gross 
3644 
the bottom face of the domain where we will set the heat influx 
132 


\footnote{In some applications, eg. 
133 


when dealing with influxes, 
134 


it is required to have the surface beeing meshed without the need of explicitly 
135 


name the surface. In this case the line forming the surface of the domain need to be 
136 


added to the \texttt{Design} 
137 


using \texttt{d.addItem(rec, l01, l12, l23, l30)}}: 
138 
gross 
2948 
\begin{python} 
139 


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


d.addItem(ps) 
141 


\end{python} 
142 


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


\begin{python} 
144 


from esys.finley import MakeDomain 
145 


domain=MakeDomain(d) 
146 


\end{python} 
147 
ahallam 
2979 
The \verbdomain object can now be used in the same way like the return object 
148 
caltinay 
2982 
of the \verbRectangle object we have used previously to generate a mesh. It 
149 


is common practice to separate the mesh generation from the PDE solution. 
150 


The main reason for this is that mesh generation can be computationally very 
151 


expensive in particular in 3D. So it is more efficient to generate the mesh 
152 


once and write it to a file. The mesh can then be read in every time a new 
153 


simulation is run. \FINLEY supports this in the following 
154 


way\footnote{An alternative is using the \texttt{dump} and \texttt{load} 
155 


functions. They work with a binary format and tend to be much smaller.}: 
156 
gross 
2948 
\begin{python} 
157 
caltinay 
2982 
# write domain to a text file 
158 
gross 
2951 
domain.write("example04.fly") 
159 
gross 
2948 
\end{python} 
160 
caltinay 
2982 
and then for reading in another script: 
161 
gross 
2948 
\begin{python} 
162 


# read domain from text file 
163 


from esys.finley import ReadMesh 
164 
gross 
2951 
domain =ReadMesh("example04.fly") 
165 
gross 
2948 
\end{python} 
166 



167 
caltinay 
2982 
Before we discuss how to solve the PDE for this problem, it is useful to 
168 


present two additional options of the \verbDesign class. 
169 


These allow the user to access the script which is used by \gmsh to generate 
170 


the mesh as well as the generated mesh itself. This is done by setting specific 
171 
ahallam 
2979 
names for these files: 
172 
gross 
2948 
\begin{python} 
173 
gross 
2951 
d.setScriptFileName("example04.geo") 
174 


d.setMeshFileName("example04.msh") 
175 
gross 
2948 
\end{python} 
176 
caltinay 
2982 
Conventionally the extension \texttt{geo} is used for the script file of the 
177 


\gmsh geometry and the extension \texttt{msh} for the mesh file. Normally these 
178 


files are deleted after usage. 
179 
ahallam 
2979 
Accessing these files can be helpful to debug the generation of more complex 
180 


geometries. The geometry and the mesh can be visualised from the command line 
181 


using 
182 
gross 
2948 
\begin{verbatim} 
183 
gross 
2951 
gmsh example04.geo # show geometry 
184 


gmsh example04.msh # show mesh 
185 
gross 
2948 
\end{verbatim} 
186 


The mesh is shown in \reffig{fig:pycad rec mesh}. 
187 
ahallam 
3153 

188 
gross 
2948 
\begin{figure}[ht] 
189 
ahallam 
3153 
\centerline{\includegraphics[width=4.in]{figures/simplemesh}} 
190 


\caption{Example 4a: Mesh over rectangular domain, see \reffig{fig:pycad rec}} 
191 


\label{fig:pycad rec mesh} 
192 
gross 
2948 
\end{figure} 
193 
ahallam 
3153 
\clearpage 
194 
gross 
2948 

195 


\section{The Steadystate Heat Equation} 
196 
gross 
2951 
\sslist{example04b.py, cblib} 
197 
ahallam 
2979 
A temperature equilibrium or steady state is reached when the temperature 
198 
caltinay 
2982 
distribution in the model does not change with time. To calculate the steady 
199 


state solution the time derivative term in \refEq{eqn:Tform nabla} needs to be 
200 


set to zero; 
201 
gross 
2948 
\begin{equation}\label{eqn:Tform nabla steady} 
202 
jfenwick 
3308 
\nabla \cdot \kappa \nabla T = q_H 
203 
gross 
2948 
\end{equation} 
204 
ahallam 
2979 
This PDE is easier to solve than the PDE in \refEq{eqn:hdgenf2}, as no time 
205 


steps (iterations) are required. The \verbD term from \refEq{eqn:hdgenf2} is 
206 


simply dropped in this case. 
207 
gross 
2948 
\begin{python} 
208 


mypde=LinearPDE(domain) 
209 


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


\end{python} 
211 
caltinay 
2982 
The temperature at the top face of the domain is known as \verbTtop~($=20 C$). 
212 


In \refSec{Sec:1DHDv0} we have already discussed how this constraint is added 
213 


to the PDE: 
214 
gross 
2948 
\begin{python} 
215 


x=Solution(domain).getX() 
216 


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


\end{python} 
218 
ahallam 
2979 
Notice that we use the \verbsup function to calculate the maximum of $y$ 
219 


coordinates of the relevant sample points. 
220 
gross 
2948 

221 
ahallam 
2979 
In all cases so far we have assumed that the domain is insulated which 
222 
caltinay 
2982 
translates into a zero normal flux $n \cdot \kappa \nabla T$, see 
223 


\refEq{eq:hom flux}. In the modelling setup of this chapter we want to set 
224 


the normal heat flux at the bottom to \verbqin while still maintaining 
225 


insulation at the left and right face. Mathematically we can express this as 
226 
gross 
2948 
\begin{equation} 
227 
jfenwick 
3308 
n \cdot \kappa \nabla T = q_{S} 
228 
gross 
2948 
\label{eq:inhom flux} 
229 


\end{equation} 
230 
jfenwick 
3308 
where $q_{S}$ is a function of its location on the boundary. Its value 
231 
caltinay 
2982 
becomes zero for locations on the left or right face of the domain while it has 
232 


the value \verbqin at the bottom face. 
233 
jfenwick 
3308 
Notice that the value of $q_{S}$ at the top face is not relevant as we 
234 
ahallam 
2979 
prescribe the temperature here. 
235 
jfenwick 
3308 
We could define $q_{S}$ by using the masking techniques demonstrated 
236 
ahallam 
2979 
earlier. The tagging mechanism provides an alternative and in many cases more 
237 
caltinay 
2982 
convenient way of defining piecewise constant functions such as 
238 
jfenwick 
3308 
$q_{S}$. Recall now that the bottom face was denoted with the name 
239 
caltinay 
2982 
\verblinebottom when we defined the domain. 
240 
jfenwick 
3308 
We can use this now to create $q_{S}$; 
241 
gross 
2948 
\begin{python} 
242 


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


qS.setTaggedValue("linebottom",qin) 
244 


\end{python} 
245 
ahallam 
2979 
In the first line \verbqS is defined as a scalar value over the sample points 
246 
caltinay 
2982 
on the boundary of the domain. It is initialised to zero for all sample points. 
247 


In the second statement the values for those sample points which are located on 
248 


the line marked by \verblinebottom are set to \verbqin. 
249 
gross 
2948 

250 
caltinay 
2982 
The Neumann boundary condition assumed by \esc has the form 
251 
gross 
2948 
\begin{equation}\label{NEUMAN 2b} 
252 


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


\end{equation} 
254 
ahallam 
2979 
In comparison to the version in \refEq{NEUMAN 2} we have used so far the right 
255 
caltinay 
2982 
hand side is now the new PDE coefficient $y$. As we have not specified $y$ in 
256 


our previous examples, \esc has assumed the value zero for $y$. A comparison of 
257 


\refEq{NEUMAN 2b} and \refEq{eq:inhom flux} reveals that one needs to choose 
258 
jfenwick 
3308 
$y=q_{S}$; 
259 
gross 
2948 
\begin{python} 
260 


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


qS.setTaggedValue("linebottom",qin) 
262 


mypde.setValue(y=qS) 
263 


\end{python} 
264 
caltinay 
2982 
To plot the results we use the \modmpl library as shown in 
265 


\refSec{Sec:2DHD plot}. For convenience the interpolation of the temperature to 
266 


a rectangular grid for contour plotting is made available via the 
267 


\verbtoRegGrid function in the \verbcblib module. Your result should look 
268 


similar to \reffig{fig:steady state heat}. 
269 
gross 
2948 

270 


\begin{figure}[ht] 
271 
ahallam 
3153 
\centerline{\includegraphics[width=4.in]{figures/simpleheat}} 
272 


\caption{Example 4b: Result of simple steady state heat problem} 
273 


\label{fig:steady state heat} 
274 
gross 
2948 
\end{figure} 
275 
ahallam 
3153 
\clearpage 
276 
gross 
2948 
