14 
We will start by examining a simple one dimensional heat diffusion example. This problem will provide a good launch pad to build our knowledge of \esc and demonstrate how to solve simple partial differential equations (PDEs)\footnote{Wikipedia provides an excellent and comprehensive introduction to \textit{Partial Differential Equations} \url{http://en.wikipedia.org/wiki/Partial_differential_equation}, however their relevance to \esc and implementation should become a clearer as we develop our understanding further into the cookbook.} 
We will start by examining a simple one dimensional heat diffusion example. This problem will provide a good launch pad to build our knowledge of \esc and demonstrate how to solve simple partial differential equations (PDEs)\footnote{Wikipedia provides an excellent and comprehensive introduction to \textit{Partial Differential Equations} \url{http://en.wikipedia.org/wiki/Partial_differential_equation}, however their relevance to \esc and implementation should become a clearer as we develop our understanding further into the cookbook.} 
15 


16 
\section{One Dimensional Heat Diffusion in an Iron Rod} 
\section{One Dimensional Heat Diffusion in an Iron Rod} 
17 
\sslist{onedheatdiff001.py and cblib.py} 

18 
%\label{Sec:1DHDv0} 
%\label{Sec:1DHDv0} 
19 
The first model consists of a simple cold iron bar at a constant temperature of zero \reffig{fig:onedhdmodel}. The bar is perfectly insulated on all sides with a heating element at one end. Intuition tells us that as heat is applied; energy will disperse along the bar via conduction. With time the bar will reach a constant temperature equivalent to that of the heat source. 
The first model consists of a simple cold iron bar at a constant temperature of zero \reffig{fig:onedhdmodel}. The bar is perfectly insulated on all sides with a heating element at one end. Intuition tells us that as heat is applied; energy will disperse along the bar via conduction. With time the bar will reach a constant temperature equivalent to that of the heat source. 
20 
\begin{figure}[h!] 
\begin{figure}[h!] 
119 
% for all $x$ in the domain. 
% for all $x$ in the domain. 
120 


121 
\subsection{Boundary Conditions} 
\subsection{Boundary Conditions} 
122 

\label{SEC BOUNDARY COND} 
123 
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 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. 
124 
A \textbf{Dirichlet boundary condition} is conceptually simpler and is used to prescribe a known value to the unknown  in our example the temperature  on parts of or 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  in our example the temperature  on parts of or the entire boundary of the region of interest. 
125 
For our model problem we want to keep the initial temperature setting on the left side of the 
For our model problem we want to keep the initial temperature setting on the left side of the 
126 
iron bar over time. This defines a Dirichlet boundary condition for the PDE \refEq{eqn:hddisc} to be solved at each time step. 
iron bar over time. This defines a Dirichlet boundary condition for the PDE \refEq{eqn:hddisc} to be solved at each time step. 
127 


128 
On the other end of the iron rod we want to add an appropriate boundary condition to define insolation to prevent 
On the other end of the iron rod we want to add an appropriate boundary condition to define insulation to prevent 
129 
any loss or inflow of energy at the right end of the rod. Mathematically this is expressed by prescribing 
any loss or inflow of energy at the right end of the rod. Mathematically this is expressed by prescribing 
130 
the heat flux $\kappa \frac{\partial T}{\partial x}$ to zero on the right end of the rod 
the heat flux $\kappa \frac{\partial T}{\partial x}$ to zero on the right end of the rod 
131 
In our simplified one dimensional model this is expressed 
In our simplified one dimensional model this is expressed 
155 
It is important to keep in mind that boundary conditions need to be complete and consistent in the sense that 
It is important to keep in mind that boundary conditions need to be complete and consistent in the sense that 
156 
at any point on the boundary either a Dirichlet or a Neuman boundary condition must be set. 
at any point on the boundary either a Dirichlet or a Neuman boundary condition must be set. 
157 


158 
Conviniently, \esc makes default assumption on the boundary conditions which the user may modify where appropriate. 
Conveniently, \esc makes default assumption on the boundary conditions which the user may modify where appropriate. 
159 
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 
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 
160 
$n\hackscore{j}A\hackscore{jl} u\hackscore{,l}=0$.} is; 
$n\hackscore{j}A\hackscore{jl} u\hackscore{,l}=0$.} is; 
161 
\begin{equation}\label{NEUMAN 2} 
\begin{equation}\label{NEUMAN 2} 
174 
\label{sec:outline} 
\label{sec:outline} 
175 
To solve the heat diffusion equation (equation \refEq{eqn:hd}) we will 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 will discuss later in details will have four major steps. Firstly we need to define the domain where we want to 
To solve the heat diffusion equation (equation \refEq{eqn:hd}) we will 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 will discuss later in details will have four major steps. Firstly we need to define the domain where we want to 
176 
calculate the temperature. For our problem this is the iron rod which has a rectangular shape. Secondly we need to define the PDE 
calculate the temperature. For our problem this is the iron rod which has a rectangular shape. Secondly we need to define the PDE 
177 
we need to solve in each time step to get the updated temperature. Thirdly we need to define the 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. 
we need to solve in each time step to get the updated temperature. Thirdly we need to define the 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. As a work flow this takes the form; 
178 

\begin{enumerate} 
179 

\item create domain 
180 

\item create PDE 
181 

\item while end time not reached: 
182 

\begin{enumerate} 
183 

\item set PDE coefficients 
184 

\item solve PDE 
185 

\item update time marker 
186 

\end{enumerate} 
187 

\item end of calculation 
188 

\end{enumerate} 
189 
In the terminology of \pyt the domain and PDE are represented by \textbf{objects}. The nice feature of an object is that it defined by it 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 defined by it usage and features 
190 
rather than its actual representation. So we will create a domain object to decribe the geometry of our iron rod. The main feature 
rather than its actual representation. So we will create a domain object to describe the geometry of our iron rod. The main feature 
191 
of the object we will use is the fact that we can define PDEs and spatially distributed values such as the temperature 
of the object we will use is the fact that we can define PDEs and spatially distributed values such as the temperature 
192 
on a domain. In fact the domain object has many more features  most of them you will 
on a domain. In fact the domain object has many more features  most of them you will 
193 
never use and do not need to understand. Similar a PDE object is defined by the fact that we can define the coefficients of the PDE and solve the PDE. At a 
never use and do not need to understand. Similar a PDE object is defined by the fact that we can define the coefficients of the PDE and solve the PDE. At a 
194 
later stage you may use more advanved features of the PDE class but you need to worry about them only at the point when you use them. 
later stage you may use more advanced features of the PDE class but you need to worry about them only at the point when you use them. 
195 


196 


197 
\begin{figure}[t] 
\begin{figure}[t] 
208 


209 
There are various ways to construct domain objects. The simplest form is as rectangular shaped region with a length and height. There is 
There are various ways to construct domain objects. The simplest form is as rectangular shaped region with a length and height. There is 
210 
a ready to use function call for this. Besides the spatial dimensions the function call will require you to specify the number 
a ready to use function call for this. Besides the spatial dimensions the function call will require you to specify the number 
211 
elements or cells to be used along the length and height, see Figure~\ref{fig:fs}. Any spacially distributed value 
elements or cells to be used along the length and height, see Figure~\ref{fig:fs}. Any specially distributed value 
212 
and the PDE is represented in discrete form using this element representation\footnote{We will use the finite element method (FEM), see \url{http://en.wikipedia.org/wiki/Finite_element_method}i 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 will use the finite element method (FEM), see \url{http://en.wikipedia.org/wiki/Finite_element_method}i for details.}. Therefore we will have access to an approximation of the true PDE solution only. 
213 
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 
214 
approximation becomes better the more elements are used. However, computational costs and compute time grow with the number of 
approximation becomes better the more elements are used. However, computational costs and compute time grow with the number of 
215 
elements being used. It therefore important that you find the right balance between the demand in accuarcy and acceptable resource usage. 
elements being used. It therefore important that you find the right balance between the demand in accuracy and acceptable resource usage. 
216 


217 
In general, one can thinks about a domain object as a composition of nodes and elements. 
In general, one can thinks about a domain object as a composition of nodes and elements. 
218 
As shown in Figure~\ref{fig:fs}, an element is defined by the nodes used to describe its vertices. 
As shown in Figure~\ref{fig:fs}, an element is defined by the nodes used to describe its vertices. 
232 
the \verbFunction(domain).getX() returns the coordinates of numerical integration points within elements, see 
the \verbFunction(domain).getX() returns the coordinates of numerical integration points within elements, see 
233 
Figure~\ref{fig:fs}. 
Figure~\ref{fig:fs}. 
234 


235 
This distiction between different representations of spatial distributed values 
This distinction between different representations of spatial distributed values 
236 
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. 
237 
The coefficients of a PDE need not be continuous thus this qualifies as a \verbFunction() type. 
The coefficients of a PDE need not be continuous thus this qualifies as a \verbFunction() type. 
238 
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. 
239 
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. 
240 
\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 
or \verbFunction(). On the other hand there is not enough information in a \verb FunctionOnBoundary() to transform it to a \verb ContinuousFunction() . 
or \verbFunction(). On the other hand there is not enough information in a \verb FunctionOnBoundary() to transform it to a \verb ContinuousFunction() . 
242 
These transformations, which ar ecalled \textbf{interpolation} are invoked autmatically by \esc if needed. 
These transformations, which are called \textbf{interpolation} are invoked automatically by \esc if needed. 
243 


244 
Later in this introduction we will discuss how 
Later in this introduction we will discuss how 
245 
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 
246 
thermal conductivities $kappa$. A very powerful techique to define these types of PDE 
thermal conductivities $kappa$. A very powerful technique to define these types of PDE 
247 
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. 
248 
This is simplifing PDE coefficient and flux definitions. It makes for much easier scripting. We will discuss this technique in Section~\ref{STEADYSTATE HEAT REFRACTION}. 
This is simplifying PDE coefficient and flux definitions. It makes for much easier scripting. We will discuss this technique in Section~\ref{STEADYSTATE HEAT REFRACTION}. 
249 


250 


251 
\subsection{A Clarification for the 1D Case} 
\subsection{A Clarification for the 1D Case} 
269 


270 
\subsection{Developing a PDE Solution Script} 
\subsection{Developing a PDE Solution Script} 
271 
\label{sec:key} 
\label{sec:key} 
272 

\sslist{onedheatdiffbase.py} 
273 
We will write a simple \pyt script which uses the \modescript, \modfinley and \modmpl modules. 
We will write a simple \pyt script which uses the \modescript, \modfinley and \modmpl modules. 
274 
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.} 
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.} 
275 
that we will require. 
that we will require. 
282 
# A useful unit handling package which will make sure all our units 
# A useful unit handling package which will make sure all our units 
283 
# match up in the equations under SI. 
# match up in the equations under SI. 
284 
from esys.escript.unitsSI import * 
from esys.escript.unitsSI import * 

import pylab as pl #Plotting package. 


import numpy as np #Array package. 


import os #This package is necessary to handle saving our data. 

285 
\end{python} 
\end{python} 
286 
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 model. The module \verb unitsSI provides support for SI unit definitions with our variables; and the \verbos module is needed to handle file outputs once our PDE has been solved. \verb pylab and \verb numpy are modules developed independently of \esc. They are used because they have efficient plotting and array handling capabilities. 
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 model. The module \verb unitsSI provides support for SI unit definitions with our variables. 
287 


288 
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 will need values. Firstly, the model upon which we wish to solve our problem needs to be defined. There are many different types of models in \modescript which we will demonstrate in later tutorials but for our iron rod, we will simply use a rectangular model. 
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 will need values. Firstly, the model upon which we wish to solve our problem needs to be defined. There are many different types of models in \modescript which we will demonstrate in later tutorials but for our iron rod, we will simply use a rectangular model. 
289 


298 
The material constants and the temperature variables must also be defined. For the iron rod in the model they are defined as: 
The material constants and the temperature variables must also be defined. For the iron rod in the model they are defined as: 
299 
\begin{python} 
\begin{python} 
300 
#PDE related 
#PDE related 

q=200. * Celsius #Kelvin  our heat source temperature 


Tref = 0. * Celsius #Kelvin  starting temp of iron bar 

301 
rho = 7874. *kg/m**3 #kg/m^{3} density of iron 
rho = 7874. *kg/m**3 #kg/m^{3} density of iron 
302 
cp = 449.*J/(kg*K) #j/Kg.K thermal capacity 
cp = 449.*J/(kg*K) # J/Kg.K thermal capacity 
303 
rhocp = rho*cp 
rhocp = rho*cp 
304 
kappa = 80.*W/m/K #watts/m.Kthermal conductivity 
kappa = 80.*W/m/K # watts/m.Kthermal conductivity 
305 

qH=0 * J/(sec*m**3) # J/(sec.m^{3}) no heat source 
306 

T0=100 * Celsius # initial temperature at left end of rod 
307 

Tref=20 * Celsius # base temperature 
308 
\end{python} 
\end{python} 
309 
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 where we would like to save the output from the solver. This is simple enough: 
310 
\begin{python} 
\begin{python} 
311 
t=0 #our start time, usually zero 
t=0 * day #our start time, usually zero 
312 
tend=5.*minute #seconds  time to end simulation 
tend=1. * day #  time to end simulation 
313 
outputs = 200 # number of time steps required. 
outputs = 200 # number of time steps required. 
314 
h=(tendt)/outputs #size of time step 
h=(tendt)/outputs #size of time step 
315 
#user warning statement 
#user warning statement 
316 
print "Expected Number of time outputs is: ", (tendt)/h 
print "Expected Number of time outputs is: ", (tendt)/h 
317 
i=0 #loop counter 
i=0 #loop counter 

#the folder to put our outputs in, leave blank "" for script path 


save_path="data/onedheatdiff001" 

318 
\end{python} 
\end{python} 
319 
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 rod as: 
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 rod as: 
320 
\begin{python} 
\begin{python} 
321 
#generate domain using rectangle 
#generate domain using rectangle 
322 
rod = Rectangle(l0=mx,l1=my,n0=ndx, n1=ndy) 
rod = Rectangle(l0=mx,l1=my,n0=ndx, n1=ndy) 
323 
\end{python} 
\end{python} 
324 
\verb rod now describes a domain in the manner of Section \ref{ss:domcon}. As we define our variables, various function spaces will be created to accommodate them. There is an easy way to extract finite points from the domain \verbrod using the domain property function \verbgetX() . This function sets the vertices of each cell as finite points to solve in the solution. If we let \verbx be these finite points, then; 
\verb rod now describes a domain in the manner of Section \ref{ss:domcon}. There is an easy way to extract 
325 

the coordinates of the nodes used to describe the domain \verbrod using the 
326 

domain property function \verbgetX() . This function sets the vertices of each cell as finite points to solve in the solution. If we let \verbx be these finite points, then; 
327 
\begin{python} 
\begin{python} 
328 
#extract finite points  the solution points 
#extract data points  the solution points 
329 
x=rod.getX() 
x=rod.getX() 
330 
\end{python} 
\end{python} 
331 
The data locations of specific function spaces can be returned in a similar manner by extracting the relevant function space from the domain followed by the \verb .getX() operator. 
The data locations of specific function spaces can be returned in a similar manner by extracting the relevant function space from the domain followed by the \verb.getX() method. 
332 


333 
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 will be discussed 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 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 will be discussed later.}. We also need to state the values of our general form variables. 
334 
\begin{python} 
\begin{python} 
335 
mypde=LinearSinglePDE(rod) 
mypde=LinearSinglePDE(rod) 
336 
mypde.setValue(A=kappa*kronecker(rod),D=rhocp/h) 
A=zeros((2,2))) 
337 
\end{python} 
A[0,0]=kappa 
338 

q=whereZero(x[0]) 
339 

mypde.setValue(A=A, D=rhocp/h, q=q, r=T0) 
340 

\end{python} 
341 

The argument \verbq has not been discussed yet: In fact the arguments \verbq and \verbr are used to define 
342 

Dirichlet boundary condition as discussed in Section~\ref{SEC BOUNDARY COND}. In the \esc 
343 

PDE from the argument \verbq indicates by a positive value for which nodes we want to apply a 
344 

Dirichlet boundary condition, ie. where we want to prescribe the value of the PDE solution 
345 

rather then using the PDE. The actually value for the solution to be taken is set by the argument \verbr. 
346 

In our case we want to keep the initial temperature $T0$ on the left face of the rode for all times. Notice, 
347 

that as set to a constant value \verbr is assumed to have the same value 
348 

at all nodes, however only the value at those nodes marked by a positive value by \verbq are actually used. 
349 


350 

In order to set \verbq we use 
351 

\verbwhereZero function. The function returns the value (positive) one for those data points (=nodes) where the argument is equal to zero and otherwise returns (nonpositive) value zero. 
352 

As \verbx[0] given the $x$coordinates of the nodes for the domain, 
353 

\verbwhereZero(x[0]) gives the value $1$ for the nodes at the left end of the rod $x=x_0=0$ and 
354 

zero elsewhere which is exactly what we need. 
355 


356 
In a few special cases it may be possible to decrease the computational time of the solver if our PDE is symmetric. Symmetry of a PDE is defined by; 
In a many cases it may be possible to decrease the computational time of the solver if the PDE is symmetric. 
357 

Symmetry of a PDE is defined by; 
358 
\begin{equation}\label{eqn:symm} 
\begin{equation}\label{eqn:symm} 
359 
A\hackscore{jl}=A\hackscore{lj} 
A\hackscore{jl}=A\hackscore{lj} 
360 
\end{equation} 
\end{equation} 
361 
Symmetry is only dependent on the $A$ coefficient in the general form and the other coefficients $D$ and $d$ as well as the RHS $Y$ and $y$ may take any value. 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 will enable symmetry via; 
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$ may take any value. 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 will enable symmetry via; 
362 
\begin{python} 
\begin{python} 
363 
myPDE.setSymmetryOn() 
myPDE.setSymmetryOn() 
364 
\end{python} 
\end{python} 
365 

Next we need to establish the initial temperature distribution \verbT. We want to have this initial 
366 
We now need to specify our boundary conditions and initial values. The initial values required to solve this PDE are temperatures for each discrete point in our model. We will set our bar to: 
value to be \verbTref except at the left end of the rod $x=0$ where we have the temperature \verbT0. We use; 

\begin{python} 


T = Tref 


\end{python} 


Boundary conditions are a little more difficult. Fortunately the \esc solver will handle our insulated boundary conditions by default with a zero flux operator. However, we will need to apply our heat source $q_{H}$ to the end of the bar at $x=0$ . \esc makes this easy by letting us define areas in our model. The finite points in the model were previously defined as \verb x and it is possible to set all of points that satisfy $x=0$ to \verb q via the \verb whereZero() function. There are a few \verb where functions available in \esc. They will return a value \verb 1 where they are satisfied and \verb 0 where they are not. In this case our \verb qH is only applied to the far LHS of our model as required. 

367 
\begin{python} 
\begin{python} 
368 
# ... set heat source: .... 
# ... set initial temperature .... 
369 
qH=q*whereZero(x[0]) 
T = T0*whereZero(x[0])+Tref*(1whereZero(x[0])) 
370 
\end{python} 
\end{python} 
371 

Finally we will initialize 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 will initialise an iteration loop to solve our PDE for all the time steps we specified in the variable section. As the RHS 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. 

372 
\begin{python} 
\begin{python} 
373 
while t<=tend: 
while t < tend: 
374 
i+=1 #increment the counter 
i+=1 #increment the counter 
375 
t+=h #increment the current time 
t+=h #increment the current time 
376 
mypde.setValue(Y=qH+rhocp/h*T) #set variable PDE coefficients 
mypde.setValue(Y=qH+rhocp/h*T) # set variable PDE coefficients 
377 
T=mypde.getSolution() #get the PDE solution 
T=mypde.getSolution() #get the PDE solution 
378 
totT = rhocp*T #get the total heat solution in the system 
totE = integrate(rhocp*T) #get the total heat (energy) in the system 
379 
\end{python} 
\end{python} 
380 

The last statement in this script calculates the total energy in the system as volume integral 
381 

of $\rho \c_p T$ over the rod. 
382 


383 

\subsection{Plotting the Total Energy} 
384 

\sslist{onedheatdiff001.py} 
385 


386 
\subsection{Plotting the heat solutions} 
\esc does not include its own plotting capabilities. However, it is possible to use a variety of free \pyt packages for visualization. 
387 
Visualisation of the solution can be achieved using \mpl a module contained within \pylab. We start by modifying our solution script from before. Prior to the \verb while loop we will need to extract our finite solution points to a data object that is compatible with \mpl. First it is necessary to convert \verb x to a list of tuples. These are then converted to a \numpy array and the $x$ locations extracted via an array slice to the variable \verb plx . 
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. 
388 

The \mpl package is a component of SciPy\footnote{\url{http://www.scipy.org}} and is good for basic graphs and plots. 
389 

For more complex visualisation tasks in particular when it comes to two and three dimensional problems it is recommended to us more advanced tools for instance \mayavi \footnote{\url{http://code.enthought.com/projects/mayavi/}} 
390 

which bases the \verbVTK toolkit. We will discuss the usage of \verbVTK based 
391 

visualization in Chapter~\ref{Sec:2DHD} where will discuss a two dimensional PDE. 
392 


393 

For our simple problem we have two plotting tasks: Firstly we are interested in showing the 
394 

behavior of the total energy over time and secondly in how the temperature distribution within the rod is 
395 

developing over time. Lets start with the first task. 
396 


397 

\begin{figure} 
398 

\begin{center} 
399 

\includegraphics[width=4in]{figures/ttrodpyplot150} 
400 

\caption{Total Energy in Rod over Time (in seconds).} 
401 

\label{fig:onedheatout1} 
402 

\end{center} 
403 

\end{figure} 
404 


405 

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

\pyt provides the concept of lists for this. Before 
407 

the time loop is opened we create empty lists for the time marks \verbt_list and the total energies \verbE_list. 
408 

After the new temperature as been calculated by solving the PDE we append the new time marker and total energy 
409 

to the corresponding list using the \verbappend method. With these modifications the script looks as follows: 
410 

\begin{python} 
411 

t_list=[] 
412 

E_list=[] 
413 

# ... start iteration: 
414 

while t<tend: 
415 

t+=h 
416 

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

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

totE=integrate(rhocp*T) 
419 

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

E_list.append(totE) # add current total energy to record 
421 

\end{python} 
422 

To plot $t$ over $totE$ we use the \mpl a module contained within \pylab which needs to be loaded before used; 
423 

\begin{python} 
424 

import pylab as pl # plotting package. 
425 

\end{python} 
426 

Here we are not using the \verbfrom pylab import * in order to name clashes for function names 
427 

with \esc. 
428 


429 

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

\begin{python} 
431 

pl.plot(t_list,E_list) 
432 

pl.title("Total Energy") 
433 

pl.savefig("totE.png") 
434 

\end{python} 
435 

The first statement hands over the time marks and corresponding total energies to the plotter. 
436 

The second statement is setting the title for the plot. The last statement renders the plot and writes the 
437 

result into the file \verbtotE.png which can be displayed by (almost) any image viewer. Your result should look 
438 

similar to Figure~\ref{fig:onedheatout1}. 
439 


440 

\subsection{Plotting the Temperature Distribution} 
441 

\sslist{onedheatdiff001b.py} 
442 

For plotting the spatial distribution of the temperature we need to modify the strategy we have used 
443 

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

picture at each time step which can be browsed as slide show or composed to a movie. 
445 

The first problem we encounter is that if we produce an image in each time step we need 
446 

to make sure that the images previously generated are not overwritten. 
447 


448 

To develop an incrementing file name we can use the following convention. It is convenient to 
449 

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

into a separate directory. As part of the \verbos module\footnote{The \texttt{os} module provides 
451 

a powerful interface to interact with the operating system, see \url{http://docs.python.org/library/os.html}.} \pyt 
452 

provides the \verbos.path.join command to build file and 
453 

directory names in a platform independent way. Assuming that 
454 

\verbsave_path is name of directory we want to put the results the command is; 
455 

\begin{python} 
456 

import os 
457 

os.path.join(save_path, "tempT%03d.png"%i ) 
458 

\end{python} 
459 

where \verbi is the time step counter. 
460 

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 in, for example a single subfolder called \verb data would be defined by; 
461 

\begin{verbatim} 
462 

save_path = "data" 
463 

\end{verbatim} 
464 

while a subfolder of \verb data called \verb onedheatdiff001 would be defined by; 
465 

\begin{verbatim} 
466 

save_path = os.path.join("data","onedheatdiff001") 
467 

\end{verbatim} 
468 

The second argument of \verb join \xspace contains a string which is the filename or subdirectory name. We can use the operator \verb% to increment our file names with the value \verbi denoting a incrementing counter. The substring \verb %03d does this by defining the following parameters; 
469 

\begin{itemize} 
470 

\item \verb 0 becomes the padding number; 
471 

\item \verb 3 tells us the amount of padding numbers that are required; and 
472 

\item \verb d indicates the end of the \verb % operator. 
473 

\end{itemize} 
474 

To increment the file name a \verb %i is required directly after the operation the string is involved in. When correctly implemented the output files from this command would be place in the directory defined by \verb save_path as; 
475 

\begin{verbatim} 
476 

rodpyplot.png 
477 

rodpyplot.png 
478 

rodpyplot.png 
479 

... 
480 

\end{verbatim} 
481 

and so on. 
482 


483 

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

\begin{verbatim} 
485 

mkDir(save_path) 
486 

\end{verbatim} 
487 

will check for the existence of \verb save_path and if missing, make the required directories. 
488 


489 

We start by modifying our solution script from before. 
490 

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 data points used to represent 
491 

the temperature as a \pyt list of tuples. As a solution of a PDE 
492 

the temperature has the \verbSolution(rod) function space attribute. We use 
493 

the \verbgetX() method to get the coordinates of the data points as an \esc object 
494 

which is then converted to a \numpy array. The $x$ component is then extracted via an array slice to the variable \verbplx; 
495 
\begin{python} 
\begin{python} 
496 

import numpy as np # array package. 
497 
#convert solution points for plotting 
#convert solution points for plotting 
498 
plx = x.toListOfTuples() 
plx = Solution(rod).getX().toListOfTuples() 
499 
plx = np.array(plx) #convert to tuple to numpy array 
plx = np.array(plx) # convert to tuple to numpy array 
500 
plx = plx[:,0] #extract x locations 
plx = plx[:,0] # extract x locations 

\end{python} 


As there are two solution outputs, we will generate two plots and save each to a file for every time step in the solution. The following is appended to the end of the \verb while loop and creates two figures. The first figure is for the temperature distribution, and the second the total temperature in the bar. Both cases are similar with a few minor changes for scale and labelling. We start by converting the solution to a tuple and then plotting this against our \textit{x coordinates} \verb plx from before. The axis is then standardised and a title applied. Finally, the figure is saved to a *.png file and cleared for the following iteration. 


\begin{python} 


#establish figure 1 for temperature vs x plots 


tempT = T.toListOfTuples(scalarastuple=False) 


pl.figure(1) #current figure 


pl.plot(plx,tempT) #plot solution 


#define axis extents and title 


pl.axis([0,1.0,273.14990+0.00008,0.004+273.1499]) 


pl.title("Temperature accross Rod") 


#save figure to file 


pl.savefig(os.path.join(save_path+"/tempT","rodpyplot%03d.png") %i) 


pl.clf() #clear figure 





#establish figure 2 for total temperature vs x plots and repeat 


tottempT = totT.toListOfTuples(scalarastuple=False) 


pl.figure(2) 


pl.plot(plx,tottempT) 


pl.axis([0,1.0,9.657E08,12000+9.657E08]) 


pl.title("Total temperature across Rod") 


pl.savefig(os.path.join(save_path+"/totT","ttrodpyplot%03d.png")%i) 


pl.clf() 

501 
\end{python} 
\end{python} 
502 


503 
\begin{figure} 
\begin{figure} 
504 
\begin{center} 
\begin{center} 
505 
\includegraphics[width=4in]{figures/ttrodpyplot150} 
\includegraphics[width=4in]{figures/rodpyplot001} 
506 
\caption{Total temperature ($T$) distribution in rod at $t=150$} 
\includegraphics[width=4in]{figures/rodpyplot050} 
507 

\includegraphics[width=4in]{figures/rodpyplot200} 
508 

\caption{Temperature ($T$) distribution in rod at time steps $1$, $50$ and $200$.} 
509 
\label{fig:onedheatout} 
\label{fig:onedheatout} 
510 
\end{center} 
\end{center} 
511 
\end{figure} 
\end{figure} 
512 


513 
\subsubsection{Parallel scripts (MPI)} 
For each time step we will generate a plot of the temperature distribution and save each to a file. We use the same 
514 
In some of the example files for this cookbook the plotting commands are a little different. 
techniques provided by \mpl as we have used to plot the total energy over time. 
515 
For example, 
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. 
516 
\begin{python} 
Finally, the figure is saved to a \verb*.png file and cleared for the following iteration. 
517 
if getMPIRankWorld() == 0: 
\begin{python} 
518 
pl.savefig(os.path.join(save_path+"/totT","ttrodpyplot%03d.png")%i) 
# ... start iteration: 
519 
pl.clf() 
while t<tend: 
520 

.... 
521 

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

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

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

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

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

# add title 
527 

pl.title("Temperature across rod at time %e minutes"%(t/minutes)) 
528 

#save figure to file 
529 

pl.savefig(os.path.join(save_path,"tempT","rodpyplot%03d.png") %i) 
530 
\end{python} 
\end{python} 
531 

Some results are shown in Figure~\ref{fig:onedheatout}. 

The additional \verb if statement is not necessary for normal desktop use. 


It becomes important for scripts run on parallel computers. 


Its purpose is to ensure that only one copy of the file is written. 


For more details on writing scripts for parallel computing please consult the \emph{user's guide}. 

532 


533 
\subsection{Make a video} 
\subsection{Make a video} 
534 
Our saved plots from the previous section can be cast into a video using the following command appended to the end of the script. \verb mencoder is linux only however, and other platform users will need to use an alternative video encoder. 
Our saved plots from the previous section can be cast into a video using the following command appended to the end of the script. \verb mencoder is Linux only however, and other platform users will need to use an alternative video encoder. 
535 
\begin{python} 
\begin{python} 
536 
# compile the *.png files to create two *.avi videos that show T change 
# compile the *.png files to create a *.avi videos that show T change 
537 
# with time. This operation uses linux mencoder. For other operating 
# with time. This operation uses Linux mencoder. For other operating 
538 
# systems it is possible to use your favourite video compiler to 
# systems it is possible to use your favorite video compiler to 
539 
# convert image files to videos. 
# convert image files to videos. 
540 


541 
os.system("mencoder mf://"+save_path+"/tempT"+"/*.png mf type=png:\ 
os.system("mencoder mf://"+save_path+"/tempT"+"/*.png mf type=png:\ 
542 
w=800:h=600:fps=25 ovc lavc lavcopts vcodec=mpeg4 oac copy o \ 
w=800:h=600:fps=25 ovc lavc lavcopts vcodec=mpeg4 oac copy o \ 
543 
onedheatdiff001tempT.avi") 
onedheatdiff001tempT.avi") 




os.system("mencoder mf://"+save_path+"/totT"+"/*.png mf type=png:\ 


w=800:h=600:fps=25 ovc lavc lavcopts vcodec=mpeg4 oac copy o \ 


onedheatdiff001totT.avi") 

544 
\end{python} 
\end{python} 
545 

