1 

2 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
3 
% 
4 
% Copyright (c) 20032012 by University of Queensland 
5 
% 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 
\section{Example 1: One Dimensional Heat Diffusion in Granite} 
15 
\label{Sec:1DHDv00} 
16 

17 
The first model consists of two blocks of isotropic material, for instance 
18 
granite, sitting next to each other (\autoref{fig:onedgbmodel}). 
19 
Initial temperature in \textit{Block 1} is \verbT1 and in \textit{Block 2} is 
20 
\verbT2. 
21 
We assume that the system is insulated. 
22 
What would happen to the temperature distribution in each block over time? 
23 
Intuition tells us that heat will be transported from the hotter block to the 
24 
cooler one until both 
25 
blocks have the same temperature. 
26 

27 
\begin{figure}[ht] 
28 
\centerline{\includegraphics[width=4.in]{figures/onedheatdiff001}} 
29 
\caption{Example 1: Temperature differential along a single interface between 
30 
two granite blocks.} 
31 
\label{fig:onedgbmodel} 
32 
\end{figure} 
33 

34 
\subsection{1D Heat Diffusion Equation} 
35 
We can model the heat distribution of this problem over time using the 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 
which is defined as: 
41 
\begin{equation} 
42 
\rho c_p \frac{\partial T}{\partial t}  \kappa \frac{\partial^{2} 
43 
T}{\partial x^{2}} = q_H 
44 
\label{eqn:hd} 
45 
\end{equation} 
46 
where $\rho$ is the material density, $c_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 
parameters are \textbf{constant}. 
53 
The heat source is defined by the right hand side of \refEq{eqn:hd} as 
54 
$q_{H}$; this can take the form of a constant or a function of time and 
55 
space. For example $q_{H} = q_{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 blockblock interface $x$. 
62 

63 
\subsection{PDEs and the General Form} 
64 
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 

74 
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 

78 
For time discretisation we use the Backward Euler approximation 
79 
scheme\footnote{see \url{http://en.wikipedia.org/wiki/Euler_method}}. It is 
80 
based on the approximation 
81 
\begin{equation} 
82 
\frac{\partial T(t)}{\partial t} \approx \frac{T(t)T(th)}{h} 
83 
\label{eqn:beuler} 
84 
\end{equation} 
85 
for $\frac{\partial T}{\partial t}$ at time $t$ 
86 
where $h$ is the time step size. This can also be written as; 
87 
\begin{equation} 
88 
\frac{\partial T}{\partial t}(t^{(n)}) \approx \frac{T^{(n)}  T^{(n1)}}{h} 
89 
\label{eqn:Tbeuler} 
90 
\end{equation} 
91 
where the upper index $n$ denotes the n\textsuperscript{th} time step. So one 
92 
has 
93 
\begin{equation} 
94 
\begin{array}{rcl} 
95 
t^{(n)} & = & t^{(n1)}+h \\ 
96 
T^{(n)} & = & T(t^{(n1)}) \\ 
97 
\end{array} 
98 
\label{eqn:Neuler} 
99 
\end{equation} 
100 
Substituting \refEq{eqn:Tbeuler} into \refEq{eqn:hd} we get; 
101 
\begin{equation} 
102 
\frac{\rho c_p}{h} (T^{(n)}  T^{(n1)})  \kappa \frac{\partial^{2} 
103 
T^{(n)}}{\partial x^{2}} = q_H 
104 
\label{eqn:hddisc} 
105 
\end{equation} 
106 
Notice that we evaluate the spatial derivative term at the current time 
107 
$t^{(n)}$  therefore the name \textbf{backward Euler} scheme. Alternatively, 
108 
one can evaluate the spatial derivative term at the previous time $t^{(n1)}$. 
109 
This approach is called the \textbf{forward Euler} scheme. This scheme can 
110 
provide some computational advantages, which 
111 
are not discussed here. However, the \textbf{forward Euler} scheme has a major 
112 
disadvantage. Namely, depending on the 
113 
material parameters as well as the domain discretization of the spatial 
114 
derivative term, the time step size $h$ needs to be chosen sufficiently small to 
115 
achieve a stable temperature when progressing in time. Stability is achieved if 
116 
the temperature does not grow beyond its initial bounds and becomes 
117 
nonphysical. 
118 
The backward Euler scheme, which we use here, is unconditionally stable meaning 
119 
that under the assumption of a 
120 
physically correct problem setup the temperature approximation remains physical 
121 
for all time steps. 
122 
The user needs to keep in mind that the discretisation error introduced by 
123 
\refEq{eqn:beuler} 
124 
is sufficiently small, thus a good approximation of the true temperature is 
125 
computed. It is therefore very important that any results are viewed with 
126 
caution. For example, one may compare the results for different time and 
127 
spatial step sizes. 
128 

129 
To get the temperature $T^{(n)}$ at time $t^{(n)}$ we need to solve the linear 
130 
differential equation \refEq{eqn:hddisc} which only includes spatial 
131 
derivatives. To solve this problem we want to use \esc. 
132 

133 
In \esc any given PDE can be described by the general form. For the purpose of 
134 
this introduction we illustrate a simpler version of the general form for full 
135 
linear PDEs which is available in the \esc user's guide. A simplified form that 
136 
suits our heat diffusion problem\footnote{The form in the \esc users guide which 
137 
uses the Einstein convention is written as 
138 
$(A_{jl} u_{,l})_{,j}+D u =Y$} 
139 
is described by; 
140 
\begin{equation}\label{eqn:commonform nabla} 
141 
\nabla\cdot(A\cdot\nabla u) + Du = f 
142 
\end{equation} 
143 
where $A$, $D$ and $f$ are known values and $u$ is the unknown solution. The 
144 
symbol $\nabla$ which is called the \textit{Nabla operator} or \textit{del 
145 
operator} represents 
146 
the spatial derivative of its subject  in this case $u$. Lets assume for a 
147 
moment that we deal with a onedimensional problem then ; 
148 
\begin{equation} 
149 
\nabla = \frac{\partial}{\partial x} 
150 
\end{equation} 
151 
and we can write \refEq{eqn:commonform nabla} as; 
152 
\begin{equation}\label{eqn:commonform} 
153 
A\frac{\partial^{2}u}{\partial x^{2}} + Du = f 
154 
\end{equation} 
155 
if $A$ is constant. To match this simplified general form to our problem 
156 
\refEq{eqn:hddisc} 
157 
we rearrange \refEq{eqn:hddisc}; 
158 
\begin{equation} 
159 
\frac{\rho c_p}{h} T^{(n)}  \kappa \frac{\partial^2 T^{(n)}}{\partial 
160 
x^2} = q_H + \frac{\rho c_p}{h} T^{(n1)} 
161 
\label{eqn:hdgenf} 
162 
\end{equation} 
163 
The PDE is now in a form that satisfies \refEq{eqn:commonform nabla} which is 
164 
required for \esc to solve our PDE. This can be done by generating a solution 
165 
for successive increments in the time nodes $t^{(n)}$ where 
166 
$t^{(0)}=0$ and $t^{(n)}=t^{(n1)}+h$ where $h>0$ is the step size and assumed 
167 
to be constant. 
168 
In the following the upper index ${(n)}$ refers to a value at time $t^{(n)}$. 
169 
Finally, by comparing \refEq{eqn:hdgenf} with \refEq{eqn:commonform} one can see 
170 
that; 
171 
\begin{equation}\label{ESCRIPT SET} 
172 
u=T^{(n)}; 
173 
A = \kappa; D = \frac{\rho c _{p}}{h}; f = q _{H} + \frac{\rho 
174 
c_p}{h} T^{(n1)} 
175 
\end{equation} 
176 

177 
\subsection{Boundary Conditions} 
178 
\label{SEC BOUNDARY COND} 
179 
With the PDE sufficiently modified, consideration must now be given to the 
180 
boundary conditions of our model. Typically there are two main types of boundary 
181 
conditions known as \textbf{Neumann} and \textbf{Dirichlet} boundary 
182 
conditions\footnote{More information on Boundary Conditions is available at 
183 
Wikipedia \url{http://en.wikipedia.org/wiki/Boundary_conditions}}, 
184 
respectively. 
185 
A \textbf{Dirichlet boundary condition} is conceptually simpler and is used to 
186 
prescribe a known value to the unknown solution (in our example the temperature) 
187 
on parts of the boundary or on the entire boundary of the region of interest. 
188 
We discuss the Dirichlet boundary condition in our second example presented in 
189 
Section~\ref{Sec:1DHDv0}. 
190 

191 
However, for this example we have made the model assumption that the system is 
192 
insulated, so we need to add an appropriate boundary condition to prevent 
193 
any loss or inflow of energy at the boundary of our domain. Mathematically this 
194 
is expressed by prescribing 
195 
the heat flux $\kappa \frac{\partial T}{\partial x}$ to zero. In our simplified 
196 
one dimensional model this is expressed 
197 
in the form; 
198 
\begin{equation} 
199 
\kappa \frac{\partial T}{\partial x} = 0 
200 
\end{equation} 
201 
or in a more general case as 
202 
\begin{equation}\label{NEUMAN 1} 
203 
\kappa \nabla T \cdot n = 0 
204 
\end{equation} 
205 
where $n$ is the outer normal field \index{outer normal field} at the surface 
206 
of the domain. 
207 
The $\cdot$ (dot) refers to the dot product of the vectors $\nabla T$ and $n$. 
208 
In fact, the term $\nabla T \cdot n$ is the normal derivative of 
209 
the temperature $T$. Other notations used here are\footnote{The \esc notation 
210 
for the normal 
211 
derivative is $T_{,i} n_i$.}; 
212 
\begin{equation} 
213 
\nabla T \cdot n = \frac{\partial T}{\partial n} \; . 
214 
\end{equation} 
215 
A condition of the type \refEq{NEUMAN 1} defines a \textbf{Neumann boundary 
216 
condition} for the PDE. 
217 

218 
The PDE \refEq{eqn:hdgenf} 
219 
and the Neumann boundary condition~\ref{eqn:hdgenf} (potentially together with 
220 
the Dirichlet boundary conditions) define a \textbf{boundary value problem}. 
221 
It is the nature of a boundary value problem to allow making statements about 
222 
the solution in the 
223 
interior of the domain from information known on the boundary only. In most 
224 
cases we use the term partial differential equation but in fact it is a 
225 
boundary value problem. 
226 
It is important to keep in mind that boundary conditions need to be complete and 
227 
consistent in the sense that 
228 
at any point on the boundary either a Dirichlet or a Neumann boundary condition 
229 
must be set. 
230 

231 
Conveniently, \esc makes a default assumption on the boundary conditions which 
232 
the user may modify where appropriate. 
233 
For a problem of the form in~\refEq{eqn:commonform nabla} the default 
234 
condition\footnote{In the \esc user guide which uses the Einstein convention 
235 
this is written as 
236 
$n_{j}A_{jl} u_{,l}=0$.} is; 
237 
\begin{equation}\label{NEUMAN 2} 
238 
n\cdot A \cdot\nabla u = 0 
239 
\end{equation} 
240 
which is used everywhere on the boundary. Again $n$ denotes the outer normal 
241 
field. 
242 
Notice that the coefficient $A$ is the same as in the \esc 
243 
PDE~\ref{eqn:commonform nabla}. 
244 
With the settings for the coefficients we have already identified in 
245 
\refEq{ESCRIPT SET} this 
246 
condition translates into 
247 
\begin{equation}\label{NEUMAN 2b} 
248 
\kappa \frac{\partial T}{\partial x} = 0 
249 
\end{equation} 
250 
for the boundary of the domain. This is identical to the Neumann boundary 
251 
condition we want to set. \esc will take care of this condition for us. We 
252 
discuss the Dirichlet boundary condition later. 
253 

254 
\subsection{Outline of the Implementation} 
255 
\label{sec:outline} 
256 
To solve the heat diffusion equation (\refEq{eqn:hd}) we write a simple \pyt 
257 
script. At this point we assume that you have some basic understanding of the 
258 
\pyt programming language. If not, there are some pointers and links available 
259 
in Section \ref{sec:escpybas}. The script (discussed in \refSec{sec:key}) has 
260 
four major steps. Firstly, we need to define the domain where we want to 
261 
calculate the temperature. For our problem this is the joint blocks of granite 
262 
which has a rectangular shape. Secondly, we need to define the PDE to solve in 
263 
each time step to get the updated temperature. Thirdly, we need to define the 
264 
coefficients of the PDE and finally we need to solve the PDE. The last two steps 
265 
need to be repeated until the final time marker has been reached. The work flow 
266 
is described in \reffig{fig:wf}. 
267 
% \begin{enumerate} 
268 
% \item create domain 
269 
% \item create PDE 
270 
% \item while end time not reached: 
271 
% \begin{enumerate} 
272 
% \item set PDE coefficients 
273 
% \item solve PDE 
274 
% \item update time marker 
275 
% \end{enumerate} 
276 
% \item end of calculation 
277 
% \end{enumerate} 
278 

279 
\begin{figure}[h!] 
280 
\centering 
281 
\includegraphics[width=1in]{figures/workflow.png} 
282 
\caption{Workflow for developing an \esc model and solution} 
283 
\label{fig:wf} 
284 
\end{figure} 
285 

286 
In the terminology of \pyt, the domain and PDE are represented by 
287 
\textbf{objects}. The nice feature of an object is that it is defined by its 
288 
usage and features 
289 
rather than its actual representation. So we will create a domain object to 
290 
describe the geometry of the two 
291 
granite blocks. Then we define PDEs and spatially distributed values such as the 
292 
temperature 
293 
on this domain. Similarly, to define a PDE object we use the fact that one needs 
294 
only to define the coefficients of the PDE and solve the PDE. The PDE object has 
295 
advanced features, but these are not required in simple cases. 
296 

297 

298 
\begin{figure}[htp] 
299 
\centering 
300 
\includegraphics[width=6in]{figures/functionspace.pdf} 
301 
\caption{\esc domain construction overview} 
302 
\label{fig:fs} 
303 
\end{figure} 
304 

305 
\subsection{The Domain Constructor in \esc} 
306 
\label{ss:domcon} 
307 
Whilst it is not strictly relevant or necessary, a better understanding of 
308 
how values are spatially distributed (\textit{e.g.} Temperature) and how PDE 
309 
coefficients are interpreted in \esc can be helpful. 
310 

311 
There are various ways to construct domain objects. The simplest form is a 
312 
rectangular shaped region with a length and height. There is 
313 
a ready to use function for this named \verb rectangle(). Besides the spatial 
314 
dimensions this function requires to specify the number of 
315 
elements or cells to be used along the length and height, see \reffig{fig:fs}. 
316 
Any spatially distributed value 
317 
and the PDE is represented in discrete form using this element 
318 
representation\footnote{We use the finite element method (FEM), see 
319 
\url{http://en.wikipedia.org/wiki/Finite_element_method} for details.}. 
320 
Therefore we will have access to an approximation of the true PDE solution 
321 
only. 
322 
The quality of the approximation depends  besides other factors  mainly on the 
323 
number of elements being used. In fact, the 
324 
approximation becomes better when more elements are used. However, computational 
325 
cost grows with the number of 
326 
elements being used. It is therefore important that you find the right balance 
327 
between the demand in accuracy and acceptable resource usage. 
328 

329 
In general, one can think about a domain object as a composition of nodes and 
330 
elements. 
331 
As shown in \reffig{fig:fs}, an element is defined by the nodes that are used to 
332 
describe its vertices. 
333 
To represent spatially distributed values the user can use 
334 
the values at the nodes, at the elements in the interior of the domain or at the 
335 
elements located on the surface of the domain. 
336 
The different approach used to represent values is called \textbf{function 
337 
space} and is attached to all objects 
338 
in \esc representing a spatially distributed value such as the solution of 
339 
a PDE. The three function spaces we use at the moment are; 
340 
\begin{enumerate} 
341 
\item the nodes, called by \verbContinuousFunction(domain) ; 
342 
\item the elements/cells, called by \verbFunction(domain) ; and 
343 
\item the boundary, called by \verbFunctionOnBoundary(domain). 
344 
\end{enumerate} 
345 
A function space object such as \verbContinuousFunction(domain) has the method 
346 
\verbgetX attached to it. This method returns the 
347 
location of the socalled \textbf{sample points} used to represent values of the 
348 
particular function space. So the 
349 
call \verbContinuousFunction(domain).getX() will return the coordinates of the 
350 
nodes used to describe the domain while 
351 
\verbFunction(domain).getX() returns the coordinates of numerical 
352 
integration points within elements, see \reffig{fig:fs}. 
353 

354 
This distinction between different representations of spatially distributed 
355 
values 
356 
is important in order to be able to vary the degrees of smoothness in a PDE 
357 
problem. 
358 
The coefficients of a PDE do not need to be continuous, thus this qualifies as a 
359 
\verbFunction() type. 
360 
On the other hand a temperature distribution must be continuous and needs to be 
361 
represented with a \verbContinuousFunction() function space. 
362 
An influx may only be defined at the boundary and is therefore a 
363 
\verbFunctionOnBoundary() object. 
364 
\esc allows certain transformations of the function spaces. A 
365 
\verbContinuousFunction() can be transformed into a 
366 
\verbFunctionOnBoundary() or \verbFunction(). On the other hand there is 
367 
not enough information in a \verbFunctionOnBoundary() to transform it to a 
368 
\verbContinuousFunction(). 
369 
These transformations, which are called \textbf{interpolation} are invoked 
370 
automatically by \esc if needed. 
371 

372 
Later in this introduction we discuss how 
373 
to define specific areas of geometry with different materials which are 
374 
represented by different material coefficients such as the 
375 
thermal conductivities $\kappa$. A very powerful technique to define these types 
376 
of PDE 
377 
coefficients is tagging. Blocks of materials and boundaries can be named and 
378 
values can be defined on subregions based on their names. 
379 
This is a method for simplifying PDE coefficient and flux definitions. It makes 
380 
scripting much easier and we will discuss this technique in 
381 
Section~\ref{STEADYSTATE HEAT REFRACTION}. 
382 

383 

384 
\subsection{A Clarification for the 1D Case} 
385 
\label{SEC: 1D CLARIFICATION} 
386 
It is necessary for clarification that we revisit our general PDE from 
387 
\refeq{eqn:commonform nabla} for a two dimensional domain. \esc is inherently 
388 
designed to solve problems that are multidimensional and so 
389 
\refEq{eqn:commonform nabla} needs to be read as a higher dimensional problem. 
390 
In the case of two spatial dimensions the \textit{Nabla operator} has in fact 
391 
two components $\nabla = (\frac{\partial}{\partial x}, \frac{\partial}{\partial 
392 
y})$. Assuming the coefficient $A$ is constant, the \refEq{eqn:commonform nabla} 
393 
takes the following form; 
394 
\begin{equation}\label{eqn:commonform2D} 
395 
A_{00}\frac{\partial^{2}u}{\partial x^{2}} 
396 
A_{01}\frac{\partial^{2}u}{\partial x\partial y} 
397 
A_{10}\frac{\partial^{2}u}{\partial y\partial x} 
398 
A_{11}\frac{\partial^{2}u}{\partial y^{2}} 
399 
+ Du = f 
400 
\end{equation} 
401 
Notice that for the higher dimensional case $A$ becomes a matrix. It is also 
402 
important to notice that the usage of the Nabla operator creates 
403 
a compact formulation which is also independent from the spatial dimension. 
404 
To make the general PDE \refEq{eqn:commonform2D} one dimensional as 
405 
shown in \refEq{eqn:commonform} we need to set 
406 
\begin{equation} 
407 
A_{00}=A; A_{01}=A_{10}=A_{11}=0 
408 
\end{equation} 
409 

410 

411 
\subsection{Developing a PDE Solution Script} 
412 
\label{sec:key} 
413 
\sslist{example01a.py} 
414 
We write a simple \pyt script which uses the \modescript, \modfinley and \modmpl 
415 
modules. 
416 
By developing a script for \esc, the heat diffusion equation can be solved at 
417 
successive time steps for a predefined period using our general form 
418 
\refEq{eqn:hdgenf}. Firstly it is necessary to import all the 
419 
libraries\footnote{The libraries contain predefined scripts that are required to 
420 
solve certain problems, these can be simple like sine and cosine functions or 
421 
more complicated like those from our \esc library.} 
422 
that we will require. 
423 
\begin{python} 
424 
from esys.escript import * 
425 
# This defines the LinearPDE module as LinearPDE 
426 
from esys.escript.linearPDEs import LinearPDE 
427 
# This imports the rectangle domain function from finley. 
428 
from esys.finley import Rectangle 
429 
# A useful unit handling package which will make sure all our units 
430 
# match up in the equations under SI. 
431 
from esys.escript.unitsSI import * 
432 
\end{python} 
433 
It is generally a good idea to import all of the \modescript library, although 
434 
if the functions and classes required are known they can be specified 
435 
individually. The function \verbLinearPDE has been imported explicitly for 
436 
ease of use later in the script. \verbRectangle is going to be our type of 
437 
domain. The module \verbunitsSI provides support for SI unit definitions with 
438 
our variables. 
439 

440 
Once our library dependencies have been established, defining the problem 
441 
specific variables is the next step. In general the number of variables needed 
442 
will vary between problems. These variables belong to two categories. They are 
443 
either directly related to the PDE and can be used as inputs into the \esc 
444 
solver, or they are script variables used to control internal functions and 
445 
iterations in our problem. For this PDE there are a number of constants which 
446 
need values. Firstly, the domain upon which we wish to solve our problem needs 
447 
to be defined. There are different types of domains in \modescript which we 
448 
demonstrate in later tutorials but for our granite blocks, we simply use a 
449 
rectangular domain. 
450 

451 
Using a rectangular domain simplifies our granite blocks (which would in reality 
452 
be a \textit{3D} object) into a single dimension. The granite blocks will have a 
453 
lengthways cross section that looks like a rectangle. As a result we do not 
454 
need to model the volume of the block due to symmetry. There are four arguments 
455 
we must consider when we decide to create a rectangular domain, the domain 
456 
\textit{length}, \textit{width} and \textit{step size} in each direction. When 
457 
defining the size of our problem it will help us determine appropriate values 
458 
for our model arguments. If we make our dimensions large but our step sizes very 
459 
small we increase the accuracy of our solution. Unfortunately we also increase 
460 
the number of calculations that must be solved per time step. This means more 
461 
computational time is required to produce a solution. In this \textit{1D} 
462 
problem, the bar is defined as being 1 metre long. An appropriate step size 
463 
\verbndx would be 1 to 10\% of the length. Our \verbndy needs only be 1, 
464 
this is because our problem stipulates no partial derivatives in the $y$ 
465 
direction. 
466 
Thus the temperature does not vary with $y$. Hence, the model parameters can be 
467 
defined as follows; note we have used the \verbunitsSI convention to make sure 
468 
all our input units are converted to SI. 
469 
\begin{python} 
470 
mx = 500.*m #meters  model length 
471 
my = 100.*m #meters  model width 
472 
ndx = 50 # mesh steps in x direction 
473 
ndy = 1 # mesh steps in y direction 
474 
boundloc = mx/2 # location of boundary between the two blocks 
475 
\end{python} 
476 
The material constants and the temperature variables must also be defined. For 
477 
the granite in the model they are defined as: 
478 
\begin{python} 
479 
#PDE related 
480 
rho = 2750. *kg/m**3 #kg/m^{3} density of iron 
481 
cp = 790.*J/(kg*K) # J/Kg.K thermal capacity 
482 
rhocp = rho*cp 
483 
kappa = 2.2*W/m/K # watts/m.Kthermal conductivity 
484 
qH=0 * J/(sec*m**3) # J/(sec.m^{3}) no heat source 
485 
T1=20 * Celsius # initial temperature at Block 1 
486 
T2=2273. * Celsius # base temperature at Block 2 
487 
\end{python} 
488 
Finally, to control our script we will have to specify our timing controls and 
489 
where we would like to save the output from the solver. This is simple enough: 
490 
\begin{python} 
491 
t=0 * day # our start time, usually zero 
492 
tend=50 * yr #  time to end simulation 
493 
outputs = 200 # number of time steps required. 
494 
h=(tendt)/outputs #size of time step 
495 
#user warning statement 
496 
print "Expected Number of time outputs is: ", (tendt)/h 
497 
i=0 #loop counter 
498 
\end{python} 
499 
Now that we know our inputs we will build a domain using the 
500 
\verbRectangle() function from \FINLEY. The four arguments allow us to 
501 
define our domain \verbmodel as: 
502 
\begin{python} 
503 
#generate domain using rectangle 
504 
blocks = Rectangle(l0=mx,l1=my,n0=ndx, n1=ndy) 
505 
\end{python} 
506 
\verbblocks now describes a domain in the manner of Section \ref{ss:domcon}. 
507 

508 
With a domain and all the required variables established, it is now possible to 
509 
set up our PDE so that it can be solved by \esc. The first step is to define the 
510 
type of PDE that we are trying to solve in each time step. In this example it is 
511 
a single linear PDE\footnote{in contrast to a system of PDEs which we discuss 
512 
later.}. We also need to state the values of our general form variables. 
513 
\begin{python} 
514 
mypde=LinearPDE(blocks) 
515 
A=zeros((2,2))) 
516 
A[0,0]=kappa 
517 
mypde.setValue(A=A, D=rhocp/h) 
518 
\end{python} 
519 
In many cases it may be possible to decrease the computational time of the 
520 
solver if the PDE is symmetric. 
521 
Symmetry of a PDE is defined by; 
522 
\begin{equation}\label{eqn:symm} 
523 
A_{jl}=A_{lj} 
524 
\end{equation} 
525 
Symmetry is only dependent on the $A$ coefficient in the general form and the 
526 
other coefficients $D$ as well as the right hand side $Y$. From the above 
527 
definition we can see that our PDE is symmetric. The \verbLinearPDE class 
528 
provides the method \method{checkSymmetry} to check if the given PDE is 
529 
symmetric. As our PDE is symmetrical we enable symmetry via; 
530 
\begin{python} 
531 
myPDE.setSymmetryOn() 
532 
\end{python} 
533 
Next we need to establish the initial temperature distribution \verbT. We need 
534 
to 
535 
assign the value \verbT1 to all sample points left to the contact interface at 
536 
$x_{0}=\frac{mx}{2}$ 
537 
and the value \verbT2 right to the contact interface. \esc 
538 
provides the \verbwhereNegative function to construct this. More 
539 
specifically, \verbwhereNegative returns the value $1$ at those sample points 
540 
where the argument has a negative value. Otherwise zero is returned. 
541 
If \verbx are the $x_{0}$ 
542 
coordinates of the sample points used to represent the temperature distribution 
543 
then \verbx[0]boundloc gives us a negative value for 
544 
all sample points left to the interface and nonnegative value to 
545 
the right of the interface. So with; 
546 
\begin{python} 
547 
# ... set initial temperature .... 
548 
T= T1*whereNegative(x[0]boundloc)+T2*(1whereNegative(x[0]boundloc)) 
549 
\end{python} 
550 
we get the desired temperature distribution. To get the actual sample points 
551 
\verbx we use the \verbgetX() method of the function space 
552 
\verbSolution(blocks) which is used to represent the solution of a PDE; 
553 
\begin{python} 
554 
x=Solution(blocks).getX() 
555 
\end{python} 
556 
As \verbx are the sample points for the function space 
557 
\verbSolution(blocks) 
558 
the initial temperature \verbT is using these sample points for 
559 
representation. 
560 
Although \esc is trying to be forgiving with the choice of sample points and to 
561 
convert 
562 
where necessary the adjustment of the function space is not always possible. So 
563 
it is advisable to make a careful choice on the function space used. 
564 

565 
Finally we initialise an iteration loop to solve our PDE for all the time steps 
566 
we specified in the variable section. As the right hand side of the general form 
567 
is dependent on the previous values for temperature \verb T across the bar this 
568 
must be updated in the loop. Our output at each time step is \verb T the heat 
569 
distribution and \verb totT the total heat in the system. 
570 
\begin{python} 
571 
while t < tend: 
572 
i+=1 #increment the counter 
573 
t+=h #increment the current time 
574 
mypde.setValue(Y=qH+rhocp/h*T) # set variable PDE coefficients 
575 
T=mypde.getSolution() #get the PDE solution 
576 
totE = integrate(rhocp*T) #get the total heat (energy) in the system 
577 
\end{python} 
578 
The last statement in this script calculates the total energy in the system as 
579 
the volume integral of $\rho c_{p} T$ over the block. 
580 
As the blocks are insulated no energy should be lost or added. 
581 
The total energy should stay constant for the example discussed here. 
582 

583 
\subsection{Running the Script} 
584 
The script presented so far is available under 
585 
\verbexample01a.py. You can edit this file with your favourite text editor. 
586 
On most operating systems\footnote{The \texttt{runescript} launcher is not 
587 
supported under {\it MS Windows} yet.} you can use the 
588 
\program{runescript} command 
589 
to launch {\it escript} scripts. For the example script use; 
590 
\begin{verbatim} 
591 
runescript example01a.py 
592 
\end{verbatim} 
593 
The program will print a progress report. Alternatively, you can use 
594 
the python interpreter directly; 
595 
\begin{verbatim} 
596 
python example01a.py 
597 
\end{verbatim} 
598 
if the system is configured correctly (please talk to your system 
599 
administrator). 
600 

601 
\subsection{Plotting the Total Energy} 
602 
\sslist{example01b.py} 
603 

604 
\esc does not include its own plotting capabilities. However, it is possible to 
605 
use a variety of free \pyt packages for visualisation. 
606 
Two types will be demonstrated in this cookbook; 
607 
\mpl\footnote{\url{http://matplotlib.sourceforge.net/}} and 
608 
\verbVTK\footnote{\url{http://www.vtk.org/}}. 
609 
The \mpl package is a component of SciPy\footnote{\url{http://www.scipy.org}} 
610 
and is good for basic graphs and plots. 
611 
For more complex visualisation tasks, in particular two and three dimensional 
612 
problems we recommend the use of more advanced tools. For instance, \mayavi 
613 
\footnote{\url{http://code.enthought.com/projects/mayavi/}} 
614 
which is based upon the \verbVTK toolkit. The usage of \verbVTK based 
615 
visualisation is discussed in Chapter~\ref{Sec:2DHD} which focuses on a two 
616 
dimensional PDE. 
617 

618 
For our simple granite block problem, we have two plotting tasks. Firstly, we 
619 
are interested in showing the 
620 
behaviour of the total energy over time and secondly, how the temperature 
621 
distribution within the block is developing over time. 
622 
Let us start with the first task. 
623 

624 
The idea is to create a record of the time marks and the corresponding total 
625 
energies observed. 
626 
\pyt provides the concept of lists for this. Before 
627 
the time loop is opened we create empty lists for the time marks \verbt_list 
628 
and the total energies \verbE_list. 
629 
After the new temperature has been calculated by solving the PDE we append the 
630 
new time marker and the total energy value for that time 
631 
to the corresponding list using the \verbappend method. With these 
632 
modifications our script looks as follows: 
633 
\begin{python} 
634 
t_list=[] 
635 
E_list=[] 
636 
# ... start iteration: 
637 
while t<tend: 
638 
t+=h 
639 
mypde.setValue(Y=qH+rhocp/h*T) # set variable PDE coefficients 
640 
T=mypde.getSolution() #get the PDE solution 
641 
totE=integrate(rhocp*T) 
642 
t_list.append(t) # add current time mark to record 
643 
E_list.append(totE) # add current total energy to record 
644 
\end{python} 
645 
To plot $t$ over $totE$ we use \mpl a module contained within \pylab which needs 
646 
to be loaded before use; 
647 
\begin{python} 
648 
import pylab as pl # plotting package. 
649 
\end{python} 
650 
Here we are not using \verbfrom pylab import * in order to avoid name 
651 
clashes for function names within \esc. 
652 

653 
The following statements are added to the script after the time loop has been 
654 
completed; 
655 
\begin{python} 
656 
pl.plot(t_list,E_list) 
657 
pl.title("Total Energy") 
658 
pl.axis([0,max(t_list),0,max(E_list)*1.1]) 
659 
pl.savefig("totE.png") 
660 
\end{python} 
661 
The first statement hands over the time marks and corresponding total energies 
662 
to the plotter. 
663 
The second statement sets the title for the plot. The third statement 
664 
sets the axis ranges. In most cases these are set appropriately by the plotter. 
665 

666 
The last statement generates the plot and writes the result into the file 
667 
\verbtotE.png which can be displayed by (almost) any image viewer. 
668 
As expected the total energy is constant over time, see 
669 
\reffig{fig:onedheatout1}. 
670 

671 
\begin{figure}[ht] 
672 
\begin{center} 
673 
\includegraphics[width=4in]{figures/ttblockspyplot150} 
674 
\caption{Example 1b: Total Energy in the Blocks over Time (in seconds)} 
675 
\label{fig:onedheatout1} 
676 
\end{center} 
677 
\end{figure} 
678 
\clearpage 
679 

680 
\subsection{Plotting the Temperature Distribution} 
681 
\label{sec: plot T} 
682 
\sslist{example01c.py} 
683 
For plotting the spatial distribution of the temperature we need to modify the 
684 
strategy we have used for the total energy. 
685 
Instead of producing a final plot at the end we will generate a 
686 
picture at each time step which can be browsed as a slide show or composed into 
687 
a movie. 
688 
The first problem we encounter is that if we produce an image at each time step 
689 
we need to make sure that the images previously generated are not overwritten. 
690 

691 
To develop an incrementing file name we can use the following convention. It is 
692 
convenient to put all image files showing the same variable  in our case the 
693 
temperature distribution  into a separate directory. 
694 
As part of the \verbos module\footnote{The \texttt{os} module provides 
695 
a powerful interface to interact with the operating system, see 
696 
\url{http://docs.python.org/library/os.html}.} \pyt 
697 
provides the \verbos.path.join command to build file and directory names in a 
698 
platform independent way. Assuming that 
699 
\verbsave_path is the name of the directory we want to put the results in the 
700 
command is; 
701 
\begin{python} 
702 
import os 
703 
os.path.join(save_path, "tempT%03d.png"%i ) 
704 
\end{python} 
705 
where \verbi is the time step counter. 
706 
There are two arguments to the \verbjoin command. The \verbsave_path 
707 
variable is a predefined string pointing to the directory we want to save our 
708 
data, for example a single subfolder called \verbdata would be defined by; 
709 
\begin{verbatim} 
710 
save_path = "data" 
711 
\end{verbatim} 
712 
while a subfolder of \verbdata called \verbexample01 would be defined by; 
713 
\begin{verbatim} 
714 
save_path = os.path.join("data","example01") 
715 
\end{verbatim} 
716 
The second argument of \verbjoin contains a string which is the file 
717 
name or subdirectory name. We can use the operator \verb% to use the value of 
718 
\verbi as part of our filename. The substring \verb%03d indicates that we 
719 
want to substitute a value into the name; 
720 
\begin{itemize} 
721 
\item \verb 0 means that small numbers should have leading zeroes; 
722 
\item \verb 3 means that numbers should be written using at least 3 digits; 
723 
and 
724 
\item \verb d means that the value to substitute will be a decimal integer. 
725 
\end{itemize} 
726 

727 
To actually substitute the value of \verbi into the name write \verb%i after 
728 
the string. 
729 
When done correctly, the output files from this command will be placed in the 
730 
directory defined by \verb save_path as; 
731 
\begin{verbatim} 
732 
blockspyplot001.png 
733 
blockspyplot002.png 
734 
blockspyplot003.png 
735 
... 
736 
\end{verbatim} 
737 
and so on. 
738 

739 
A subfolder check/constructor is available in \esc. The command; 
740 
\begin{verbatim} 
741 
mkDir(save_path) 
742 
\end{verbatim} 
743 
will check for the existence of \verb save_path and if missing, create the 
744 
required directories. 
745 

746 
We start by modifying our solution script. 
747 
Prior to the \verbwhile loop we need to extract our finite solution 
748 
points to a data object that is compatible with \mpl. First we create the node 
749 
coordinates of the sample points used to represent 
750 
the temperature as a \pyt list of tuples or a \numpy array as requested by the 
751 
plotting function. 
752 
We need to convert the array \verbx previously set as 
753 
\verbSolution(blocks).getX() into a \pyt list 
754 
and then to a \numpy array. The $x_{0}$ component is then extracted via 
755 
an array slice to the variable \verbplx; 
756 
\begin{python} 
757 
import numpy as np # array package. 
758 
#convert solution points for plotting 
759 
plx = x.toListOfTuples() 
760 
plx = np.array(plx) # convert to tuple to numpy array 
761 
plx = plx[:,0] # extract x locations 
762 
\end{python} 
763 

764 
\begin{figure} 
765 
\begin{center} 
766 
\includegraphics[width=4in]{figures/blockspyplot001} 
767 
\includegraphics[width=4in]{figures/blockspyplot050} 
768 
\includegraphics[width=4in]{figures/blockspyplot200} 
769 
\caption{Example 1c: Temperature ($T$) distribution in the blocks at time steps 
770 
$1$, $50$ and $200$} 
771 
\label{fig:onedheatout} 
772 
\end{center} 
773 
\end{figure} 
774 
\clearpage 
775 

776 
We use the same techniques provided by \mpl as we have used to plot the total 
777 
energy over time. 
778 
For each time step we generate a plot of the temperature distribution and save 
779 
each to a file. 
780 
The following is appended to the end of the \verbwhile loop and creates one 
781 
figure of the temperature distribution. We start by converting the solution to a 
782 
tuple and then plotting this against our \textit{x coordinates} \verbplx we 
783 
have generated before. We add a title to the diagram before it is rendered into 
784 
a file. 
785 
Finally, the figure is saved to a \verb*.png file and cleared for the 
786 
following iteration. 
787 
\begin{python} 
788 
# ... start iteration: 
789 
while t<tend: 
790 
i+=1 
791 
t+=h 
792 
mypde.setValue(Y=qH+rhocp/h*T) 
793 
T=mypde.getSolution() 
794 
totE=integrate(rhocp*T) 
795 
print "time step %s at t=%e days completed. total energy = %e."%(i,t/day,totE) 
796 
t_list.append(t) 
797 
E_list.append(totE) 
798 

799 
#establish figure 1 for temperature vs x plots 
800 
tempT = T.toListOfTuples() 
801 
pl.figure(1) #current figure 
802 
pl.plot(plx,tempT) #plot solution 
803 
# add title 
804 
pl.axis([0,mx,T1*.9,T2*1.1]) 
805 
pl.title("Temperature across blocks at time %d days"%(t/day)) 
806 
#save figure to file 
807 
pl.savefig(os.path.join(save_path,"tempT", "blockspyplot%03d.png"%i)) 
808 
pl.clf() #clear figure 
809 
\end{python} 
810 
Some results are shown in \reffig{fig:onedheatout}. 
811 

812 
\subsection{Making a Video} 
813 
Our saved plots from the previous section can be cast into a video using the 
814 
following command appended to the end of the script. The \verb mencoder command 
815 
is not available on every platform, so some users need to use an alternative 
816 
video encoder. 
817 
\begin{python} 
818 
# compile the *.png files to create a *.avi video that shows T change 
819 
# with time. This operation uses Linux mencoder. For other operating 
820 
# systems it is possible to use your favourite video compiler to 
821 
# convert image files to videos. 
822 

823 
os.system("mencoder mf://"+save_path+"/tempT"+"/*.png mf type=png:\ 
824 
w=800:h=600:fps=25 ovc lavc lavcopts vcodec=mpeg4 oac copy o \ 
825 
example01tempT.avi") 
826 
\end{python} 
827 
