1 
ahallam 
2401 

2 


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


% 
4 
jfenwick 
2881 
% Copyright (c) 20032010 by University of Queensland 
5 
ahallam 
2401 
% 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 
2905 
\begin{figure}[h!] 
15 


\centerline{\includegraphics[width=4.in]{figures/onedheatdiff001}} 
16 
ahallam 
2979 
\caption{Example 1: Temperature differential along a single interface between 
17 


two granite blocks.} 
18 
gross 
2905 
\label{fig:onedgbmodel} 
19 


\end{figure} 
20 
ahallam 
2801 

21 
gross 
2949 
\section{Example 1: One Dimensional Heat Diffusion in Granite} 
22 
gross 
2905 
\label{Sec:1DHDv00} 
23 
gross 
2878 

24 
ahallam 
2979 
The first model consists of two blocks of isotropic material, for instance 
25 


granite, sitting next to each other. 
26 


Initial temperature in \textit{Block 1} is \verbT1 and in \textit{Block 2} is 
27 


\verbT2. 
28 
gross 
2905 
We assume that the system is insulated. 
29 


What would happen to the temperature distribution in each block over time? 
30 
ahallam 
2979 
Intuition tells us that heat will be transported from the hotter block to the 
31 


cooler one until both 
32 
gross 
2905 
blocks have the same temperature. 
33 



34 
ahallam 
2495 
\subsection{1D Heat Diffusion Equation} 
35 
ahallam 
2979 
We can model the heat distribution of this problem over time using one 
36 


dimensional heat diffusion equation\footnote{A detailed discussion on how the 
37 


heat diffusion equation is derived can be found at 
38 


\url{ 
39 


http://online.redwoods.edu/instruct/darnold/DEProj/sp02/AbeRichards/paper.pdf}}; 
40 
ahallam 
2494 
which is defined as: 
41 
ahallam 
2401 
\begin{equation} 
42 
ahallam 
2979 
\rho c\hackscore p \frac{\partial T}{\partial t}  \kappa \frac{\partial^{2} 
43 


T}{\partial x^{2}} = q\hackscore H 
44 
ahallam 
2401 
\label{eqn:hd} 
45 


\end{equation} 
46 
ahallam 
2979 
where $\rho$ is the material density, $c\hackscore p$ is the specific heat and 
47 


$\kappa$ is the thermal 
48 


conductivity\footnote{A list of some common thermal conductivities is available 
49 


from Wikipedia 
50 


\url{http://en.wikipedia.org/wiki/List_of_thermal_conductivities}}. Here we 
51 


assume that these material 
52 
gross 
2861 
parameters are \textbf{constant}. 
53 
ahallam 
2979 
The heat source is defined by the right hand side of \refEq{eqn:hd} as 
54 


$q\hackscore{H}$; this can take the form of a constant or a function of time and 
55 


space. For example $q\hackscore{H} = q\hackscore{0}e^{\gamma t}$ where we have 
56 


the output of our heat source decaying with time. There are also two partial 
57 


derivatives in \refEq{eqn:hd}; $\frac{\partial T}{\partial t}$ describes the 
58 


change in temperature with time while $\frac{\partial ^2 T}{\partial x^2}$ is 
59 


the spatial change of temperature. As there is only a single spatial dimension 
60 


to our problem, our temperature solution $T$ is only dependent on the time $t$ 
61 


and our signed distance from the the blockblock interface $x$. 
62 
ahallam 
2401 

63 
gross 
2861 
\subsection{PDEs and the General Form} 
64 
ahallam 
2979 
It is possible to solve PDE \refEq{eqn:hd} analytically and obtain an exact 
65 


solution to our problem. However, it is not always practical to solve the 
66 


problem this way. Alternatively, computers can be used to find the solution. To 
67 


do this, a numerical approach is required to discretise 
68 


the PDE \refEq{eqn:hd} across time and space, this reduces the problem to a 
69 


finite number of equations for a finite number of spatial points and time steps. 
70 


These parameters together define the model. While discretisation introduces 
71 


approximations and a degree of error, a sufficiently sampled model is generally 
72 


accurate enough to satisfy the accuracy requirements for the final solution. 
73 
gross 
2477 

74 
ahallam 
2979 
Firstly, we discretise the PDE \refEq{eqn:hd} in time. This leaves us with a 
75 


steady linear PDE which involves spatial derivatives only and needs to be solved 
76 


in each time step to progress in time. \esc can help us here. 
77 
gross 
2861 

78 
ahallam 
2979 
For time discretization we use the Backwards Euler approximation 
79 


scheme\footnote{see \url{http://en.wikipedia.org/wiki/Euler_method}}. It is 
80 


based on the 
81 
gross 
2861 
approximation 
82 


\begin{equation} 
83 


\frac{\partial T(t)}{\partial t} \approx \frac{T(t)T(th)}{h} 
84 


\label{eqn:beuler} 
85 


\end{equation} 
86 


for $\frac{\partial T}{\partial t}$ at time $t$ 
87 


where $h$ is the time step size. This can also be written as; 
88 


\begin{equation} 
89 


\frac{\partial T}{\partial t}(t^{(n)}) \approx \frac{T^{(n)}  T^{(n1)}}{h} 
90 


\label{eqn:Tbeuler} 
91 


\end{equation} 
92 
ahallam 
2979 
where the upper index $n$ denotes the n\textsuperscript{th} time step. So one 
93 


has 
94 
gross 
2861 
\begin{equation} 
95 


\begin{array}{rcl} 
96 


t^{(n)} & = & t^{(n1)}+h \\ 
97 


T^{(n)} & = & T(t^{(n1)}) \\ 
98 


\end{array} 
99 


\label{eqn:Neuler} 
100 


\end{equation} 
101 


Substituting \refEq{eqn:Tbeuler} into \refEq{eqn:hd} we get; 
102 


\begin{equation} 
103 
ahallam 
2979 
\frac{\rho c\hackscore p}{h} (T^{(n)}  T^{(n1)})  \kappa \frac{\partial^{2} 
104 


T^{(n)}}{\partial x^{2}} = q\hackscore H 
105 
gross 
2861 
\label{eqn:hddisc} 
106 


\end{equation} 
107 
ahallam 
2979 
Notice that we evaluate the spatial derivative term at the current time 
108 


$t^{(n)}$  therefore the name \textbf{backward Euler} scheme. Alternatively, 
109 


one can evaluate the spatial derivative term at the previous time $t^{(n1)}$. 
110 


This 
111 


approach is called the \textbf{forward Euler} scheme. This scheme can provide 
112 


some computational advantages, which 
113 


are not discussed here. However, the \textbf{forward Euler} scheme has a major 
114 


disadvantage. Namely, depending on the 
115 


material parameters as well as the domain discretization of the spatial 
116 


derivative term, the time step size $h$ needs to be chosen sufficiently small to 
117 


achieve a stable temperature when progressing in time. Stabiliy is achieved if 
118 


the temperature does not grow beyond its initial bounds and become 
119 


nonphysical. 
120 


The backward Euler scheme, which we use here, is unconditionally stable meaning 
121 


that under the assumption of a 
122 


physically correct problem setup the temperature approximation remains physical 
123 


for all time steps. 
124 


The user needs to keep in mind that the discretization error introduced by 
125 


\refEq{eqn:beuler} 
126 


is sufficiently small, thus a good approximation of the true temperature is 
127 


computed. It is 
128 


therefore very important that any results are viewed with caution. For example, 
129 


one may compare the results for different time and spatial step sizes. 
130 
gross 
2861 

131 


To get the temperature $T^{(n)}$ at time $t^{(n)}$ we need to solve the linear 
132 
ahallam 
2979 
differential equation \refEq{eqn:hddisc} which only includes spatial 
133 


derivatives. To solve this problem 
134 
artak 
2958 
we want to use \esc. 
135 
gross 
2861 

136 
ahallam 
2979 
In \esc any given PDE can be described by the general form. For the purpose of 
137 


this introduction we illustrate a simpler version of the general form for full 
138 


linear PDEs which is available in the \esc user's guide. A simplified form that 
139 


suits our heat diffusion problem\footnote{The form in the \esc users guide which 
140 


uses the Einstein convention is written as 
141 
gross 
2477 
$(A\hackscore{jl} u\hackscore{,l})\hackscore{,j}+D u =Y$} 
142 
ahallam 
2606 
is described by; 
143 
gross 
2477 
\begin{equation}\label{eqn:commonform nabla} 
144 
jfenwick 
2657 
\nabla\cdot(A\cdot\nabla u) + Du = f 
145 
ahallam 
2411 
\end{equation} 
146 
ahallam 
2979 
where $A$, $D$ and $f$ are known values and $u$ is the unknown solution. The 
147 


symbol $\nabla$ which is called the \textit{Nabla operator} or \textit{del 
148 


operator} represents 
149 


the spatial derivative of its subject  in this case $u$. Lets assume for a 
150 


moment that we deal with a onedimensional problem then ; 
151 
gross 
2477 
\begin{equation} 
152 


\nabla = \frac{\partial}{\partial x} 
153 


\end{equation} 
154 
gross 
2861 
and we can write \refEq{eqn:commonform nabla} as; 
155 
ahallam 
2411 
\begin{equation}\label{eqn:commonform} 
156 
gross 
2477 
A\frac{\partial^{2}u}{\partial x^{2}} + Du = f 
157 
ahallam 
2411 
\end{equation} 
158 
ahallam 
2979 
if $A$ is constant. To match this simplified general form to our problem 
159 


\refEq{eqn:hddisc} 
160 
gross 
2861 
we rearrange \refEq{eqn:hddisc}; 
161 
ahallam 
2411 
\begin{equation} 
162 
ahallam 
2979 
\frac{\rho c\hackscore p}{h} T^{(n)}  \kappa \frac{\partial^2 T^{(n)}}{\partial 
163 


x^2} = q\hackscore H + \frac{\rho c\hackscore p}{h} T^{(n1)} 
164 
ahallam 
2494 
\label{eqn:hdgenf} 
165 


\end{equation} 
166 
ahallam 
2979 
The PDE is now in a form that satisfies \refEq{eqn:commonform nabla} which is 
167 


required for \esc to solve our PDE. This can be done by generating a solution 
168 


for successive increments in the time nodes $t^{(n)}$ where 
169 


$t^{(0)}=0$ and $t^{(n)}=t^{(n1)}+h$ where $h>0$ is the step size and assumed 
170 


to be constant. 
171 


In the following the upper index ${(n)}$ refers to a value at time $t^{(n)}$. 
172 


Finally, by comparing \refEq{eqn:hdgenf} with \refEq{eqn:commonform} one can see 
173 


that; 
174 
gross 
2862 
\begin{equation}\label{ESCRIPT SET} 
175 
gross 
2861 
u=T^{(n)}; 
176 
ahallam 
2979 
A = \kappa; D = \frac{\rho c \hackscore{p}}{h}; f = q \hackscore{H} + \frac{\rho 
177 


c\hackscore p}{h} T^{(n1)} 
178 
ahallam 
2494 
\end{equation} 
179 



180 
ahallam 
2495 
\subsection{Boundary Conditions} 
181 
gross 
2878 
\label{SEC BOUNDARY COND} 
182 
ahallam 
2979 
With the PDE sufficiently modified, consideration must now be given to the 
183 


boundary conditions of our model. Typically there are two main types of boundary 
184 


conditions known as \textbf{Neumann} and \textbf{Dirichlet} boundary 
185 


conditions\footnote{More information on Boundary Conditions is available at 
186 


Wikipedia \url{http://en.wikipedia.org/wiki/Boundary_conditions}}, 
187 


respectively. 
188 


A \textbf{Dirichlet boundary condition} is conceptually simpler and is used to 
189 


prescribe a known value to the unknown solution (in our example the temperature) 
190 


on parts of the boundary or on the entire boundary of the region of interest. 
191 


We discuss Dirichlet boundary condition in our second example presented in 
192 


Section~\ref{Sec:1DHDv0}. 
193 
ahallam 
2495 

194 
ahallam 
2979 
However, for this example we have made the model assumption that the system is 
195 


insulated, so we need 
196 
gross 
2905 
to add an appropriate boundary condition to prevent 
197 
ahallam 
2979 
any loss or inflow of energy at the boundary of our domain. Mathematically this 
198 


is expressed by prescribing 
199 


the heat flux $\kappa \frac{\partial T}{\partial x}$ to zero. In our simplified 
200 


one dimensional model this is expressed 
201 
gross 
2862 
in the form; 
202 
ahallam 
2494 
\begin{equation} 
203 
gross 
2862 
\kappa \frac{\partial T}{\partial x} = 0 
204 
ahallam 
2494 
\end{equation} 
205 
gross 
2862 
or in a more general case as 
206 


\begin{equation}\label{NEUMAN 1} 
207 


\kappa \nabla T \cdot n = 0 
208 


\end{equation} 
209 
ahallam 
2979 
where $n$ is the outer normal field \index{outer normal field} at the surface 
210 


of the domain. 
211 


The $\cdot$ (dot) refers to the dot product of the vectors $\nabla T$ and $n$. 
212 


In fact, the term $\nabla T \cdot n$ is the normal derivative of 
213 


the temperature $T$. Other notations used here are\footnote{The \esc notation 
214 


for the normal 
215 
gross 
2862 
derivative is $T\hackscore{,i} n\hackscore i$.}; 
216 
ahallam 
2645 
\begin{equation} 
217 
gross 
2862 
\nabla T \cdot n = \frac{\partial T}{\partial n} \; . 
218 
ahallam 
2645 
\end{equation} 
219 
ahallam 
2979 
A condition of the type \refEq{NEUMAN 1} defines a \textbf{Neuman boundary 
220 


condition} for the PDE. 
221 
ahallam 
2494 

222 
gross 
2905 
The PDE \refEq{eqn:hdgenf} 
223 
ahallam 
2979 
and the Neuman boundary condition~\ref{eqn:hdgenf} (potentially together with 
224 


the Dirichlet boundary conditions) define a \textbf{boundary value problem}. 
225 


It is the nature of a boundary value problem to allow making statements about 
226 


the solution in the 
227 


interior of the domain from information known on the boundary only. In most 
228 


cases 
229 


we use the term partial differential equation but in fact it is a boundary value 
230 


problem. 
231 


It is important to keep in mind that boundary conditions need to be complete and 
232 


consistent in the sense that 
233 


at any point on the boundary either a Dirichlet or a Neuman boundary condition 
234 


must be set. 
235 
gross 
2862 

236 
ahallam 
2979 
Conveniently, \esc makes default assumption on the boundary conditions which the 
237 


user may modify where appropriate. 
238 


For a problem of the form in~\refEq{eqn:commonform nabla} the default 
239 


condition\footnote{In the form of the \esc users guide which is using the 
240 


Einstein convention is written as 
241 
gross 
2862 
$n\hackscore{j}A\hackscore{jl} u\hackscore{,l}=0$.} is; 
242 


\begin{equation}\label{NEUMAN 2} 
243 
gross 
2948 
n\cdot A \cdot\nabla u = 0 
244 
gross 
2862 
\end{equation} 
245 
ahallam 
2979 
which is used everywhere on the boundary. Again $n$ denotes the outer normal 
246 


field. 
247 


Notice that the coefficient $A$ is the same as in the \esc 
248 


PDE~\ref{eqn:commonform nabla}. 
249 


With the settings for the coefficients we have already identified in 
250 


\refEq{ESCRIPT SET} this 
251 
gross 
2862 
condition translates into 
252 
gross 
2867 
\begin{equation}\label{NEUMAN 2b} 
253 
gross 
2862 
\kappa \frac{\partial T}{\partial x} = 0 
254 


\end{equation} 
255 
ahallam 
2979 
for the boundary of the domain. This is identical to the Neuman boundary 
256 


condition we want to set. \esc will take care of this condition for us. We 
257 


discuss the Dirichlet boundary condition later. 
258 
gross 
2862 

259 
gross 
2870 
\subsection{Outline of the Implementation} 
260 


\label{sec:outline} 
261 
ahallam 
2979 
To solve the heat diffusion equation (\refEq{eqn:hd}) we write a simple \pyt 
262 


script. At this point we assume that you have some basic understanding of the 
263 


\pyt programming language. If not, there are some pointers and links available 
264 


in Section \ref{sec:escpybas}. The script (discussed in \refSec{sec:key}) has 
265 


four major steps. Firstly, we need to define the domain where we want to 
266 


calculate the temperature. For our problem this is the joint blocks of granite 
267 


which has a rectangular shape. Secondly, we need to define the PDE to solve in 
268 


each time step to get the updated temperature. Thirdly, we need to define the 
269 


coefficients of the PDE and finally we need to solve the PDE. The last two steps 
270 


need to be repeated until the final time marker has been reached. The work flow 
271 


described in \reffig{fig:wf}. 
272 
ahallam 
2975 
% \begin{enumerate} 
273 


% \item create domain 
274 


% \item create PDE 
275 


% \item while end time not reached: 
276 


% \begin{enumerate} 
277 


% \item set PDE coefficients 
278 


% \item solve PDE 
279 


% \item update time marker 
280 


% \end{enumerate} 
281 


% \item end of calculation 
282 


% \end{enumerate} 
283 



284 


\begin{figure}[h!] 
285 


\centering 
286 


\includegraphics[width=1in]{figures/workflow.png} 
287 


\caption{Workflow for developing an \esc model and solution.} 
288 


\label{fig:wf} 
289 


\end{figure} 
290 



291 
ahallam 
2979 
In the terminology of \pyt, the domain and PDE are represented by 
292 


\textbf{objects}. The nice feature of an object is that it is defined by its 
293 


usage and features 
294 


rather than its actual representation. So we will create a domain object to 
295 


describe the geometry of the two 
296 


granite blocks. Then we define PDEs and spatially distributed values such as the 
297 


temperature 
298 


on this domain. Similarly, to define a PDE object we use the fact that one needs 
299 


only to define the coefficients of the PDE and solve the PDE. The PDE object has 
300 


advanced features, but these are not required in simple cases. 
301 
gross 
2870 

302 



303 


\begin{figure}[t] 
304 


\centering 
305 


\includegraphics[width=6in]{figures/functionspace.pdf} 
306 
ahallam 
2975 
\caption{\esc domain construction overview} 
307 
gross 
2870 
\label{fig:fs} 
308 


\end{figure} 
309 



310 


\subsection{The Domain Constructor in \esc} 
311 


\label{ss:domcon} 
312 
ahallam 
2979 
Whilst it is not strictly relevant or necessary, a better understanding of 
313 


how values are spatially distributed (\textit{e.g.} Temperature) and how PDE 
314 


coefficients are interpreted in \esc can be helpful. 
315 
gross 
2870 

316 
ahallam 
2979 
There are various ways to construct domain objects. The simplest form is a 
317 


rectangular shaped region with a length and height. There is 
318 


a ready to use function for this named \verb rectangle(). Besides the spatial 
319 


dimensions this function requires to specify the number of 
320 


elements or cells to be used along the length and height, see \reffig{fig:fs}. 
321 


Any spatially distributed value 
322 


and the PDE is represented in discrete form using this element 
323 


representation\footnote{We use the finite element method (FEM), see 
324 


\url{http://en.wikipedia.org/wiki/Finite_element_method} for details.}. 
325 


Therefore we will have access to an approximation of the true PDE solution 
326 


only. 
327 


The quality of the approximation depends  besides other factors mainly on the 
328 


number of elements being used. In fact, the 
329 


approximation becomes better when more elements are used. However, computational 
330 


cost grows with the number of 
331 


elements being used. It is therefore important that you find the right balance 
332 


between the demand in accuracy and acceptable resource usage. 
333 
gross 
2870 

334 
ahallam 
2979 
In general, one can think about a domain object as a composition of nodes and 
335 


elements. 
336 


As shown in \reffig{fig:fs}, an element is defined by the nodes that are used to 
337 


describe its vertices. 
338 
gross 
2870 
To represent spatial distributed values the user can use 
339 
ahallam 
2979 
the values at the nodes, at the elements in the interior of the domain or at the 
340 


elements located at the surface of the domain. 
341 


The different approach used to represent values is called \textbf{function 
342 


space} and is attached to all objects 
343 


in \esc representing a spatial distributed value such as the solution of a PDE. 
344 


The three 
345 
artak 
2963 
function spaces we use at the moment are; 
346 
gross 
2870 
\begin{enumerate} 
347 


\item the nodes, called by \verbContinuousFunction(domain) ; 
348 


\item the elements/cells, called by \verbFunction(domain) ; and 
349 


\item the boundary, called by \verbFunctionOnBoundary(domain) . 
350 


\end{enumerate} 
351 
ahallam 
2979 
A function space object such as \verbContinuousFunction(domain) has the method 
352 


\verbgetX attached to it. This method returns the 
353 


location of the socalled \textbf{sample points} used to represent values of the 
354 


particular function space. So the 
355 


call \verbContinuousFunction(domain).getX() will return the coordinates of the 
356 


nodes used to describe the domain while 
357 


the \verbFunction(domain).getX() returns the coordinates of numerical 
358 


integration points within elements, see 
359 
gross 
2905 
\reffig{fig:fs}. 
360 
gross 
2870 

361 
ahallam 
2979 
This distinction between different representations of spatially distributed 
362 


values 
363 


is important in order to be able to vary the degrees of smoothness in a PDE 
364 


problem. 
365 


The coefficients of a PDE do not need to be continuous, thus this qualifies as a 
366 


\verbFunction() type. 
367 


On the other hand a temperature distribution must be continuous and needs to be 
368 


represented with a \verbContinuousFunction() function space. 
369 


An influx may only be defined at the boundary and is therefore a \verb 
370 


FunctionOnBoundary() object. 
371 


\esc allows certain transformations of the function spaces. A \verb 
372 


ContinuousFunction() can be transformed into a \verbFunctionOnBoundary() 
373 


or \verbFunction(). On the other hand there is not enough information in a 
374 


\verb FunctionOnBoundary() to transform it to a \verb ContinuousFunction() . 
375 


These transformations, which are called \textbf{interpolation} are invoked 
376 


automatically by \esc if needed. 
377 
gross 
2870 

378 
artak 
2963 
Later in this introduction we discuss how 
379 
ahallam 
2979 
to define specific areas of geometry with different materials which are 
380 


represented by different material coefficients such the 
381 


thermal conductivities $\kappa$. A very powerful technique to define these types 
382 


of PDE 
383 


coefficients is tagging. Blocks of materials and boundaries can be named and 
384 


values can be defined on subregions based on their names. 
385 


This is a method for simplifying PDE coefficient and flux definitions. It makes 
386 


scripting much easier and we will discuss this technique in 
387 


Section~\ref{STEADYSTATE HEAT REFRACTION}. 
388 
gross 
2870 

389 



390 


\subsection{A Clarification for the 1D Case} 
391 
gross 
2931 
\label{SEC: 1D CLARIFICATION} 
392 
ahallam 
2979 
It is necessary for clarification that we revisit our general PDE from 
393 


\refeq{eqn:commonform nabla} for two dimensional domain. \esc is inherently 
394 


designed to solve problems that are greater than one dimension and so 
395 


\refEq{eqn:commonform nabla} needs to be read as a higher dimensional problem. 
396 


In the case of two spatial dimensions the \textit{Nabla operator} has in fact 
397 


two components $\nabla = (\frac{\partial}{\partial x}, \frac{\partial}{\partial 
398 


y})$. Assuming the coefficient $A$ is constant, the \refEq{eqn:commonform nabla} 
399 


takes the following form; 
400 
gross 
2477 
\begin{equation}\label{eqn:commonform2D} 
401 


A\hackscore{00}\frac{\partial^{2}u}{\partial x^{2}} 
402 


A\hackscore{01}\frac{\partial^{2}u}{\partial x\partial y} 
403 


A\hackscore{10}\frac{\partial^{2}u}{\partial y\partial x} 
404 


A\hackscore{11}\frac{\partial^{2}u}{\partial y^{2}} 
405 


+ Du = f 
406 


\end{equation} 
407 
ahallam 
2606 
Notice that for the higher dimensional case $A$ becomes a matrix. It is also 
408 
ahallam 
2495 
important to notice that the usage of the Nabla operator creates 
409 


a compact formulation which is also independent from the spatial dimension. 
410 
artak 
2963 
To make the general PDE \refEq{eqn:commonform2D} one dimensional as 
411 
gross 
2861 
shown in \refEq{eqn:commonform} we need to set 
412 
ahallam 
2606 
\begin{equation} 
413 
ahallam 
2494 
A\hackscore{00}=A; A\hackscore{01}=A\hackscore{10}=A\hackscore{11}=0 
414 
gross 
2477 
\end{equation} 
415 



416 
gross 
2867 

417 
ahallam 
2495 
\subsection{Developing a PDE Solution Script} 
418 
ahallam 
2801 
\label{sec:key} 
419 
gross 
2949 
\sslist{example01a.py} 
420 
ahallam 
2979 
We write a simple \pyt script which uses the \modescript, \modfinley and \modmpl 
421 


modules. 
422 


By developing a script for \esc, the heat diffusion equation can be solved at 
423 


successive time steps for a predefined period using our general form 
424 


\refEq{eqn:hdgenf}. Firstly it is necessary to import all the 
425 


libraries\footnote{The libraries contain predefined scripts that are required to 
426 


solve certain problems, these can be simple like sine and cosine functions or 
427 


more complicated like those from our \esc library.} 
428 
ahallam 
2495 
that we will require. 
429 
ahallam 
2775 
\begin{python} 
430 
ahallam 
2495 
from esys.escript import * 
431 
ahallam 
2606 
# This defines the LinearPDE module as LinearPDE 
432 


from esys.escript.linearPDEs import LinearPDE 
433 


# This imports the rectangle domain function from finley. 
434 


from esys.finley import Rectangle 
435 


# A useful unit handling package which will make sure all our units 
436 


# match up in the equations under SI. 
437 


from esys.escript.unitsSI import * 
438 
ahallam 
2775 
\end{python} 
439 
ahallam 
2979 
It is generally a good idea to import all of the \modescript library, although 
440 


if the functions and classes required are known they can be specified 
441 


individually. The function \verbLinearPDE has been imported explicitly for 
442 


ease of use later in the script. \verbRectangle is going to be our type of 
443 


domain. The module \verb unitsSI provides support for SI unit definitions with 
444 


our variables. 
445 
gross 
2477 

446 
ahallam 
2979 
Once our library dependencies have been established, defining the problem 
447 


specific variables is the next step. In general the number of variables needed 
448 


will vary between problems. These variables belong to two categories. They are 
449 


either directly related to the PDE and can be used as inputs into the \esc 
450 


solver, or they are script variables used to control internal functions and 
451 


iterations in our problem. For this PDE there are a number of constants which 
452 


need values. Firstly, the domain upon which we wish to solve our problem needs 
453 


to be defined. There are different types of domains in \modescript which we 
454 


demonstrate in later tutorials but for our granite blocks, we simply use a 
455 


rectangular domain. 
456 
ahallam 
2401 

457 
ahallam 
2979 
Using a rectangular domain simplifies our granite blocks (which would in reality 
458 


be a \textit{3D} object) into a single dimension. The granite blocks will have a 
459 


lengthways cross section that looks like a rectangle. As a result we do not 
460 


need to model the volume of the block due to symmetry. There are four arguments 
461 


we must consider when we decide to create a rectangular domain, the domain 
462 


\textit{length}, \textit{width} and \textit{step size} in each direction. When 
463 


defining the size of our problem it will help us determine appropriate values 
464 


for our model arguments. If we make our dimensions large but our step sizes very 
465 


small we increase the accuracy of our solution. Unfortunately we also increase 
466 


the number of calculations that must be solved per time step. This means more 
467 


computational time is required to produce a solution. In this \textit{1D} 
468 


problem, the bar is defined as being 1 metre long. An appropriate step size 
469 


\verbndx would be 1 to 10\% of the length. Our \verbndy need only be 1, this 
470 


is because our problem stipulates no partial derivatives in the $y$ direction. 
471 


Thus the temperature does not vary with $y$. Hence, the model parameters can be 
472 


defined as follows; note we have used the \verb unitsSI convention to make sure 
473 


all our input units are converted to SI. 
474 
ahallam 
2775 
\begin{python} 
475 
gross 
2905 
mx = 500.*m #meters  model length 
476 


my = 100.*m #meters  model width 
477 


ndx = 50 # mesh steps in x direction 
478 


ndy = 1 # mesh steps in y direction 
479 


boundloc = mx/2 # location of boundary between the two blocks 
480 
ahallam 
2775 
\end{python} 
481 
ahallam 
2979 
The material constants and the temperature variables must also be defined. For 
482 


the granite in the model they are defined as: 
483 
ahallam 
2775 
\begin{python} 
484 
ahallam 
2495 
#PDE related 
485 
gross 
2905 
rho = 2750. *kg/m**3 #kg/m^{3} density of iron 
486 


cp = 790.*J/(kg*K) # J/Kg.K thermal capacity 
487 
gross 
2878 
rhocp = rho*cp 
488 
gross 
2905 
kappa = 2.2*W/m/K # watts/m.Kthermal conductivity 
489 
gross 
2878 
qH=0 * J/(sec*m**3) # J/(sec.m^{3}) no heat source 
490 
gross 
2905 
T1=20 * Celsius # initial temperature at Block 1 
491 


T2=2273. * Celsius # base temperature at Block 2 
492 
ahallam 
2775 
\end{python} 
493 
ahallam 
2979 
Finally, to control our script we will have to specify our timing controls and 
494 


where we would like to save the output from the solver. This is simple enough: 
495 
ahallam 
2775 
\begin{python} 
496 
gross 
2878 
t=0 * day #our start time, usually zero 
497 


tend=1. * day #  time to end simulation 
498 
ahallam 
2495 
outputs = 200 # number of time steps required. 
499 


h=(tendt)/outputs #size of time step 
500 
ahallam 
2606 
#user warning statement 
501 


print "Expected Number of time outputs is: ", (tendt)/h 
502 


i=0 #loop counter 
503 
ahallam 
2775 
\end{python} 
504 
ahallam 
2979 
Now that we know our inputs we will build a domain using the 
505 


\verb Rectangle() function from \verb finley . The four arguments allow us to define our domain 
506 


\verb model as: 
507 
ahallam 
2775 
\begin{python} 
508 
ahallam 
2606 
#generate domain using rectangle 
509 
gross 
2905 
blocks = Rectangle(l0=mx,l1=my,n0=ndx, n1=ndy) 
510 
ahallam 
2775 
\end{python} 
511 
artak 
2963 
\verb blocks now describes a domain in the manner of Section \ref{ss:domcon}. 
512 
ahallam 
2658 

513 
ahallam 
2979 
With a domain and all our required variables established, it is now possible to 
514 


set up our PDE so that it can be solved by \esc. The first step is to define the 
515 


type of PDE that we are trying to solve in each time step. In this example it is 
516 


a single linear PDE\footnote{in comparison to a system of PDEs which we discuss 
517 


later.}. We also need to state the values of our general form variables. 
518 
ahallam 
2775 
\begin{python} 
519 
gross 
2905 
mypde=LinearPDE(blocks) 
520 
gross 
2878 
A=zeros((2,2))) 
521 


A[0,0]=kappa 
522 
gross 
2905 
mypde.setValue(A=A, D=rhocp/h) 
523 
ahallam 
2775 
\end{python} 
524 
ahallam 
2979 
In a many cases it may be possible to decrease the computational time of the 
525 


solver if the PDE is symmetric. 
526 
gross 
2878 
Symmetry of a PDE is defined by; 
527 
ahallam 
2495 
\begin{equation}\label{eqn:symm} 
528 


A\hackscore{jl}=A\hackscore{lj} 
529 


\end{equation} 
530 
ahallam 
2979 
Symmetry is only dependent on the $A$ coefficient in the general form and the 
531 


other coefficients $D$ as well as the right hand side $Y$. From the above 
532 


definition we can see that our PDE is symmetric. The \verb LinearPDE class 
533 


provides the method \method{checkSymmetry} to check if the given PDE is 
534 


symmetric. As our PDE is symmetrical we enable symmetry via; 
535 
ahallam 
2775 
\begin{python} 
536 
ahallam 
2495 
myPDE.setSymmetryOn() 
537 
ahallam 
2775 
\end{python} 
538 
ahallam 
2979 
Next we need to establish the initial temperature distribution \verbT. We need 
539 


to 
540 


assign the value \verbT1 to all sample points left to the contact interface at 
541 


$x\hackscore{0}=\frac{mx}{2}$ 
542 
gross 
2905 
and the value \verbT2 right to the contact interface. \esc 
543 


provides the \verbwhereNegative function to construct this. In fact, 
544 
ahallam 
2979 
\verbwhereNegative returns the value $1$ at those sample points where the 
545 


argument 
546 


has a negative value. Otherwise zero is returned. If \verbx are the 
547 


$x\hackscore{0}$ 
548 
gross 
2905 
coordinates of the sample points used to represent the temperature distribution 
549 


then \verbx[0]boundloc gives us a negative value for 
550 


all sample points left to the interface and nonnegative value to 
551 


the right of the interface. So with; 
552 
gross 
2878 
\begin{python} 
553 


# ... set initial temperature .... 
554 
gross 
2905 
T= T1*whereNegative(x[0]boundloc)+T2*(1whereNegative(x[0]boundloc)) 
555 
gross 
2878 
\end{python} 
556 
ahallam 
2979 
we get the desired temperature distribution. To get the actual sample points 
557 


\verbx we use 
558 
gross 
2905 
the \verbgetX() method of the function space \verbSolution(blocks) 
559 


which is used to represent the solution of a PDE; 
560 
gross 
2878 
\begin{python} 
561 
gross 
2905 
x=Solution(blocks).getX() 
562 


\end{python} 
563 
ahallam 
2979 
As \verbx are the sample points for the function space 
564 


\verbSolution(blocks) 
565 


the initial temperature \verbT is using these sample points for 
566 


representation. 
567 


Although \esc is trying to be forgiving with the choice of sample points and to 
568 


convert 
569 


where necessary the adjustment of the function space is not always possible. So 
570 


it is 
571 
gross 
2905 
advisable to make a careful choice on the function space used. 
572 



573 
ahallam 
2979 
Finally we initialise an iteration loop to solve our PDE for all the time steps 
574 


we specified in the variable section. As the right hand side of the general form 
575 


is dependent on the previous values for temperature \verb T across the bar this 
576 


must be updated in the loop. Our output at each time step is \verb T the heat 
577 


distribution and \verb totT the total heat in the system. 
578 
gross 
2905 
\begin{python} 
579 
gross 
2878 
while t < tend: 
580 


i+=1 #increment the counter 
581 


t+=h #increment the current time 
582 


mypde.setValue(Y=qH+rhocp/h*T) # set variable PDE coefficients 
583 


T=mypde.getSolution() #get the PDE solution 
584 


totE = integrate(rhocp*T) #get the total heat (energy) in the system 
585 


\end{python} 
586 
ahallam 
2979 
The last statement in this script calculates the total energy in the system as 
587 


volume integral 
588 


of $\rho c\hackscore{p} T$ over the block. As the blocks are insulated no energy 
589 


should be get lost or added. 
590 
gross 
2905 
The total energy should stay constant for the example discussed here. 
591 
ahallam 
2401 

592 
gross 
2905 
\subsection{Running the Script} 
593 
ahallam 
2975 
The script presented so far is available under 
594 
gross 
2949 
\verbexample01a.py. You can edit this file with your favourite text editor. 
595 
ahallam 
2979 
On most operating systems\footnote{The you can use \texttt{runescript} launcher 
596 


is not supported under {\it MS Windows} yet.} you can use the 
597 


\program{runescript} command 
598 
gross 
2905 
to launch {\it escript} scripts. For the example script use; 
599 


\begin{verbatim} 
600 
gross 
2949 
runescript example01a.py 
601 
gross 
2905 
\end{verbatim} 
602 


The program will print a progress report. Alternatively, you can use 
603 


the python interpreter directly; 
604 


\begin{verbatim} 
605 
gross 
2949 
python example01a.py 
606 
gross 
2905 
\end{verbatim} 
607 
ahallam 
2979 
if the system is configured correctly (Please talk to your system 
608 


administrator). 
609 
gross 
2905 

610 


\begin{figure} 
611 


\begin{center} 
612 


\includegraphics[width=4in]{figures/ttblockspyplot150} 
613 
gross 
2952 
\caption{Example 1b: Total Energy in the Blocks over Time (in seconds).} 
614 
gross 
2905 
\label{fig:onedheatout1} 
615 


\end{center} 
616 


\end{figure} 
617 



618 
gross 
2878 
\subsection{Plotting the Total Energy} 
619 
gross 
2949 
\sslist{example01b.py} 
620 
gross 
2878 

621 
ahallam 
2979 
\esc does not include its own plotting capabilities. However, it is possible to 
622 


use a variety of free \pyt packages for visualisation. 
623 


Two types will be demonstrated in this cookbook; 
624 


\mpl\footnote{\url{http://matplotlib.sourceforge.net/}} and 
625 


\verb VTK \footnote{\url{http://www.vtk.org/}} visualisation. 
626 


The \mpl package is a component of SciPy\footnote{\url{http://www.scipy.org}} 
627 


and is good for basic graphs and plots. 
628 


For more complex visualisation tasks in particular, two and three dimensional 
629 


problems we recommend the use of more advanced tools. For instance, \mayavi 
630 


\footnote{\url{http://code.enthought.com/projects/mayavi/}} 
631 
ahallam 
2975 
which is based upon the \verbVTK toolkit. The usage of \verbVTK based 
632 
ahallam 
2979 
visualization is discussed in Chapter~\ref{Sec:2DHD} which focusses on a two 
633 


dimensional PDE. 
634 
gross 
2878 

635 
ahallam 
2979 
For our simple granite block problem, we have two plotting tasks. Firstly, we 
636 


are interested in showing the 
637 


behaviour of the total energy over time and secondly, how the temperature 
638 


distribution within the block is 
639 
gross 
2878 
developing over time. Lets start with the first task. 
640 



641 
ahallam 
2979 
The trick is to create a record of the time marks and the corresponding total 
642 


energies observed. 
643 
gross 
2878 
\pyt provides the concept of lists for this. Before 
644 
ahallam 
2979 
the time loop is opened we create empty lists for the time marks \verbt_list 
645 


and the total energies \verbE_list. 
646 


After the new temperature has been calculated by solving the PDE we append the 
647 


new time marker and the total energy value for that time 
648 


to the corresponding list using the \verbappend method. With these 
649 


modifications our script looks as follows: 
650 
ahallam 
2775 
\begin{python} 
651 
gross 
2878 
t_list=[] 
652 


E_list=[] 
653 


# ... start iteration: 
654 


while t<tend: 
655 


t+=h 
656 


mypde.setValue(Y=qH+rhocp/h*T) # set variable PDE coefficients 
657 


T=mypde.getSolution() #get the PDE solution 
658 


totE=integrate(rhocp*T) 
659 


t_list.append(t) # add current time mark to record 
660 


E_list.append(totE) # add current total energy to record 
661 
ahallam 
2775 
\end{python} 
662 
ahallam 
2979 
To plot $t$ over $totE$ we use \mpl a module contained within \pylab which needs 
663 


to be loaded before used; 
664 
ahallam 
2775 
\begin{python} 
665 
gross 
2878 
import pylab as pl # plotting package. 
666 
ahallam 
2775 
\end{python} 
667 
ahallam 
2979 
Here we are not using the \verbfrom pylab import * in order to avoid name 
668 


clashes for function names 
669 
gross 
2905 
within \esc. 
670 
ahallam 
2401 

671 
ahallam 
2979 
The following statements are added to the script after the time loop has been 
672 


completed; 
673 
ahallam 
2775 
\begin{python} 
674 
gross 
2878 
pl.plot(t_list,E_list) 
675 


pl.title("Total Energy") 
676 
gross 
2905 
pl.axis([0,max(t_list),0,max(E_list)*1.1]) 
677 
gross 
2878 
pl.savefig("totE.png") 
678 
ahallam 
2775 
\end{python} 
679 
ahallam 
2979 
The first statement hands over the time marks and corresponding total energies 
680 


to the plotter. 
681 
gross 
2905 
The second statment is setting the title for the plot. The third statement 
682 
ahallam 
2979 
sets the axis ranges. In most cases these are set appropriately by the plotter. 
683 



684 
gross 
2905 
The last statement renders the plot and writes the 
685 
ahallam 
2979 
result into the file \verbtotE.png which can be displayed by (almost) any 
686 


image viewer. 
687 


As expected the total energy is constant over time, see 
688 


\reffig{fig:onedheatout1}. 
689 
ahallam 
2401 

690 
gross 
2878 
\subsection{Plotting the Temperature Distribution} 
691 
gross 
2931 
\label{sec: plot T} 
692 
gross 
2949 
\sslist{example01c.py} 
693 
ahallam 
2979 
For plotting the spatial distribution of the temperature we need to modify the 
694 


strategy we have used 
695 


for the total energy. Instead of producing a final plot at the end we will 
696 


generate a 
697 


picture at each time step which can be browsed as a slide show or composed into 
698 


a movie. 
699 


The first problem we encounter is that if we produce an image at each time step 
700 


we need 
701 
gross 
2878 
to make sure that the images previously generated are not overwritten. 
702 



703 
ahallam 
2979 
To develop an incrementing file name we can use the following convention. It is 
704 


convenient to 
705 


put all image file showing the same variable  in our case the temperature 
706 


distribution  
707 


into a separate directory. As part of the \verbos module\footnote{The 
708 


\texttt{os} module provides 
709 


a powerful interface to interact with the operating system, see 
710 


\url{http://docs.python.org/library/os.html}.} \pyt 
711 
gross 
2878 
provides the \verbos.path.join command to build file and 
712 


directory names in a platform independent way. Assuming that 
713 
ahallam 
2979 
\verbsave_path is name of directory we want to put the results the command 
714 


is; 
715 
ahallam 
2775 
\begin{python} 
716 
gross 
2878 
import os 
717 


os.path.join(save_path, "tempT%03d.png"%i ) 
718 
ahallam 
2775 
\end{python} 
719 
gross 
2878 
where \verbi is the time step counter. 
720 
ahallam 
2979 
There are two arguments to the \verb join command. The 
721 


\verb save_path variable is a predefined string pointing to the directory we want to save our 
722 


data, for example a single subfolder called \verb data would be defined by; 
723 
gross 
2878 
\begin{verbatim} 
724 


save_path = "data" 
725 


\end{verbatim} 
726 
gross 
2949 
while a subfolder of \verb data called \verb example01 would be defined by; 
727 
gross 
2878 
\begin{verbatim} 
728 
gross 
2949 
save_path = os.path.join("data","example01") 
729 
gross 
2878 
\end{verbatim} 
730 
ahallam 
2979 
The second argument of \verb join \xspace contains a string which is the file 
731 


name or subdirectory name. We can use the operator \verb% to use the value of 
732 


\verbi as part of our filename. The substring \verb %03d indicates that we 
733 


want to substitude a value into the name; 
734 
gross 
2878 
\begin{itemize} 
735 
jfenwick 
2961 
\item \verb 0 means that small numbers should have leading zeroes; 
736 
ahallam 
2979 
\item \verb 3 means that numbers should be written using at least 3 digits; 
737 


and 
738 
jfenwick 
2961 
\item \verb d means that the value to substitute will be an integer. 
739 
gross 
2878 
\end{itemize} 
740 
ahallam 
2979 
To actually substitute the value of \verbi into the name write \verb %i after 
741 


the string. 
742 


When done correctly, the output files from this command would be place in the 
743 


directory defined by \verb save_path as; 
744 
gross 
2878 
\begin{verbatim} 
745 
jfenwick 
2961 
blockspyplot001.png 
746 


blockspyplot002.png 
747 


blockspyplot003.png 
748 
gross 
2878 
... 
749 


\end{verbatim} 
750 


and so on. 
751 



752 


A subfolder check/constructor is available in \esc. The command; 
753 


\begin{verbatim} 
754 


mkDir(save_path) 
755 


\end{verbatim} 
756 
ahallam 
2979 
will check for the existence of \verb save_path and if missing, make the 
757 


required directories. 
758 
gross 
2878 

759 
artak 
2964 
We start by modifying our solution script. 
760 
ahallam 
2979 
Prior to the \verbwhile loop we will need to extract our finite solution 
761 


points to a data object that is compatible with \mpl. First we create the node 
762 


coordinates of the sample points used to represent 
763 


the temperature as a \pyt list of tuples or a \numpy array as requested by the 
764 


plotting function. 
765 


We need to convert the array \verbx previously set as 
766 


\verbSolution(blocks).getX() into a \pyt list 
767 


and then to a \numpy array. The $x\hackscore{0}$ component is then extracted via 
768 


an array slice to the variable \verbplx; 
769 
ahallam 
2775 
\begin{python} 
770 
gross 
2878 
import numpy as np # array package. 
771 


#convert solution points for plotting 
772 
gross 
2905 
plx = x.toListOfTuples() 
773 
gross 
2878 
plx = np.array(plx) # convert to tuple to numpy array 
774 


plx = plx[:,0] # extract x locations 
775 
ahallam 
2775 
\end{python} 
776 
gross 
2878 

777 
ahallam 
2645 
\begin{figure} 
778 


\begin{center} 
779 
gross 
2905 
\includegraphics[width=4in]{figures/blockspyplot001} 
780 


\includegraphics[width=4in]{figures/blockspyplot050} 
781 


\includegraphics[width=4in]{figures/blockspyplot200} 
782 
ahallam 
2979 
\caption{Example 1c: Temperature ($T$) distribution in the blocks at time steps 
783 


$1$, $50$ and $200$.} 
784 
ahallam 
2645 
\label{fig:onedheatout} 
785 


\end{center} 
786 


\end{figure} 
787 



788 
ahallam 
2979 
We use the same techniques provided by \mpl as we have used to plot the total 
789 


energy over time. 
790 


For each time step we generate a plot of the temperature distribution and save 
791 


each to a file. 
792 


The following is appended to the end of the \verb while loop and creates one 
793 


figure of the temperature distribution. We start by converting the solution to a 
794 


tuple and then plotting this against our \textit{x coordinates} \verb plx we 
795 


have generated before. We add a title to the diagram before it is rendered into 
796 


a file. 
797 


Finally, the figure is saved to a \verb*.png file and cleared for the 
798 


following iteration. 
799 
ahallam 
2775 
\begin{python} 
800 
gross 
2878 
# ... start iteration: 
801 


while t<tend: 
802 


.... 
803 


T=mypde.getSolution() #get the PDE solution 
804 


tempT = T.toListOfTuples() # convert to a tuple 
805 


pl.plot(plx,tempT) # plot solution 
806 


# set scale (Temperature should be between Tref and T0) 
807 


pl.axis([0,mx,Tref*.9,T0*1.1]) 
808 


# add title 
809 
gross 
2905 
pl.title("Temperature across the blocks at time %e minutes"%(t/day)) 
810 
gross 
2878 
#save figure to file 
811 
gross 
2905 
pl.savefig(os.path.join(save_path,"tempT","blockspyplot%03d.png") %i) 
812 
ahallam 
2775 
\end{python} 
813 
gross 
2905 
Some results are shown in \reffig{fig:onedheatout}. 
814 
jfenwick 
2657 

815 
ahallam 
2606 
\subsection{Make a video} 
816 
ahallam 
2979 
Our saved plots from the previous section can be cast into a video using the 
817 


following command appended to the end of the script. The \verb mencoder command 
818 


is Linux only, so other platform users need to use an alternative video encoder. 
819 
ahallam 
2775 
\begin{python} 
820 
gross 
2878 
# compile the *.png files to create a *.avi videos that show T change 
821 


# with time. This operation uses Linux mencoder. For other operating 
822 
gross 
2905 
# systems it is possible to use your favourite video compiler to 
823 
ahallam 
2606 
# convert image files to videos. 
824 
gross 
2477 

825 
ahallam 
2606 
os.system("mencoder mf://"+save_path+"/tempT"+"/*.png mf type=png:\ 
826 
gross 
2878 
w=800:h=600:fps=25 ovc lavc lavcopts vcodec=mpeg4 oac copy o \ 
827 
gross 
2949 
example01tempT.avi") 
828 
ahallam 
2775 
\end{python} 
829 
gross 
2477 
