37 
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 $\kappa$ is the thermal 
38 
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 
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 
39 
parameters are \textbf{constant}. 
parameters are \textbf{constant}. 
40 
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 position along the iron bar $x$. 
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$. 
41 


42 
\subsection{PDEs and the General Form} 
\subsection{PDEs and the General Form} 
43 
Potentially, 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 a problem this way. Alternatively, computers can be used to solve these kinds of problems. To do this, a numerical approach is required to discretise 
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 
44 
the PDE \refEq{eqn:hd} in time and space, then we left with a finite number of equations for a finite number of spatial points and time steps in the model. While discretisation introduces approximations and a degree of error, a sufficiently sampled model is generally accurate enough to satisfy the requirements of the modeller. 
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. 
45 


46 
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. 
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. 
47 


72 
\end{equation} 
\end{equation} 
73 
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 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 
74 
approach is called the \textbf{forward Euler} scheme. This scheme can provide some computational advantages, which 
approach is called the \textbf{forward Euler} scheme. This scheme can provide some computational advantages, which 
75 
are not discussed here. However, this scheme has a major disadvantage, namely depending on the 
are not discussed here. However, the \textbf{forward Euler} scheme has a major disadvantage. Depending on the 
76 
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. The term \textit{stable} means 
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. 

that the approximation of the temperature will not grow beyond its initial bounds and become nonphysical. 

77 
The backward Euler scheme, which we use here, is unconditionally stable meaning that under the assumption of 
The backward Euler scheme, which we use here, is unconditionally stable meaning that under the assumption of 
78 
physically correct problem setup the temperature approximation remains physical for all time steps. 
physically correct problem setup the temperature approximation remains physical for all time steps. 
79 
The user needs to keep in mind that the discretization error introduced by \refEq{eqn:beuler} 
The user needs to keep in mind that the discretization error introduced by \refEq{eqn:beuler} 
80 
is sufficiently small, thus a good approximation of the true temperature is computed. It is 
is sufficiently small, thus a good approximation of the true temperature is computed. It is 
81 
therefore crucial that the user remains critical about his/her results and for instance compares 
therefore crucial that the user remains sceptical about their results and for instance compares 
82 
the results for different time and spatial step sizes. 
the results for different time and spatial step sizes for correlation. 
83 


84 
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 
85 
differential equation \refEq{eqn:hddisc} which only includes spatial derivatives. To solve this problem 
differential equation \refEq{eqn:hddisc} which only includes spatial derivatives. To solve this problem 
86 
we want to use \esc. 
we want to use \esc. 
87 


88 
In \esc any given PDE described by 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 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 
89 
$(A\hackscore{jl} u\hackscore{,l})\hackscore{,j}+D u =Y$} 
$(A\hackscore{jl} u\hackscore{,l})\hackscore{,j}+D u =Y$} 
90 
is described by; 
is described by; 
91 
\begin{equation}\label{eqn:commonform nabla} 
\begin{equation}\label{eqn:commonform nabla} 
120 
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. 
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. 
121 
We discuss Dirichlet boundary condition in our second example presented in Section~\ref{Sec:1DHDv0}. 
We discuss Dirichlet boundary condition in our second example presented in Section~\ref{Sec:1DHDv0}. 
122 


123 
We make 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 insulated, so we need 
124 
to add an appropriate boundary condition to prevent 
to add an appropriate boundary condition to prevent 
125 
any loss or inflow of energy at boundary of our domain. Mathematically this is expressed by prescribing 
any loss or inflow of energy at the boundary of our domain. Mathematically this is expressed by prescribing 
126 
the heat flux $\kappa \frac{\partial T}{\partial x}$ to zero. In our simplified one dimensional model this is expressed 
the heat flux $\kappa \frac{\partial T}{\partial x}$ to zero. In our simplified one dimensional model this is expressed 
127 
in the form; 
in the form; 
128 
\begin{equation} 
\begin{equation} 
166 


167 
\subsection{Outline of the Implementation} 
\subsection{Outline of the Implementation} 
168 
\label{sec:outline} 
\label{sec:outline} 
169 
To solve the heat diffusion equation (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 we discuss later in details have 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 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 
170 
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 is as follows: 
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}. 
171 
\begin{enumerate} 
% \begin{enumerate} 
172 
\item create domain 
% \item create domain 
173 
\item create PDE 
% \item create PDE 
174 
\item while end time not reached: 
% \item while end time not reached: 
175 
\begin{enumerate} 
% \begin{enumerate} 
176 
\item set PDE coefficients 
% \item set PDE coefficients 
177 
\item solve PDE 
% \item solve PDE 
178 
\item update time marker 
% \item update time marker 
179 
\end{enumerate} 
% \end{enumerate} 
180 
\item end of calculation 
% \item end of calculation 
181 
\end{enumerate} 
% \end{enumerate} 
182 


183 

\begin{figure}[h!] 
184 

\centering 
185 

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

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

\label{fig:wf} 
188 

\end{figure} 
189 


190 
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 \textbf{objects}. The nice feature of an object is that it is defined by its usage and features 
191 
rather than its actual representation. So we will create a domain object to describe the geometry of the two 
rather than its actual representation. So we will create a domain object to describe the geometry of the two 
192 
granite blocks. Then we define PDEs and spatially distributed values such as the temperature 
granite blocks. Then we define PDEs and spatially distributed values such as the temperature 
196 
\begin{figure}[t] 
\begin{figure}[t] 
197 
\centering 
\centering 
198 
\includegraphics[width=6in]{figures/functionspace.pdf} 
\includegraphics[width=6in]{figures/functionspace.pdf} 

\label{fig:fs} 

199 
\caption{\esc domain construction overview} 
\caption{\esc domain construction overview} 
200 

\label{fig:fs} 
201 
\end{figure} 
\end{figure} 
202 


203 
\subsection{The Domain Constructor in \esc} 
\subsection{The Domain Constructor in \esc} 
204 
\label{ss:domcon} 
\label{ss:domcon} 
205 
It is helpful to have a better understanding how spatially distributed value such as temperature or PDE coefficients are interpreted in \esc. Again 
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. 

from the user's point of view the representation of these spatially distributed values is not relevant. 

206 


207 
There are various ways to construct domain objects. The simplest form is a rectangular shaped region with a length and height. There is 
There are various ways to construct domain objects. The simplest form is a rectangular shaped region with a length and height. There is 
208 
a ready to use function for this named \verb rectangle(). Besides the spatial dimensions this function requires to specify the number of 
a ready to use function for this named \verb rectangle(). Besides the spatial dimensions this function requires to specify the number of 
209 
elements or cells to be used along the length and height, see \reffig{fig:fs}. Any spatially distributed value 
elements or cells to be used along the length and height, see \reffig{fig:fs}. Any spatially distributed value 
210 
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. 
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. 
211 
The quality of the approximation depends  besides other factors mainly on the number of elements being used. In fact, the 
The quality of the approximation depends  besides other factors mainly on the number of elements being used. In fact, the 
212 
approximation becomes better when more elements are used. However, computational costs and compute time grow with the number of 
approximation becomes better when more elements are used. However, computational costs and time grow with the number of 
213 
elements being used. It is therefore important that you find the right balance between the demand in accuracy and acceptable resource usage. 
elements being used. It is therefore important that you find the right balance between the demand in accuracy and acceptable resource usage. 
214 


215 
In general, one can think about a domain object as a composition of nodes and elements. 
In general, one can think about a domain object as a composition of nodes and elements. 
232 


233 
This distinction between different representations of spatially distributed values 
This distinction between different representations of spatially distributed values 
234 
is important in order to be able to vary the degrees of smoothness in a PDE problem. 
is important in order to be able to vary the degrees of smoothness in a PDE problem. 
235 
The coefficients of a PDE do not need be continuous, thus this qualifies as a \verbFunction() type. 
The coefficients of a PDE do not need to be continuous, thus this qualifies as a \verbFunction() type. 
236 
On the other hand a temperature distribution must be continuous and needs to be represented with a \verbContinuousFunction() function space. 
On the other hand a temperature distribution must be continuous and needs to be represented with a \verbContinuousFunction() function space. 
237 
An influx may only be defined at the boundary and is therefore a \verb FunctionOnBoundary() object. 
An influx may only be defined at the boundary and is therefore a \verb FunctionOnBoundary() object. 
238 
\esc allows certain transformations of the function spaces. A \verb ContinuousFunction() can be transformed into a \verbFunctionOnBoundary() 
\esc allows certain transformations of the function spaces. A \verb ContinuousFunction() can be transformed into a \verbFunctionOnBoundary() 
241 


242 
Later in this introduction we discuss how 
Later in this introduction we discuss how 
243 
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 represented by different material coefficients such the 
244 
thermal conductivities $kappa$. A very powerful technique to define these types of PDE 
thermal conductivities $\kappa$. A very powerful technique to define these types of PDE 
245 
coefficients is tagging. Blocks of materials and boundaries can be named and values can be defined on subregions based on their names. 
coefficients is tagging. Blocks of materials and boundaries can be named and values can be defined on subregions based on their names. 
246 
This is simplifying PDE coefficient and flux definitions. It makes scripting much easier. We will discuss this technique in Section~\ref{STEADYSTATE HEAT REFRACTION}. 
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}. 
247 


248 


249 
\subsection{A Clarification for the 1D Case} 
\subsection{A Clarification for the 1D Case} 
286 


287 
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. 
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. 
288 


289 
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. 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. 
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. 
290 
\begin{python} 
\begin{python} 
291 
mx = 500.*m #meters  model length 
mx = 500.*m #meters  model length 
292 
my = 100.*m #meters  model width 
my = 100.*m #meters  model width 
378 
The total energy should stay constant for the example discussed here. 
The total energy should stay constant for the example discussed here. 
379 


380 
\subsection{Running the Script} 
\subsection{Running the Script} 
381 
The script presented so for is available under 
The script presented so far is available under 
382 
\verbexample01a.py. You can edit this file with your favourite text editor. 
\verbexample01a.py. You can edit this file with your favourite text editor. 
383 
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 is not supported under {\it MS Windows} yet.} you can use the \program{runescript} command 
384 
to launch {\it escript} scripts. For the example script use; 
to launch {\it escript} scripts. For the example script use; 
406 
\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 use a variety of free \pyt packages for visualisation. 
407 
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. 
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. 
408 
The \mpl package is a component of SciPy\footnote{\url{http://www.scipy.org}} and is good for basic graphs and plots. 
The \mpl package is a component of SciPy\footnote{\url{http://www.scipy.org}} and is good for basic graphs and plots. 
409 
For more complex visualisation tasks in particular when it comes to two and three dimensional problems we recommend to use more advanced tools for instance \mayavi \footnote{\url{http://code.enthought.com/projects/mayavi/}} 
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/}} 
410 
which bases on the \verbVTK toolkit. We discuss the usage of \verbVTK based 
which is based upon the \verbVTK toolkit. The usage of \verbVTK based 
411 
visualization in Chapter~\ref{Sec:2DHD} where focus on a two dimensional PDE. 
visualization is discussed in Chapter~\ref{Sec:2DHD} which focusses on a two dimensional PDE. 
412 


413 
For our simple problem we have two plotting tasks: we are interested in showing the 
For our simple granite block problem, we have two plotting tasks. Firstly, we are interested in showing the 
414 
behaviour of the total energy over time and how the temperature distribution within the block is 
behaviour of the total energy over time and secondly, how the temperature distribution within the block is 
415 
developing over time. Lets start with the first task. 
developing over time. Lets start with the first task. 
416 


417 
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 energies observed. 
418 
\pyt provides the concept of lists for this. Before 
\pyt provides the concept of lists for this. Before 
419 
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 and the total energies \verbE_list. 
420 
After the new temperature as been calculated by solving the PDE we append the new time marker and total energy 
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 
421 
to the corresponding list using the \verbappend method. With these modifications our script looks as follows: 
to the corresponding list using the \verbappend method. With these modifications our script looks as follows: 
422 
\begin{python} 
\begin{python} 
423 
t_list=[] 
t_list=[] 
431 
t_list.append(t) # add current time mark to record 
t_list.append(t) # add current time mark to record 
432 
E_list.append(totE) # add current total energy to record 
E_list.append(totE) # add current total energy to record 
433 
\end{python} 
\end{python} 
434 
To plot $t$ over $totE$ we use the \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 to be loaded before used; 
435 
\begin{python} 
\begin{python} 
436 
import pylab as pl # plotting package. 
import pylab as pl # plotting package. 
437 
\end{python} 
\end{python} 
457 
\sslist{example01c.py} 
\sslist{example01c.py} 
458 
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 strategy we have used 
459 
for the total energy. Instead of producing a final plot at the end we will generate a 
for the total energy. Instead of producing a final plot at the end we will generate a 
460 
picture at each time step which can be browsed as slide show or composed to a movie. 
picture at each time step which can be browsed as a slide show or composed into a movie. 
461 
The first problem we encounter is that if we produce an image in each time step we need 
The first problem we encounter is that if we produce an image at each time step we need 
462 
to make sure that the images previously generated are not overwritten. 
to make sure that the images previously generated are not overwritten. 
463 


464 
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 convenient to 