13 


14 
\begin{figure}[h!] 
\begin{figure}[h!] 
15 
\centerline{\includegraphics[width=4.in]{figures/onedheatdiff001}} 
\centerline{\includegraphics[width=4.in]{figures/onedheatdiff001}} 
16 
\caption{Example 1: Temperature differential along a single interface between two granite blocks.} 
\caption{Example 1: Temperature differential along a single interface between 
17 

two granite blocks.} 
18 
\label{fig:onedgbmodel} 
\label{fig:onedgbmodel} 
19 
\end{figure} 
\end{figure} 
20 


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


24 
The first model consists of two blocks of isotropic material, for instance granite, sitting next to each other. 
The first model consists of two blocks of isotropic material, for instance 
25 
Initial temperature in \textit{Block 1} is \verbT1 and in \textit{Block 2} is \verbT2. 
granite, sitting next to each other. 
26 

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

\verbT2. 
28 
We assume that the system is insulated. 
We assume that the system is insulated. 
29 
What would happen to the temperature distribution in each block over time? 
What would happen to the temperature distribution in each block over time? 
30 
Intuition tells us that heat will be transported from the hotter block to the cooler one until both 
Intuition tells us that heat will be transported from the hotter block to the 
31 

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


34 
\subsection{1D Heat Diffusion Equation} 
\subsection{1D Heat Diffusion Equation} 
35 
We can model the heat distribution of this problem over time using one dimensional heat diffusion equation\footnote{A detailed discussion on how the heat diffusion equation is derived can be found at \url{http://online.redwoods.edu/instruct/darnold/DEProj/sp02/AbeRichards/paper.pdf}}; 
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 
which is defined as: 
which is defined as: 
41 
\begin{equation} 
\begin{equation} 
42 
\rho c\hackscore p \frac{\partial T}{\partial t}  \kappa \frac{\partial^{2} T}{\partial x^{2}} = q\hackscore H 
\rho c\hackscore p \frac{\partial T}{\partial t}  \kappa \frac{\partial^{2} 
43 

T}{\partial x^{2}} = q\hackscore H 
44 
\label{eqn:hd} 
\label{eqn:hd} 
45 
\end{equation} 
\end{equation} 
46 
where $\rho$ is the material density, $c\hackscore p$ is the specific heat and $\kappa$ is the thermal 
where $\rho$ is the material density, $c\hackscore p$ is the specific heat and 
47 
conductivity\footnote{A list of some common thermal conductivities is available from Wikipedia \url{http://en.wikipedia.org/wiki/List_of_thermal_conductivities}}. Here we assume that these material 
$\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}. 
parameters are \textbf{constant}. 
53 
The heat source is defined by the right hand side of \refEq{eqn:hd} as $q\hackscore{H}$; this can take the form of a constant or a function of time and space. For example $q\hackscore{H} = q\hackscore{0}e^{\gamma t}$ where we have the output of our heat source decaying with time. There are also two partial derivatives in \refEq{eqn:hd}; $\frac{\partial T}{\partial t}$ describes the change in temperature with time while $\frac{\partial ^2 T}{\partial x^2}$ is the spatial change of temperature. As there is only a single spatial dimension to our problem, our temperature solution $T$ is only dependent on the time $t$ and our signed distance from the the blockblock interface $x$. 
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 


63 
\subsection{PDEs and the General Form} 
\subsection{PDEs and the General Form} 
64 
It is possible to solve PDE \refEq{eqn:hd} analytically and obtain an exact solution to our problem. However, it is not always possible or practical to solve the problem this way. Alternatively, computers can be used to find the solution. To do this, a numerical approach is required to discretise 
It is possible to solve PDE \refEq{eqn:hd} analytically and obtain an exact 
65 
the PDE \refEq{eqn:hd} across time and space, this reduces the problem to a finite number of equations for a finite number of spatial points and time steps. These parameters together define the model. While discretisation introduces approximations and a degree of error, a sufficiently sampled model is generally accurate enough to satisfy the accuracy requirements for the final solution. 
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 
Firstly, we discretise the PDE \refEq{eqn:hd} in time. This leaves us with a steady linear PDE which involves spatial derivatives only and needs to be solved in each time step to progress in time. \esc can help us here. 
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 
For time discretization we use the Backwards Euler approximation scheme\footnote{see \url{http://en.wikipedia.org/wiki/Euler_method}}. It is based on the 
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 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 
approximation 
approximation 
82 
\begin{equation} 
\begin{equation} 
83 
\frac{\partial T(t)}{\partial t} \approx \frac{T(t)T(th)}{h} 
\frac{\partial T(t)}{\partial t} \approx \frac{T(t)T(th)}{h} 
89 
\frac{\partial T}{\partial t}(t^{(n)}) \approx \frac{T^{(n)}  T^{(n1)}}{h} 
\frac{\partial T}{\partial t}(t^{(n)}) \approx \frac{T^{(n)}  T^{(n1)}}{h} 
90 
\label{eqn:Tbeuler} 
\label{eqn:Tbeuler} 
91 
\end{equation} 
\end{equation} 
92 
where the upper index $n$ denotes the n\textsuperscript{th} time step. So one has 
where the upper index $n$ denotes the n\textsuperscript{th} time step. So one 
93 

has 
94 
\begin{equation} 
\begin{equation} 
95 
\begin{array}{rcl} 
\begin{array}{rcl} 
96 
t^{(n)} & = & t^{(n1)}+h \\ 
t^{(n)} & = & t^{(n1)}+h \\ 
100 
\end{equation} 
\end{equation} 
101 
Substituting \refEq{eqn:Tbeuler} into \refEq{eqn:hd} we get; 
Substituting \refEq{eqn:Tbeuler} into \refEq{eqn:hd} we get; 
102 
\begin{equation} 
\begin{equation} 
103 
\frac{\rho c\hackscore p}{h} (T^{(n)}  T^{(n1)})  \kappa \frac{\partial^{2} T^{(n)}}{\partial x^{2}} = q\hackscore H 
\frac{\rho c\hackscore p}{h} (T^{(n)}  T^{(n1)})  \kappa \frac{\partial^{2} 
104 

T^{(n)}}{\partial x^{2}} = q\hackscore H 
105 
\label{eqn:hddisc} 
\label{eqn:hddisc} 
106 
\end{equation} 
\end{equation} 
107 
Notice that we evaluate the spatial derivative term at current time $t^{(n)}$  therefore the name \textbf{backward Euler} scheme. Alternatively, one can evaluate the spatial derivative term at the previous time $t^{(n1)}$. This 
Notice that we evaluate the spatial derivative term at the current time 
108 
approach is called the \textbf{forward Euler} scheme. This scheme can provide some computational advantages, which 
$t^{(n)}$  therefore the name \textbf{backward Euler} scheme. Alternatively, 
109 
are not discussed here. However, the \textbf{forward Euler} scheme has a major disadvantage. Depending on the 
one can evaluate the spatial derivative term at the previous time $t^{(n1)}$. 
110 
material parameter as well as the discretization of the spatial derivative term, the time step size $h$ needs to be chosen sufficiently small to achieve a stable temperature when progressing in time. Stabiliy is achieved if the temperature does not grow beyond its initial bounds and become nonphysical. 
This 
111 
The backward Euler scheme, which we use here, is unconditionally stable meaning that under the assumption of 
approach is called the \textbf{forward Euler} scheme. This scheme can provide 
112 
physically correct problem setup the temperature approximation remains physical for all time steps. 
some computational advantages, which 
113 
The user needs to keep in mind that the discretization error introduced by \refEq{eqn:beuler} 
are not discussed here. However, the \textbf{forward Euler} scheme has a major 
114 
is sufficiently small, thus a good approximation of the true temperature is computed. It is 
disadvantage. Namely, depending on the 
115 
therefore crucial that the user remains sceptical about their results and for instance compares 
material parameters as well as the domain discretization of the spatial 
116 
the results for different time and spatial step sizes for correlation. 
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 


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

derivatives. To solve this problem 
134 
we want to use \esc. 
we want to use \esc. 
135 


136 
In \esc any given PDE can be described by the general form. For the purpose of this introduction we illustrate a simpler version of the general form for full linear PDEs which is available in the \esc user's guide. A simplified form that suits our heat diffusion problem\footnote{The form in the \esc users guide which uses the Einstein convention is written as 
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 
$(A\hackscore{jl} u\hackscore{,l})\hackscore{,j}+D u =Y$} 
$(A\hackscore{jl} u\hackscore{,l})\hackscore{,j}+D u =Y$} 
142 
is described by; 
is described by; 
143 
\begin{equation}\label{eqn:commonform nabla} 
\begin{equation}\label{eqn:commonform nabla} 
144 
\nabla\cdot(A\cdot\nabla u) + Du = f 
\nabla\cdot(A\cdot\nabla u) + Du = f 
145 
\end{equation} 
\end{equation} 
146 
where $A$, $D$ and $f$ are known values and $u$ is the unknown solution. The symbol $\nabla$ which is called the \textit{Nabla operator} or \textit{del operator} represents 
where $A$, $D$ and $f$ are known values and $u$ is the unknown solution. The 
147 
the spatial derivative of its subject  in this case $u$. Lets assume for a moment that we deal with a onedimensional problem then ; 
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 
\begin{equation} 
\begin{equation} 
152 
\nabla = \frac{\partial}{\partial x} 
\nabla = \frac{\partial}{\partial x} 
153 
\end{equation} 
\end{equation} 
155 
\begin{equation}\label{eqn:commonform} 
\begin{equation}\label{eqn:commonform} 
156 
A\frac{\partial^{2}u}{\partial x^{2}} + Du = f 
A\frac{\partial^{2}u}{\partial x^{2}} + Du = f 
157 
\end{equation} 
\end{equation} 
158 
if $A$ is constant. To match this simplified general form to our problem \refEq{eqn:hddisc} 
if $A$ is constant. To match this simplified general form to our problem 
159 

\refEq{eqn:hddisc} 
160 
we rearrange \refEq{eqn:hddisc}; 
we rearrange \refEq{eqn:hddisc}; 
161 
\begin{equation} 
\begin{equation} 
162 
\frac{\rho c\hackscore p}{h} T^{(n)}  \kappa \frac{\partial^2 T^{(n)}}{\partial x^2} = q\hackscore H + \frac{\rho c\hackscore p}{h} T^{(n1)} 
\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 
\label{eqn:hdgenf} 
\label{eqn:hdgenf} 
165 
\end{equation} 
\end{equation} 
166 
The PDE is now in a form that satisfies \refEq{eqn:commonform nabla} which is required for \esc to solve our PDE. This can be done by generating a solution for successive increments in the time nodes $t^{(n)}$ where 
The PDE is now in a form that satisfies \refEq{eqn:commonform nabla} which is 
167 
$t^{(0)}=0$ and $t^{(n)}=t^{(n1)}+h$ where $h>0$ is the step size and assumed to be constant. 
required for \esc to solve our PDE. This can be done by generating a solution 
168 
In the following the upper index ${(n)}$ refers to a value at time $t^{(n)}$. Finally, by comparing \refEq{eqn:hdgenf} with \refEq{eqn:commonform} one can see that; 
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 
\begin{equation}\label{ESCRIPT SET} 
\begin{equation}\label{ESCRIPT SET} 
175 
u=T^{(n)}; 
u=T^{(n)}; 
176 
A = \kappa; D = \frac{\rho c \hackscore{p}}{h}; f = q \hackscore{H} + \frac{\rho c\hackscore p}{h} T^{(n1)} 
A = \kappa; D = \frac{\rho c \hackscore{p}}{h}; f = q \hackscore{H} + \frac{\rho 
177 

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


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


194 
However, for this example we have made the model assumption that the system is insulated, so we need 
However, for this example we have made the model assumption that the system is 
195 

insulated, so we need 
196 
to add an appropriate boundary condition to prevent 
to add an appropriate boundary condition to prevent 
197 
any loss or inflow of energy at the boundary of our domain. Mathematically this is expressed by prescribing 
any loss or inflow of energy at the boundary of our domain. Mathematically this 
198 
the heat flux $\kappa \frac{\partial T}{\partial x}$ to zero. In our simplified one dimensional model this is expressed 
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 
in the form; 
in the form; 
202 
\begin{equation} 
\begin{equation} 
203 
\kappa \frac{\partial T}{\partial x} = 0 
\kappa \frac{\partial T}{\partial x} = 0 
206 
\begin{equation}\label{NEUMAN 1} 
\begin{equation}\label{NEUMAN 1} 
207 
\kappa \nabla T \cdot n = 0 
\kappa \nabla T \cdot n = 0 
208 
\end{equation} 
\end{equation} 
209 
where $n$ is the outer normal field \index{outer normal field} at the surface of the domain. 
where $n$ is the outer normal field \index{outer normal field} at the surface 
210 
The $\cdot$ (dot) refers to the dot product of the vectors $\nabla T$ and $n$. In fact, the term $\nabla T \cdot n$ is the normal derivative of 
of the domain. 
211 
the temperature $T$. Other notations used here are\footnote{The \esc notation for the normal 
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 
derivative is $T\hackscore{,i} n\hackscore i$.}; 
derivative is $T\hackscore{,i} n\hackscore i$.}; 
216 
\begin{equation} 
\begin{equation} 
217 
\nabla T \cdot n = \frac{\partial T}{\partial n} \; . 
\nabla T \cdot n = \frac{\partial T}{\partial n} \; . 
218 
\end{equation} 
\end{equation} 
219 
A condition of the type \refEq{NEUMAN 1} defines a \textbf{Neuman boundary condition} for the PDE. 
A condition of the type \refEq{NEUMAN 1} defines a \textbf{Neuman boundary 
220 

condition} for the PDE. 
221 


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

we use the term partial differential equation but in fact it is a boundary value 
230 
Conveniently, \esc makes default assumption on the boundary conditions which the user may modify where appropriate. 
problem. 
231 
For a problem of the form in~\refEq{eqn:commonform nabla} the default condition\footnote{In the form of the \esc users guide which is using the Einstein convention is written as 
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 


236 

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 
$n\hackscore{j}A\hackscore{jl} u\hackscore{,l}=0$.} is; 
$n\hackscore{j}A\hackscore{jl} u\hackscore{,l}=0$.} is; 
242 
\begin{equation}\label{NEUMAN 2} 
\begin{equation}\label{NEUMAN 2} 
243 
n\cdot A \cdot\nabla u = 0 
n\cdot A \cdot\nabla u = 0 
244 
\end{equation} 
\end{equation} 
245 
which is used everywhere on the boundary. Again $n$ denotes the outer normal field. 
which is used everywhere on the boundary. Again $n$ denotes the outer normal 
246 
Notice that the coefficient $A$ is the same as in the \esc PDE~\ref{eqn:commonform nabla}. 
field. 
247 
With the settings for the coefficients we have already identified in \refEq{ESCRIPT SET} this 
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 
condition translates into 
condition translates into 
252 
\begin{equation}\label{NEUMAN 2b} 
\begin{equation}\label{NEUMAN 2b} 
253 
\kappa \frac{\partial T}{\partial x} = 0 
\kappa \frac{\partial T}{\partial x} = 0 
254 
\end{equation} 
\end{equation} 
255 
for the boundary of the domain. This is identical to the Neuman boundary condition we want to set. \esc will take care of this condition for us. We discuss the Dirichlet boundary condition later. 
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 


259 
\subsection{Outline of the Implementation} 
\subsection{Outline of the Implementation} 
260 
\label{sec:outline} 
\label{sec:outline} 
261 
To solve the heat diffusion equation (\refEq{eqn:hd}) we write a simple \pyt script. At this point we assume that you have some basic understanding of the \pyt programming language. If not, there are some pointers and links available in Section \ref{sec:escpybas}. The script (discussed in \refSec{sec:key}) has four major steps. Firstly, we need to define the domain where we want to 
To solve the heat diffusion equation (\refEq{eqn:hd}) we write a simple \pyt 
262 
calculate the temperature. For our problem this is the joint blocks of granite which has a rectangular shape. Secondly, we need to define the PDE to solve in each time step to get the updated temperature. Thirdly, we need to define the coefficients of the PDE and finally we need to solve the PDE. The last two steps need to be repeated until the final time marker has been reached. The work flow described in \reffig{fig:wf}. 
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 
% \begin{enumerate} 
% \begin{enumerate} 
273 
% \item create domain 
% \item create domain 
274 
% \item create PDE 
% \item create PDE 
288 
\label{fig:wf} 
\label{fig:wf} 
289 
\end{figure} 
\end{figure} 
290 


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


302 


303 
\begin{figure}[t] 
\begin{figure}[t] 
309 


310 
\subsection{The Domain Constructor in \esc} 
\subsection{The Domain Constructor in \esc} 
311 
\label{ss:domcon} 
\label{ss:domcon} 
312 
While it is not strictly relevant or necessary, it can be helpful to have a better understanding of how values are spatially distributed and how PDE coefficients are interpreted in \esc. The example in this case would be temperature. 
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 
There are various ways to construct domain objects. The simplest form is a rectangular shaped region with a length and height. There is 
coefficients are interpreted in \esc can be helpful. 
315 
a ready to use function for this named \verb rectangle(). Besides the spatial dimensions this function requires to specify the number of 

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

Any spatially distributed value 
322 
In general, one can think about a domain object as a composition of nodes and elements. 
and the PDE is represented in discrete form using this element 
323 
As shown in \reffig{fig:fs}, an element is defined by the nodes that are used to describe its vertices. 
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 


334 

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 
To represent spatial distributed values the user can use 
To represent spatial distributed values the user can use 
339 
the values at the nodes, at the elements in the interior of the domain or at the elements located at the surface of the domain. 
the values at the nodes, at the elements in the interior of the domain or at the 
340 
The different approach used to represent values is called \textbf{function space} and is attached to all objects 
elements located at the surface of the domain. 
341 
in \esc representing a spatial distributed value such as the solution of a PDE. The three 
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 
function spaces we use at the moment are; 
function spaces we use at the moment are; 
346 
\begin{enumerate} 
\begin{enumerate} 
347 
\item the nodes, called by \verbContinuousFunction(domain) ; 
\item the nodes, called by \verbContinuousFunction(domain) ; 
348 
\item the elements/cells, called by \verbFunction(domain) ; and 
\item the elements/cells, called by \verbFunction(domain) ; and 
349 
\item the boundary, called by \verbFunctionOnBoundary(domain) . 
\item the boundary, called by \verbFunctionOnBoundary(domain) . 
350 
\end{enumerate} 
\end{enumerate} 
351 
A function space object such as \verbContinuousFunction(domain) has the method \verbgetX attached to it. This method returns the 
A function space object such as \verbContinuousFunction(domain) has the method 
352 
location of the socalled \textbf{sample points} used to represent values of the particular function space. So the 
\verbgetX attached to it. This method returns the 
353 
call \verbContinuousFunction(domain).getX() will return the coordinates of the nodes used to describe the domain while 
location of the socalled \textbf{sample points} used to represent values of the 
354 
the \verbFunction(domain).getX() returns the coordinates of numerical integration points within elements, see 
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 
\reffig{fig:fs}. 
\reffig{fig:fs}. 
360 


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


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


389 


390 
\subsection{A Clarification for the 1D Case} 
\subsection{A Clarification for the 1D Case} 
391 
\label{SEC: 1D CLARIFICATION} 
\label{SEC: 1D CLARIFICATION} 
392 
It is necessary for clarification that we revisit our general PDE from \refeq{eqn:commonform nabla} for two dimensional domain. \esc is inherently designed to solve problems that are greater than one dimension and so \refEq{eqn:commonform nabla} needs to be read as a higher dimensional problem. In the case of two spatial dimensions the \textit{Nabla operator} has in fact two components $\nabla = (\frac{\partial}{\partial x}, \frac{\partial}{\partial y})$. Assuming the coefficient $A$ is constant, the \refEq{eqn:commonform nabla} takes the following form; 
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 
\begin{equation}\label{eqn:commonform2D} 
\begin{equation}\label{eqn:commonform2D} 
401 
A\hackscore{00}\frac{\partial^{2}u}{\partial x^{2}} 
A\hackscore{00}\frac{\partial^{2}u}{\partial x^{2}} 
402 
A\hackscore{01}\frac{\partial^{2}u}{\partial x\partial y} 
A\hackscore{01}\frac{\partial^{2}u}{\partial x\partial y} 
417 
\subsection{Developing a PDE Solution Script} 
\subsection{Developing a PDE Solution Script} 
418 
\label{sec:key} 
\label{sec:key} 
419 
\sslist{example01a.py} 
\sslist{example01a.py} 
420 
We write a simple \pyt script which uses the \modescript, \modfinley and \modmpl modules. 
We write a simple \pyt script which uses the \modescript, \modfinley and \modmpl 
421 
By developing a script for \esc, the heat diffusion equation can be solved at successive time steps for a predefined period using our general form \refEq{eqn:hdgenf}. Firstly it is necessary to import all the libraries\footnote{The libraries contain predefined scripts that are required to solve certain problems, these can be simple like sine and cosine functions or more complicated like those from our \esc library.} 
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 
that we will require. 
that we will require. 
429 
\begin{python} 
\begin{python} 
430 
from esys.escript import * 
from esys.escript import * 
436 
# match up in the equations under SI. 
# match up in the equations under SI. 
437 
from esys.escript.unitsSI import * 
from esys.escript.unitsSI import * 
438 
\end{python} 
\end{python} 
439 
It is generally a good idea to import all of the \modescript library, although if the functions and classes required are known they can be specified individually. The function \verbLinearPDE has been imported explicitly for ease of use later in the script. \verbRectangle is going to be our type of domain. The module \verb unitsSI provides support for SI unit definitions with our variables. 
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 
Once our library dependencies have been established, defining the problem specific variables is the next step. In general the number of variables needed will vary between problems. These variables belong to two categories. They are either directly related to the PDE and can be used as inputs into the \esc solver, or they are script variables used to control internal functions and iterations in our problem. For this PDE there are a number of constants which need values. Firstly, the domain upon which we wish to solve our problem needs to be defined. There are different types of domains in \modescript which we demonstrate in later tutorials but for our granite blocks, we simply use a rectangular domain. 
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 
Using a rectangular domain simplifies our granite blocks (which would in reality be a \textit{3D} object) into a single dimension. The granite blocks will have a lengthways cross section that looks like a rectangle. As a result we do not need to model the volume of the block due to symmetry. There are four arguments we must consider when we decide to create a rectangular domain, the domain \textit{length}, \textit{width} and \textit{step size} in each direction. When defining the size of our problem it will help us determine appropriate values for our model arguments. If we make our dimensions large but our step sizes very small we increase the accuracy of our solution. Unfortunately we also increase the number of calculations that must be solved per time step. This means more computational time is required to produce a solution. In this \textit{1D} problem, the bar is defined as being 1 metre long. An appropriate step size \verbndx would be 1 to 10\% of the length. Our \verbndy need only be 1, this is because our problem stipulates no partial derivatives in the $y$ direction. Thus the temperature does not vary with $y$. Hence, the model parameters can be defined as follows; note we have used the \verb unitsSI convention to make sure all our input units are converted to SI. 
domain. The module \verb unitsSI provides support for SI unit definitions with 
444 

our variables. 
445 


446 

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 


457 

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 
\begin{python} 
\begin{python} 
475 
mx = 500.*m #meters  model length 
mx = 500.*m #meters  model length 
476 
my = 100.*m #meters  model width 
my = 100.*m #meters  model width 
478 
ndy = 1 # mesh steps in y direction 
ndy = 1 # mesh steps in y direction 
479 
boundloc = mx/2 # location of boundary between the two blocks 
boundloc = mx/2 # location of boundary between the two blocks 
480 
\end{python} 
\end{python} 
481 
The material constants and the temperature variables must also be defined. For the granite in the model they are defined as: 
The material constants and the temperature variables must also be defined. For 
482 

the granite in the model they are defined as: 
483 
\begin{python} 
\begin{python} 
484 
#PDE related 
#PDE related 
485 
rho = 2750. *kg/m**3 #kg/m^{3} density of iron 
rho = 2750. *kg/m**3 #kg/m^{3} density of iron 
490 
T1=20 * Celsius # initial temperature at Block 1 
T1=20 * Celsius # initial temperature at Block 1 
491 
T2=2273. * Celsius # base temperature at Block 2 
T2=2273. * Celsius # base temperature at Block 2 
492 
\end{python} 
\end{python} 
493 
Finally, to control our script we will have to specify our timing controls and where we would like to save the output from the solver. This is simple enough: 
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 
\begin{python} 
\begin{python} 
496 
t=0 * day #our start time, usually zero 
t=0 * day #our start time, usually zero 
497 
tend=1. * day #  time to end simulation 
tend=1. * day #  time to end simulation 
501 
print "Expected Number of time outputs is: ", (tendt)/h 
print "Expected Number of time outputs is: ", (tendt)/h 
502 
i=0 #loop counter 
i=0 #loop counter 
503 
\end{python} 
\end{python} 
504 
Now that we know our inputs we will build a domain using the \verb Rectangle() function from \verb finley . The four arguments allow us to define our domain \verb model as: 
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 
\begin{python} 
\begin{python} 
508 
#generate domain using rectangle 
#generate domain using rectangle 
509 
blocks = Rectangle(l0=mx,l1=my,n0=ndx, n1=ndy) 
blocks = Rectangle(l0=mx,l1=my,n0=ndx, n1=ndy) 
510 
\end{python} 
\end{python} 
511 
\verb blocks now describes a domain in the manner of Section \ref{ss:domcon}. 
\verb blocks now describes a domain in the manner of Section \ref{ss:domcon}. 
512 


513 
With a domain and all our required variables established, it is now possible to set up our PDE so that it can be solved by \esc. The first step is to define the type of PDE that we are trying to solve in each time step. In this example it is a single linear PDE\footnote{in comparison to a system of PDEs which we discuss later.}. We also need to state the values of our general form variables. 
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 
\begin{python} 
\begin{python} 
519 
mypde=LinearPDE(blocks) 
mypde=LinearPDE(blocks) 
520 
A=zeros((2,2))) 
A=zeros((2,2))) 
521 
A[0,0]=kappa 
A[0,0]=kappa 
522 
mypde.setValue(A=A, D=rhocp/h) 
mypde.setValue(A=A, D=rhocp/h) 
523 
\end{python} 
\end{python} 
524 
In a many cases it may be possible to decrease the computational time of the solver if the PDE is symmetric. 
In a many cases it may be possible to decrease the computational time of the 
525 

solver if the PDE is symmetric. 
526 
Symmetry of a PDE is defined by; 
Symmetry of a PDE is defined by; 
527 
\begin{equation}\label{eqn:symm} 
\begin{equation}\label{eqn:symm} 
528 
A\hackscore{jl}=A\hackscore{lj} 
A\hackscore{jl}=A\hackscore{lj} 
529 
\end{equation} 
\end{equation} 
530 
Symmetry is only dependent on the $A$ coefficient in the general form and the other coefficients $D$ as well as the right hand side $Y$. From the above definition we can see that our PDE is symmetric. The \verb LinearPDE class provides the method \method{checkSymmetry} to check if the given PDE is symmetric. As our PDE is symmetrical we enable symmetry via; 
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 
\begin{python} 
\begin{python} 
536 
myPDE.setSymmetryOn() 
myPDE.setSymmetryOn() 
537 
\end{python} 
\end{python} 
538 
Next we need to establish the initial temperature distribution \verbT. We need to 
Next we need to establish the initial temperature distribution \verbT. We need 
539 
assign the value \verbT1 to all sample points left to the contact interface at $x\hackscore{0}=\frac{mx}{2}$ 
to 
540 

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

$x\hackscore{0}=\frac{mx}{2}$ 
542 
and the value \verbT2 right to the contact interface. \esc 
and the value \verbT2 right to the contact interface. \esc 
543 
provides the \verbwhereNegative function to construct this. In fact, 
provides the \verbwhereNegative function to construct this. In fact, 
544 
\verbwhereNegative returns the value $1$ at those sample points where the argument 
\verbwhereNegative returns the value $1$ at those sample points where the 
545 
has a negative value. Otherwise zero is returned. If \verbx are the $x\hackscore{0}$ 
argument 
546 

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

$x\hackscore{0}$ 
548 
coordinates of the sample points used to represent the temperature distribution 
coordinates of the sample points used to represent the temperature distribution 
549 
then \verbx[0]boundloc gives us a negative value for 
then \verbx[0]boundloc gives us a negative value for 
550 
all sample points left to the interface and nonnegative value to 
all sample points left to the interface and nonnegative value to 
553 
# ... set initial temperature .... 
# ... set initial temperature .... 
554 
T= T1*whereNegative(x[0]boundloc)+T2*(1whereNegative(x[0]boundloc)) 
T= T1*whereNegative(x[0]boundloc)+T2*(1whereNegative(x[0]boundloc)) 
555 
\end{python} 
\end{python} 
556 
we get the desired temperature distribution. To get the actual sample points \verbx we use 
we get the desired temperature distribution. To get the actual sample points 
557 

\verbx we use 
558 
the \verbgetX() method of the function space \verbSolution(blocks) 
the \verbgetX() method of the function space \verbSolution(blocks) 
559 
which is used to represent the solution of a PDE; 
which is used to represent the solution of a PDE; 
560 
\begin{python} 
\begin{python} 
561 
x=Solution(blocks).getX() 
x=Solution(blocks).getX() 
562 
\end{python} 
\end{python} 
563 
As \verbx are the sample points for the function space \verbSolution(blocks) 
As \verbx are the sample points for the function space 
564 
the initial temperature \verbT is using these sample points for representation. 
\verbSolution(blocks) 
565 
Although \esc is trying to be forgiving with the choice of sample points and to convert 
the initial temperature \verbT is using these sample points for 
566 
where necessary the adjustment of the function space is not always possible. So it is 
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 
advisable to make a careful choice on the function space used. 
advisable to make a careful choice on the function space used. 
572 


573 
Finally we initialise an iteration loop to solve our PDE for all the time steps we specified in the variable section. As the right hand side of the general form is dependent on the previous values for temperature \verb T across the bar this must be updated in the loop. Our output at each time step is \verb T the heat distribution and \verb totT the total heat in the system. 
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 
\begin{python} 
\begin{python} 
579 
while t < tend: 
while t < tend: 
580 
i+=1 #increment the counter 
i+=1 #increment the counter 
583 
T=mypde.getSolution() #get the PDE solution 
T=mypde.getSolution() #get the PDE solution 
584 
totE = integrate(rhocp*T) #get the total heat (energy) in the system 
totE = integrate(rhocp*T) #get the total heat (energy) in the system 
585 
\end{python} 
\end{python} 
586 
The last statement in this script calculates the total energy in the system as volume integral 
The last statement in this script calculates the total energy in the system as 
587 
of $\rho c\hackscore{p} T$ over the block. As the blocks are insulated no energy should be get lost or added. 
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 
The total energy should stay constant for the example discussed here. 
The total energy should stay constant for the example discussed here. 
591 


592 
\subsection{Running the Script} 
\subsection{Running the Script} 
593 
The script presented so far is available under 
The script presented so far is available under 
594 
\verbexample01a.py. You can edit this file with your favourite text editor. 
\verbexample01a.py. You can edit this file with your favourite text editor. 
595 
On most operating systems\footnote{The you can use \texttt{runescript} launcher is not supported under {\it MS Windows} yet.} you can use the \program{runescript} command 
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 
to launch {\it escript} scripts. For the example script use; 
to launch {\it escript} scripts. For the example script use; 
599 
\begin{verbatim} 
\begin{verbatim} 
600 
runescript example01a.py 
runescript example01a.py 
604 
\begin{verbatim} 
\begin{verbatim} 
605 
python example01a.py 
python example01a.py 
606 
\end{verbatim} 
\end{verbatim} 
607 
if the system is configured correctly (Please talk to your system administrator). 
if the system is configured correctly (Please talk to your system 
608 

administrator). 
609 


610 
\begin{figure} 
\begin{figure} 
611 
\begin{center} 
\begin{center} 
618 
\subsection{Plotting the Total Energy} 
\subsection{Plotting the Total Energy} 
619 
\sslist{example01b.py} 
\sslist{example01b.py} 
620 


621 
\esc does not include its own plotting capabilities. However, it is possible to use a variety of free \pyt packages for visualisation. 
\esc does not include its own plotting capabilities. However, it is possible to 
622 
Two types will be demonstrated in this cookbook; \mpl\footnote{\url{http://matplotlib.sourceforge.net/}} and \verb VTK \footnote{\url{http://www.vtk.org/}} visualisation. 
use a variety of free \pyt packages for visualisation. 
623 
The \mpl package is a component of SciPy\footnote{\url{http://www.scipy.org}} and is good for basic graphs and plots. 
Two types will be demonstrated in this cookbook; 
624 
For more complex visualisation tasks in particular, two and three dimensional problems we recommend the use of more advanced tools. For instance, \mayavi \footnote{\url{http://code.enthought.com/projects/mayavi/}} 
\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 
which is based upon the \verbVTK toolkit. The usage of \verbVTK based 
which is based upon the \verbVTK toolkit. The usage of \verbVTK based 
632 
visualization is discussed in Chapter~\ref{Sec:2DHD} which focusses on a two dimensional PDE. 
visualization is discussed in Chapter~\ref{Sec:2DHD} which focusses on a two 
633 

dimensional PDE. 
634 


635 
For our simple granite block problem, we have two plotting tasks. Firstly, we are interested in showing the 
For our simple granite block problem, we have two plotting tasks. Firstly, we 
636 
behaviour of the total energy over time and secondly, how the temperature distribution within the block is 
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 
developing over time. Lets start with the first task. 
developing over time. Lets start with the first task. 
640 


641 
The trick is to create a record of the time marks and the corresponding total energies observed. 
The trick is to create a record of the time marks and the corresponding total 
642 

energies observed. 
643 
\pyt provides the concept of lists for this. Before 
\pyt provides the concept of lists for this. Before 
644 
the time loop is opened we create empty lists for the time marks \verbt_list and the total energies \verbE_list. 
the time loop is opened we create empty lists for the time marks \verbt_list 
645 
After the new temperature has been calculated by solving the PDE we append the new time marker and the total energy value for that time 
and the total energies \verbE_list. 
646 
to the corresponding list using the \verbappend method. With these modifications our script looks as follows: 
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 
\begin{python} 
\begin{python} 
651 
t_list=[] 
t_list=[] 
652 
E_list=[] 
E_list=[] 
659 
t_list.append(t) # add current time mark to record 
t_list.append(t) # add current time mark to record 
660 
E_list.append(totE) # add current total energy to record 
E_list.append(totE) # add current total energy to record 
661 
\end{python} 
\end{python} 
662 
To plot $t$ over $totE$ we use \mpl a module contained within \pylab which needs to be loaded before used; 
To plot $t$ over $totE$ we use \mpl a module contained within \pylab which needs 
663 

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

clashes for function names 
669 
within \esc. 
within \esc. 
670 


671 
The following statements are added to the script after the time loop has been completed; 
The following statements are added to the script after the time loop has been 
672 

completed; 
673 
\begin{python} 
\begin{python} 
674 
pl.plot(t_list,E_list) 
pl.plot(t_list,E_list) 
675 
pl.title("Total Energy") 
pl.title("Total Energy") 
676 
pl.axis([0,max(t_list),0,max(E_list)*1.1]) 
pl.axis([0,max(t_list),0,max(E_list)*1.1]) 
677 
pl.savefig("totE.png") 
pl.savefig("totE.png") 
678 
\end{python} 
\end{python} 
679 
The first statement hands over the time marks and corresponding total energies to the plotter. 
The first statement hands over the time marks and corresponding total energies 
680 

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


684 
The last statement renders the plot and writes the 
The last statement renders the plot and writes the 
685 
result into the file \verbtotE.png which can be displayed by (almost) any image viewer. 
result into the file \verbtotE.png which can be displayed by (almost) any 
686 
As expected the total energy is constant over time, see \reffig{fig:onedheatout1}. 
image viewer. 
687 

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

\reffig{fig:onedheatout1}. 
689 


690 
\subsection{Plotting the Temperature Distribution} 
\subsection{Plotting the Temperature Distribution} 
691 
\label{sec: plot T} 
\label{sec: plot T} 
692 
\sslist{example01c.py} 
\sslist{example01c.py} 
693 
For plotting the spatial distribution of the temperature we need to modify the strategy we have used 
For plotting the spatial distribution of the temperature we need to modify the 
694 
for the total energy. Instead of producing a final plot at the end we will generate a 
strategy we have used 
695 
picture at each time step which can be browsed as a slide show or composed into a movie. 
for the total energy. Instead of producing a final plot at the end we will 
696 
The first problem we encounter is that if we produce an image at each time step we need 
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 
to make sure that the images previously generated are not overwritten. 
to make sure that the images previously generated are not overwritten. 
702 


703 
To develop an incrementing file name we can use the following convention. It is convenient to 
To develop an incrementing file name we can use the following convention. It is 
704 
put all image file showing the same variable  in our case the temperature distribution  
convenient to 
705 
into a separate directory. As part of the \verbos module\footnote{The \texttt{os} module provides 
put all image file showing the same variable  in our case the temperature 
706 
a powerful interface to interact with the operating system, see \url{http://docs.python.org/library/os.html}.} \pyt 
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 
provides the \verbos.path.join command to build file and 
provides the \verbos.path.join command to build file and 
712 
directory names in a platform independent way. Assuming that 
directory names in a platform independent way. Assuming that 
713 
\verbsave_path is name of directory we want to put the results the command is; 
\verbsave_path is name of directory we want to put the results the command 
714 

is; 
715 
\begin{python} 
\begin{python} 
716 
import os 
import os 
717 
os.path.join(save_path, "tempT%03d.png"%i ) 
os.path.join(save_path, "tempT%03d.png"%i ) 
718 
\end{python} 
\end{python} 
719 
where \verbi is the time step counter. 
where \verbi is the time step counter. 
720 
There are two arguments to the \verb join command. The \verb save_path variable is a predefined string pointing to the directory we want to save our data, for example a single subfolder called \verb data would be defined by; 
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 
\begin{verbatim} 
\begin{verbatim} 
724 
save_path = "data" 
save_path = "data" 
725 
\end{verbatim} 
\end{verbatim} 
727 
\begin{verbatim} 
\begin{verbatim} 
728 
save_path = os.path.join("data","example01") 
save_path = os.path.join("data","example01") 
729 
\end{verbatim} 
\end{verbatim} 
730 
The second argument of \verb join \xspace contains a string which is the file name or subdirectory name. We can use the operator \verb% to use the value of \verbi as part of our filename. The substring \verb %03d indicates that we want to substitude a value into the name; 
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 
\begin{itemize} 
\begin{itemize} 
735 
\item \verb 0 means that small numbers should have leading zeroes; 
\item \verb 0 means that small numbers should have leading zeroes; 
736 
\item \verb 3 means that numbers should be written using at least 3 digits; and 
\item \verb 3 means that numbers should be written using at least 3 digits; 
737 

and 
738 
\item \verb d means that the value to substitute will be an integer. 
\item \verb d means that the value to substitute will be an integer. 
739 
\end{itemize} 
\end{itemize} 
740 
To actually substitute the value of \verbi into the name write \verb %i after the string. 
To actually substitute the value of \verbi into the name write \verb %i after 
741 
When done correctly, the output files from this command would be place in the directory defined by \verb save_path as; 
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 
\begin{verbatim} 
\begin{verbatim} 
745 
blockspyplot001.png 
blockspyplot001.png 
746 
blockspyplot002.png 
blockspyplot002.png 
753 
\begin{verbatim} 
\begin{verbatim} 
754 
mkDir(save_path) 
mkDir(save_path) 
755 
\end{verbatim} 
\end{verbatim} 
756 
will check for the existence of \verb save_path and if missing, make the required directories. 
will check for the existence of \verb save_path and if missing, make the 
757 

required directories. 
758 


759 
We start by modifying our solution script. 
We start by modifying our solution script. 
760 
Prior to the \verbwhile loop we will need to extract our finite solution points to a data object that is compatible with \mpl. First we create the node coordinates of the sample points used to represent 
Prior to the \verbwhile loop we will need to extract our finite solution 
761 
the temperature as a \pyt list of tuples or a \numpy array as requested by the plotting function. 
points to a data object that is compatible with \mpl. First we create the node 
762 
We need to convert the array \verbx previously set as \verbSolution(blocks).getX() into a \pyt list 
coordinates of the sample points used to represent 
763 
and then to a \numpy array. The $x\hackscore{0}$ component is then extracted via an array slice to the variable \verbplx; 
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 
\begin{python} 
\begin{python} 
770 
import numpy as np # array package. 
import numpy as np # array package. 
771 
#convert solution points for plotting 
#convert solution points for plotting 
779 
\includegraphics[width=4in]{figures/blockspyplot001} 
\includegraphics[width=4in]{figures/blockspyplot001} 
780 
\includegraphics[width=4in]{figures/blockspyplot050} 
\includegraphics[width=4in]{figures/blockspyplot050} 
781 
\includegraphics[width=4in]{figures/blockspyplot200} 
\includegraphics[width=4in]{figures/blockspyplot200} 
782 
\caption{Example 1c: Temperature ($T$) distribution in the blocks at time steps $1$, $50$ and $200$.} 
\caption{Example 1c: Temperature ($T$) distribution in the blocks at time steps 
783 

$1$, $50$ and $200$.} 
784 
\label{fig:onedheatout} 
\label{fig:onedheatout} 
785 
\end{center} 
\end{center} 
786 
\end{figure} 
\end{figure} 
787 


788 
We use the same techniques provided by \mpl as we have used to plot the total energy over time. 
We use the same techniques provided by \mpl as we have used to plot the total 
789 
For each time step we generate a plot of the temperature distribution and save each to a file. 
energy over time. 
790 
The following is appended to the end of the \verb while loop and creates one figure of the temperature distribution. We start by converting the solution to a tuple and then plotting this against our \textit{x coordinates} \verb plx we have generated before. We add a title to the diagram before it is rendered into a file. 
For each time step we generate a plot of the temperature distribution and save 
791 
Finally, the figure is saved to a \verb*.png file and cleared for the following iteration. 
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 
\begin{python} 
\begin{python} 
800 
# ... start iteration: 
# ... start iteration: 
801 
while t<tend: 
while t<tend: 
813 
Some results are shown in \reffig{fig:onedheatout}. 
Some results are shown in \reffig{fig:onedheatout}. 
814 


815 
\subsection{Make a video} 
\subsection{Make a video} 
816 
Our saved plots from the previous section can be cast into a video using the following command appended to the end of the script. The \verb mencoder command is Linux only, so other platform users need to use an alternative video encoder. 
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 
\begin{python} 
\begin{python} 
820 
# compile the *.png files to create a *.avi videos that show T change 
# compile the *.png files to create a *.avi videos that show T change 
821 
# with time. This operation uses Linux mencoder. For other operating 
# with time. This operation uses Linux mencoder. For other operating 