14 
\section{Steadystate Heat Refraction} 
\section{Steadystate Heat Refraction} 
15 
\label{STEADYSTATE HEAT REFRACTION} 
\label{STEADYSTATE HEAT REFRACTION} 
16 


17 
In this chapter we will learn how to handle more complex geometries. 
In this chapter we show how to handle more complex geometries. 
18 


19 
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 models. In conjunction with \gmsh we can generate finite element meshes that conform to our domain's shape providing accurate modeling 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. 
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 


21 
We will introduce in two stages. First we will look at a very simple geometry. In fact, we will solve 
We proceed in this chapter by first looking at a very simple geometry. In fact, we solve 
22 
the steadystate heat equation over a rectangular domain. This is not a very interesting problem but in a second 
the steadystate heat equation over a rectangular domain. This is not a very interesting problem but then we extend our script by introducing an internal, curved interface. 

step we will extend our script by introducing an internal, curved interface. 

23 


24 
\begin{figure}[ht] 
\begin{figure}[ht] 
25 
\centerline{\includegraphics[width=4.in]{figures/pycadrec}} 
\centerline{\includegraphics[width=4.in]{figures/pycadrec}} 
27 
\label{fig:pycad rec} 
\label{fig:pycad rec} 
28 
\end{figure} 
\end{figure} 
29 


30 
\section{Example 4: Creating the model with \pycad} 
\section{Example 4: Creating the domain with \pycad} 
31 
\sslist{example04a.py} 
\sslist{example04a.py} 
32 


33 
We will modify the example in Chapter~\ref{CHAP HEAT 2a} in two ways: firstly we will look the steady state 
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 secondly we will use a more flexible tool 
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. 
to generate to generate the geometry. Lets look at the geometry first. 
36 


37 
We want to define a rectangular domain of width $5 km$ and depth $6 km$ below the surface of the Earth. Under the assumption of a known temperature at the surface, a known heat flux at the bottom and 
We want to define a rectangular domain of width $5 km$ and depth $6 km$ below the surface of the Earth. Under the assumption of a known temperature at the surface, a known heat flux at the bottom and 
38 
insulation to both sides we want to calculate the steadystate temperature distribution. 
insulation to both sides we want to calculate the steadystate temperature distribution. 
39 


40 
In \pycad there are a few primary constructors that build upon each other to define domains and boundaries; 
In \pycad there are a few primary constructors that build upon each other to define domains and boundaries; 
41 
the ones we will use are: 
the ones we use are: 
42 
\begin{python} 
\begin{python} 
43 
from esys.pycad import * 
from esys.pycad import * 
44 
Point() #Create a point in space. 
Point() #Create a point in space. 
87 
from esys.pycad.gmsh import Design 
from esys.pycad.gmsh import Design 
88 
d=Design(dim=2, element_size=200*m) 
d=Design(dim=2, element_size=200*m) 
89 
\end{python} 
\end{python} 
90 
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 choose with care in order to avoid very large meshes if you don't want to. In our case with an element size of $200$m 
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 large meshes if you don't want to. In our case with an element size of $200$m 
91 
and a domain length of $6000$m we will end up with about $\frac{6000m}{200m}=30$ triangles in each spatial direction. So will end up with about $30 \times 30 = 900$ triangles which is a size that can be handled easily. 
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. 
92 
We can easily add our domain \verbrec to the \verbDesign; 
We can easily add our domain \verbrec to the \verbDesign; 
93 
\begin{python} 
\begin{python} 
94 
d.addItem(rec) 
d.addItem(rec) 
95 
\end{python} 
\end{python} 
96 
We have the plan to set a heat flux on the bottom of the domain. One can use the masking technique to do this 
We have the plan to set a heat flux on the bottom of the domain. One can use the masking technique to do this 
97 
but \pycad offers a more convenient technique called tagging. With this techniques items in the domain are 
but \pycad offers a more convenient technique called tagging. With this technique items in the domain are 
98 
named using the \verbPropertySet class. We can then later use this name to set values secifically for 
named using the \verbPropertySet class. We can then later use this name to set values secifically for 
99 
those sample points located on the named items. Here we name the bottom face of the 
those sample points located on the named items. Here we name the bottom face of the 
100 
domain where we will set the heat influx: 
domain where we will set the heat influx: 
116 
# write domain to an text file 
# write domain to an text file 
117 
domain.write("example04.fly") 
domain.write("example04.fly") 
118 
\end{python} 
\end{python} 
119 
and then for recovery in another script; 
and then for reading in another script; 
120 
\begin{python} 
\begin{python} 
121 
# read domain from text file 
# read domain from text file 
122 
from esys.finley import ReadMesh 
from esys.finley import ReadMesh 
131 


132 
Before we discuss how we solve the PDE for the 
Before we discuss how we solve the PDE for the 
133 
problem it is useful to present two additional options of the \verbDesign. 
problem it is useful to present two additional options of the \verbDesign. 
134 
They allow the user accessing the script which is used to by \gmsh to generate the mesh as well as 
They allow the user accessing the script which is used by \gmsh to generate the mesh as well as 
135 
the mesh as it has been generated by \gmsh. This is done by setting specific names for these files: 
the mesh as it has been generated by \gmsh. This is done by setting specific names for these files: 
136 
\begin{python} 
\begin{python} 
137 
d.setScriptFileName("example04.geo") 
d.setScriptFileName("example04.geo") 
139 
\end{python} 
\end{python} 
140 
Usually the extension \texttt{geo} is used for the script file of the \gmsh geometry and 
Usually the extension \texttt{geo} is used for the script file of the \gmsh geometry and 
141 
the extension \texttt{msh} for the mesh file. Normally these files are deleted after usage. 
the extension \texttt{msh} for the mesh file. Normally these files are deleted after usage. 
142 
Accessing these file can be helpful to debug the generation of more complex geometries. The geometry and the mesher can be visualise with \gmsh from the command line using 
Accessing these file can be helpful to debug the generation of more complex geometries. The geometry and the mesh can be visualised from the command line using 
143 
\begin{verbatim} 
\begin{verbatim} 
144 
gmsh example04.geo # show geometry 
gmsh example04.geo # show geometry 
145 
gmsh example04.msh # show mesh 
gmsh example04.msh # show mesh 
155 
\section{The Steadystate Heat Equation} 
\section{The Steadystate Heat Equation} 
156 
\sslist{example04b.py, cblib} 
\sslist{example04b.py, cblib} 
157 


158 
Steady state is reached in the temperature is not changing in time. So to calculate the 
Steady state is reached in the temperature when it is not changing in time. So to calculate the 
159 
the steady state we set the time derivative term in \refEq{eqn:Tform nabla} to zero; 
the steady state we set the time derivative term in \refEq{eqn:Tform nabla} to zero; 
160 
\begin{equation}\label{eqn:Tform nabla steady} 
\begin{equation}\label{eqn:Tform nabla steady} 
161 
\nabla \cdot \kappa \nabla T = q\hackscore H 
\nabla \cdot \kappa \nabla T = q\hackscore H 
162 
\end{equation} 
\end{equation} 
163 
This PDE is easier then the PDE in \refEq{eqn:hdgenf2} we needed to solve in each time step. We can just drop 
This PDE is easier to solve, than the PDE in \refEq{eqn:hdgenf2} in each time step. We can just drop 
164 
the \verbD term; 
the \verbD term; 
165 
\begin{python} 
\begin{python} 
166 
mypde=LinearPDE(domain) 
mypde=LinearPDE(domain) 
172 
x=Solution(domain).getX() 
x=Solution(domain).getX() 
173 
mypde.setValue(q=whereZero(x[1]sup(x[1])),r=Ttop) 
mypde.setValue(q=whereZero(x[1]sup(x[1])),r=Ttop) 
174 
\end{python} 
\end{python} 
175 
Notice that we have use the \verbsup function to calculate the maximum $y$ coordinator of the relevant sample points. 
Notice that we use the \verbsup function to calculate the maximum of $y$ coordinates of the relevant sample points. 
176 


177 
In all cases so far we have assumed that the domain is insulated which translates 
In all cases so far we have assumed that the domain is insulated which translates 
178 
into a zero normal flux $n \cdot \kappa \nabla T$, see \refEq{eq:hom flux}. In the modeling 
into a zero normal flux $n \cdot \kappa \nabla T$, see \refEq{eq:hom flux}. In the modeling 
185 
where $q\hackscore{S}$ is a function of its location on the boundary. Its value becomes zero 
where $q\hackscore{S}$ is a function of its location on the boundary. Its value becomes zero 
186 
for locations on the left or right face of the domain while it has the value \verbqin at the bottom face. 
for locations on the left or right face of the domain while it has the value \verbqin at the bottom face. 
187 
Notice that the value of $q\hackscore{S}$ at the top face is not relevant as we prescribe the temperature here. 
Notice that the value of $q\hackscore{S}$ at the top face is not relevant as we prescribe the temperature here. 
188 
One could use the masking techniques we have already extensively used to define $q\hackscore{S}$ in \esc 
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 

but the tagging mechanism provides an alternative and in many cases more convenient way of defining piecewise 

189 
constant functions such as $q\hackscore{S}$. You need to remember now that we 
constant functions such as $q\hackscore{S}$. You need to remember now that we 
190 
marked the bottom face with the name \verblinebottom when we defined the domain. 
marked the bottom face with the name \verblinebottom when we defined the domain. 
191 
We can use this now to create $q\hackscore{S}$; 
We can use this now to create $q\hackscore{S}$; 
233 
Our 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 in fact warped by the structure of the model and is heavily dependent upon the material properties and shape. 
Our 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 in fact warped by the structure of the model and is heavily dependent upon the material properties and shape. 
234 


235 


236 
We will modify the example of the previous section by subdividing the block into two parts. The curve 
We modify the example of the previous section by subdividing the block into two parts. The curve 
237 
separating the two blocks is given as a spline, see \reffig{fig:anticlinehrmodel}. The data points 
separating the two blocks is given as a spline, see \reffig{fig:anticlinehrmodel}. The data points 
238 
used to define the curve may be measurement but for simplicity we assume here that there coordinates are 
used to define the curve may be measurement but for simplicity we assume here that there coordinates are 
239 
known in analytic form. 
known in analytic form. 
242 


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


245 
The curve defining our clinal structure is located in approximately 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. 
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. 
246 
\begin{python} 
\begin{python} 
247 
# Material Boundary 
# Material Boundary 
248 
x=[ Point(i*dsp\ 
x=[ Point(i*dsp\ 
275 
#clockwise check 
#clockwise check 
276 
bblockloop=CurveLoop(mysp,Line(x2,p2),Line(p2,p1),Line(p1,x1)) 
bblockloop=CurveLoop(mysp,Line(x2,p2),Line(p2,p1),Line(p1,x1)) 
277 
\end{python} 
\end{python} 
278 
The last few steps in creating the model take the previously defined domain and subdomain points and generate a mesh that can be imported into \esc. 
The last few steps in creating the model are: take the previously defined domain and subdomain points and generate a mesh that can be imported into \esc. 
279 
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. 
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. 
280 
\begin{python} 
\begin{python} 
281 
# Create a Design which can make the mesh 
# Create a Design which can make the mesh 
313 
\sslist{example05b.py and cblib.py} 
\sslist{example05b.py and cblib.py} 
314 


315 
We want to investigate the profile of the data of the last example. 
We want to investigate the profile of the data of the last example. 
316 
We are particularly interested in the deepth profile of the heat flux which is 
We are particularly interested in the depth profile of the heat flux which is 
317 
the second component of $\kappa \nabla T$. We will extend the script developed in the 
the second component of $\kappa \nabla T$. We extend the script developed in the 
318 
previous section to show how for instance vertical profile can be plotted. 
previous section to show how for instance vertical profile can be plotted. 
319 


320 
The first important information is that \esc assumes that $\kappa \nabla T$ is not smooth and 
The first important information is that \esc assumes that $\kappa \nabla T$ is not smooth and 
322 
the flux is the product of the piecewise constant function $\kappa$ and 
the flux is the product of the piecewise constant function $\kappa$ and 
323 
the gradient of the temperature $T$ which has a kink across the rock interface. 
the gradient of the temperature $T$ which has a kink across the rock interface. 
324 
Before plotting this function we need to smooth it using the 
Before plotting this function we need to smooth it using the 
325 
\verbProjector() factory; 
\verbProjector() class; 
326 
\begin{python} 
\begin{python} 
327 
from esys.escript.pdetools import Projector 
from esys.escript.pdetools import Projector 
328 
proj=Projector(domain) 
proj=Projector(domain) 
335 
\begin{python} 
\begin{python} 
336 
xiq,yiq,ziq = toRegGrid(qu[1],200,200) 
xiq,yiq,ziq = toRegGrid(qu[1],200,200) 
337 
\end{python} 
\end{python} 
338 
using the \verbtoRegGrid function from the cookbook library which we are already using for the contour plot. 
using the \verbtoRegGrid function from the cookbook library which we are using for the contour plot. 
339 
At return \verbziq[j,i] is the value of vertical heat flux at point 
At return \verbziq[j,i] is the value of vertical heat flux at point 
340 
(\verbxiq[i],\verbyiq[j]). We can easily create deep profiles now by 
(\verbxiq[i],\verbyiq[j]). We can easily create deep profiles now by 
341 
plotting slices \verbziq[:,i] over \verbyiq. The following script 
plotting slices \verbziq[:,i] over \verbyiq. The following script 
358 


359 
\section{Arrow plots in \mpl} 
\section{Arrow plots in \mpl} 
360 
\sslist{example05c.py and cblib.py} 
\sslist{example05c.py and cblib.py} 
361 
We would like to visualize the distribution of the flux $\kappa \nabla T$ over the domain 
We would like to visualise the distribution of the flux $\kappa \nabla T$ over the domain 
362 
and produce a plot like shown in \reffig{fig:hr001qumodel}. 
and produce a plot like shown in \reffig{fig:hr001qumodel}. 
363 
The plot puts together three components. A contour plot of the temperature, 
The plot puts together three components. A contour plot of the temperature, 
364 
a colored representation of the two subdomains where colour represents the thermal conductivity 
a colored representation of the two subdomains where colour represents the thermal conductivity 