1 
ahallam 
2401 

2 


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


% 
4 
jfenwick 
2881 
% Copyright (c) 20032010 by University of Queensland 
5 
ahallam 
2401 
% Earth Systems Science Computational Center (ESSCC) 
6 


% http://www.uq.edu.au/esscc 
7 


% 
8 


% Primary Business: Queensland, Australia 
9 


% Licensed under the Open Software License version 3.0 
10 


% http://www.opensource.org/licenses/osl3.0.php 
11 


% 
12 


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
13 



14 
gross 
2949 
\section{Example 1: One Dimensional Heat Diffusion in Granite} 
15 
gross 
2905 
\label{Sec:1DHDv00} 
16 
gross 
2878 

17 
ahallam 
2979 
The first model consists of two blocks of isotropic material, for instance 
18 
ahallam 
3370 
granite, sitting next to each other (\autoref{fig:onedgbmodel}). 
19 
ahallam 
2979 
Initial temperature in \textit{Block 1} is \verbT1 and in \textit{Block 2} is 
20 


\verbT2. 
21 
gross 
2905 
We assume that the system is insulated. 
22 


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


cooler one until both 
25 
gross 
2905 
blocks have the same temperature. 
26 



27 
ahallam 
3370 
\begin{figure}[ht] 
28 


\centerline{\includegraphics[width=4.in]{figures/onedheatdiff001}} 
29 


\caption{Example 1: Temperature differential along a single interface between 
30 


two granite blocks.} 
31 


\label{fig:onedgbmodel} 
32 


\end{figure} 
33 



34 
ahallam 
2495 
\subsection{1D Heat Diffusion Equation} 
35 
caltinay 
2982 
We can model the heat distribution of this problem over time using the one 
36 
ahallam 
2979 
dimensional heat diffusion equation\footnote{A detailed discussion on how the 
37 


heat diffusion equation is derived can be found at 
38 


\url{ 
39 


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


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


\end{equation} 
46 
jfenwick 
3308 
where $\rho$ is the material density, $c_p$ is the specific heat and 
47 
ahallam 
2979 
$\kappa$ is the thermal 
48 


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


from Wikipedia 
50 


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


assume that these material 
52 
gross 
2861 
parameters are \textbf{constant}. 
53 
ahallam 
2979 
The heat source is defined by the right hand side of \refEq{eqn:hd} as 
54 
jfenwick 
3308 
$q_{H}$; this can take the form of a constant or a function of time and 
55 


space. For example $q_{H} = q_{0}e^{\gamma t}$ where we have 
56 
ahallam 
2979 
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 
caltinay 
2982 
and our signed distance from the blockblock interface $x$. 
62 
ahallam 
2401 

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


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


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


do this, a numerical approach is required to discretise 
68 


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


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


These parameters together define the model. While discretisation introduces 
71 


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


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

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


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


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

78 
caltinay 
2982 
For time discretisation we use the Backward Euler approximation 
79 
ahallam 
2979 
scheme\footnote{see \url{http://en.wikipedia.org/wiki/Euler_method}}. It is 
80 
caltinay 
2982 
based on the approximation 
81 
gross 
2861 
\begin{equation} 
82 


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


\label{eqn:beuler} 
84 


\end{equation} 
85 


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


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


\begin{equation} 
88 


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


\label{eqn:Tbeuler} 
90 


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


has 
93 
gross 
2861 
\begin{equation} 
94 


\begin{array}{rcl} 
95 


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


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


\end{array} 
98 


\label{eqn:Neuler} 
99 


\end{equation} 
100 


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


\begin{equation} 
102 
jfenwick 
3308 
\frac{\rho c_p}{h} (T^{(n)}  T^{(n1)})  \kappa \frac{\partial^{2} 
103 


T^{(n)}}{\partial x^{2}} = q_H 
104 
gross 
2861 
\label{eqn:hddisc} 
105 


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


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


one can evaluate the spatial derivative term at the previous time $t^{(n1)}$. 
109 
caltinay 
2982 
This approach is called the \textbf{forward Euler} scheme. This scheme can 
110 


provide some computational advantages, which 
111 
ahallam 
2979 
are not discussed here. However, the \textbf{forward Euler} scheme has a major 
112 


disadvantage. Namely, depending on the 
113 


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


derivative term, the time step size $h$ needs to be chosen sufficiently small to 
115 
caltinay 
2982 
achieve a stable temperature when progressing in time. Stability is achieved if 
116 


the temperature does not grow beyond its initial bounds and becomes 
117 
ahallam 
2979 
nonphysical. 
118 


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


that under the assumption of a 
120 


physically correct problem setup the temperature approximation remains physical 
121 


for all time steps. 
122 
caltinay 
2982 
The user needs to keep in mind that the discretisation error introduced by 
123 
ahallam 
2979 
\refEq{eqn:beuler} 
124 


is sufficiently small, thus a good approximation of the true temperature is 
125 
caltinay 
2982 
computed. It is therefore very important that any results are viewed with 
126 


caution. For example, one may compare the results for different time and 
127 


spatial step sizes. 
128 
gross 
2861 

129 


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

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


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


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


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


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


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


operator} represents 
146 


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


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


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


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


\refEq{eqn:hddisc} 
157 
gross 
2861 
we rearrange \refEq{eqn:hddisc}; 
158 
ahallam 
2411 
\begin{equation} 
159 
jfenwick 
3308 
\frac{\rho c_p}{h} T^{(n)}  \kappa \frac{\partial^2 T^{(n)}}{\partial 
160 


x^2} = q_H + \frac{\rho c_p}{h} T^{(n1)} 
161 
ahallam 
2494 
\label{eqn:hdgenf} 
162 


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


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


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


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


to be constant. 
168 


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


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


that; 
171 
gross 
2862 
\begin{equation}\label{ESCRIPT SET} 
172 
gross 
2861 
u=T^{(n)}; 
173 
jfenwick 
3308 
A = \kappa; D = \frac{\rho c _{p}}{h}; f = q _{H} + \frac{\rho 
174 


c_p}{h} T^{(n1)} 
175 
ahallam 
2494 
\end{equation} 
176 



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


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


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


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


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


respectively. 
185 


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


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


on parts of the boundary or on the entire boundary of the region of interest. 
188 
caltinay 
2982 
We discuss the Dirichlet boundary condition in our second example presented in 
189 
ahallam 
2979 
Section~\ref{Sec:1DHDv0}. 
190 
ahallam 
2495 

191 
ahallam 
2979 
However, for this example we have made the model assumption that the system is 
192 
caltinay 
2982 
insulated, so we need to add an appropriate boundary condition to prevent 
193 
ahallam 
2979 
any loss or inflow of energy at the boundary of our domain. Mathematically this 
194 


is expressed by prescribing 
195 


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


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


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


\kappa \nabla T \cdot n = 0 
204 


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


of the domain. 
207 
caltinay 
2982 
The $\cdot$ (dot) refers to the dot product of the vectors $\nabla T$ and $n$. 
208 
ahallam 
2979 
In fact, the term $\nabla T \cdot n$ is the normal derivative of 
209 


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


for the normal 
211 
jfenwick 
3308 
derivative is $T_{,i} n_i$.}; 
212 
ahallam 
2645 
\begin{equation} 
213 
gross 
2862 
\nabla T \cdot n = \frac{\partial T}{\partial n} \; . 
214 
ahallam 
2645 
\end{equation} 
215 
caltinay 
2982 
A condition of the type \refEq{NEUMAN 1} defines a \textbf{Neumann boundary 
216 
ahallam 
2979 
condition} for the PDE. 
217 
ahallam 
2494 

218 
gross 
2905 
The PDE \refEq{eqn:hdgenf} 
219 
caltinay 
2982 
and the Neumann boundary condition~\ref{eqn:hdgenf} (potentially together with 
220 
ahallam 
2979 
the Dirichlet boundary conditions) define a \textbf{boundary value problem}. 
221 


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


the solution in the 
223 


interior of the domain from information known on the boundary only. In most 
224 
caltinay 
2982 
cases we use the term partial differential equation but in fact it is a 
225 


boundary value problem. 
226 
ahallam 
2979 
It is important to keep in mind that boundary conditions need to be complete and 
227 


consistent in the sense that 
228 
caltinay 
2982 
at any point on the boundary either a Dirichlet or a Neumann boundary condition 
229 
ahallam 
2979 
must be set. 
230 
gross 
2862 

231 
caltinay 
2982 
Conveniently, \esc makes a default assumption on the boundary conditions which 
232 


the user may modify where appropriate. 
233 
ahallam 
2979 
For a problem of the form in~\refEq{eqn:commonform nabla} the default 
234 
caltinay 
2982 
condition\footnote{In the \esc user guide which uses the Einstein convention 
235 


this is written as 
236 
jfenwick 
3308 
$n_{j}A_{jl} u_{,l}=0$.} is; 
237 
gross 
2862 
\begin{equation}\label{NEUMAN 2} 
238 
gross 
2948 
n\cdot A \cdot\nabla u = 0 
239 
gross 
2862 
\end{equation} 
240 
ahallam 
2979 
which is used everywhere on the boundary. Again $n$ denotes the outer normal 
241 


field. 
242 


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


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


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


\refEq{ESCRIPT SET} this 
246 
gross 
2862 
condition translates into 
247 
gross 
2867 
\begin{equation}\label{NEUMAN 2b} 
248 
gross 
2862 
\kappa \frac{\partial T}{\partial x} = 0 
249 


\end{equation} 
250 
caltinay 
2982 
for the boundary of the domain. This is identical to the Neumann boundary 
251 
ahallam 
2979 
condition we want to set. \esc will take care of this condition for us. We 
252 


discuss the Dirichlet boundary condition later. 
253 
gross 
2862 

254 
gross 
2870 
\subsection{Outline of the Implementation} 
255 


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


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


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


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


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


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


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


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


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


need to be repeated until the final time marker has been reached. The work flow 
266 
caltinay 
2982 
is described in \reffig{fig:wf}. 
267 
ahallam 
2975 
% \begin{enumerate} 
268 


% \item create domain 
269 


% \item create PDE 
270 


% \item while end time not reached: 
271 


% \begin{enumerate} 
272 


% \item set PDE coefficients 
273 


% \item solve PDE 
274 


% \item update time marker 
275 


% \end{enumerate} 
276 


% \item end of calculation 
277 


% \end{enumerate} 
278 



279 


\begin{figure}[h!] 
280 


\centering 
281 


\includegraphics[width=1in]{figures/workflow.png} 
282 
caltinay 
2982 
\caption{Workflow for developing an \esc model and solution} 
283 
ahallam 
2975 
\label{fig:wf} 
284 


\end{figure} 
285 



286 
ahallam 
2979 
In the terminology of \pyt, the domain and PDE are represented by 
287 


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


usage and features 
289 


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


describe the geometry of the two 
291 


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


temperature 
293 


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


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


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

297 



298 
ahallam 
3029 
\begin{figure}[htp] 
299 
gross 
2870 
\centering 
300 


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


\end{figure} 
304 



305 


\subsection{The Domain Constructor in \esc} 
306 


\label{ss:domcon} 
307 
ahallam 
2979 
Whilst it is not strictly relevant or necessary, a better understanding of 
308 
caltinay 
2982 
how values are spatially distributed (\textit{e.g.} Temperature) and how PDE 
309 
ahallam 
2979 
coefficients are interpreted in \esc can be helpful. 
310 
gross 
2870 

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


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


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


dimensions this function requires to specify the number of 
315 


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


Any spatially distributed value 
317 


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


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


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


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


only. 
322 
caltinay 
2982 
The quality of the approximation depends  besides other factors  mainly on the 
323 
ahallam 
2979 
number of elements being used. In fact, the 
324 


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


cost grows with the number of 
326 


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


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

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


elements. 
331 


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


describe its vertices. 
333 
caltinay 
2982 
To represent spatially distributed values the user can use 
334 
ahallam 
2979 
the values at the nodes, at the elements in the interior of the domain or at the 
335 
caltinay 
2982 
elements located on the surface of the domain. 
336 
ahallam 
2979 
The different approach used to represent values is called \textbf{function 
337 


space} and is attached to all objects 
338 
caltinay 
2982 
in \esc representing a spatially distributed value such as the solution of 
339 


a PDE. The three function spaces we use at the moment are; 
340 
gross 
2870 
\begin{enumerate} 
341 


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


\item the elements/cells, called by \verbFunction(domain) ; and 
343 
caltinay 
2982 
\item the boundary, called by \verbFunctionOnBoundary(domain). 
344 
gross 
2870 
\end{enumerate} 
345 
ahallam 
2979 
A function space object such as \verbContinuousFunction(domain) has the method 
346 


\verbgetX attached to it. This method returns the 
347 


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


particular function space. So the 
349 


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


nodes used to describe the domain while 
351 
caltinay 
2982 
\verbFunction(domain).getX() returns the coordinates of numerical 
352 


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

354 
ahallam 
2979 
This distinction between different representations of spatially distributed 
355 


values 
356 


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


problem. 
358 


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


\verbFunction() type. 
360 


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


represented with a \verbContinuousFunction() function space. 
362 
caltinay 
2982 
An influx may only be defined at the boundary and is therefore a 
363 


\verbFunctionOnBoundary() object. 
364 


\esc allows certain transformations of the function spaces. A 
365 


\verbContinuousFunction() can be transformed into a 
366 


\verbFunctionOnBoundary() or \verbFunction(). On the other hand there is 
367 


not enough information in a \verbFunctionOnBoundary() to transform it to a 
368 


\verbContinuousFunction(). 
369 
ahallam 
2979 
These transformations, which are called \textbf{interpolation} are invoked 
370 


automatically by \esc if needed. 
371 
gross 
2870 

372 
artak 
2963 
Later in this introduction we discuss how 
373 
ahallam 
2979 
to define specific areas of geometry with different materials which are 
374 
caltinay 
2982 
represented by different material coefficients such as the 
375 
ahallam 
2979 
thermal conductivities $\kappa$. A very powerful technique to define these types 
376 


of PDE 
377 


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


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


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


scripting much easier and we will discuss this technique in 
381 


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

383 



384 


\subsection{A Clarification for the 1D Case} 
385 
gross 
2931 
\label{SEC: 1D CLARIFICATION} 
386 
ahallam 
2979 
It is necessary for clarification that we revisit our general PDE from 
387 
caltinay 
2982 
\refeq{eqn:commonform nabla} for a two dimensional domain. \esc is inherently 
388 


designed to solve problems that are multidimensional and so 
389 
ahallam 
2979 
\refEq{eqn:commonform nabla} needs to be read as a higher dimensional problem. 
390 


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


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


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


takes the following form; 
394 
gross 
2477 
\begin{equation}\label{eqn:commonform2D} 
395 
jfenwick 
3308 
A_{00}\frac{\partial^{2}u}{\partial x^{2}} 
396 


A_{01}\frac{\partial^{2}u}{\partial x\partial y} 
397 


A_{10}\frac{\partial^{2}u}{\partial y\partial x} 
398 


A_{11}\frac{\partial^{2}u}{\partial y^{2}} 
399 
gross 
2477 
+ Du = f 
400 


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


a compact formulation which is also independent from the spatial dimension. 
404 
artak 
2963 
To make the general PDE \refEq{eqn:commonform2D} one dimensional as 
405 
gross 
2861 
shown in \refEq{eqn:commonform} we need to set 
406 
ahallam 
2606 
\begin{equation} 
407 
jfenwick 
3308 
A_{00}=A; A_{01}=A_{10}=A_{11}=0 
408 
gross 
2477 
\end{equation} 
409 



410 
gross 
2867 

411 
ahallam 
2495 
\subsection{Developing a PDE Solution Script} 
412 
ahallam 
2801 
\label{sec:key} 
413 
gross 
2949 
\sslist{example01a.py} 
414 
ahallam 
2979 
We write a simple \pyt script which uses the \modescript, \modfinley and \modmpl 
415 


modules. 
416 


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


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


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


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


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


more complicated like those from our \esc library.} 
422 
ahallam 
2495 
that we will require. 
423 
ahallam 
2775 
\begin{python} 
424 
ahallam 
2495 
from esys.escript import * 
425 
ahallam 
2606 
# This defines the LinearPDE module as LinearPDE 
426 


from esys.escript.linearPDEs import LinearPDE 
427 


# This imports the rectangle domain function from finley. 
428 


from esys.finley import Rectangle 
429 


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


# match up in the equations under SI. 
431 


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


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


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


ease of use later in the script. \verbRectangle is going to be our type of 
437 
caltinay 
2982 
domain. The module \verbunitsSI provides support for SI unit definitions with 
438 
ahallam 
2979 
our variables. 
439 
gross 
2477 

440 
ahallam 
2979 
Once our library dependencies have been established, defining the problem 
441 


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


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


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


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


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


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


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


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


rectangular domain. 
450 
ahallam 
2401 

451 
ahallam 
2979 
Using a rectangular domain simplifies our granite blocks (which would in reality 
452 


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


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


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


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


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


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


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


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


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


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


problem, the bar is defined as being 1 metre long. An appropriate step size 
463 
caltinay 
2982 
\verbndx would be 1 to 10\% of the length. Our \verbndy needs only be 1, 
464 


this is because our problem stipulates no partial derivatives in the $y$ 
465 


direction. 
466 
ahallam 
2979 
Thus the temperature does not vary with $y$. Hence, the model parameters can be 
467 
caltinay 
2982 
defined as follows; note we have used the \verbunitsSI convention to make sure 
468 
ahallam 
2979 
all our input units are converted to SI. 
469 
ahallam 
2775 
\begin{python} 
470 
gross 
2905 
mx = 500.*m #meters  model length 
471 


my = 100.*m #meters  model width 
472 


ndx = 50 # mesh steps in x direction 
473 


ndy = 1 # mesh steps in y direction 
474 


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


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


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


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


where we would like to save the output from the solver. This is simple enough: 
490 
ahallam 
2775 
\begin{python} 
491 
ahallam 
3373 
t=0 * day # our start time, usually zero 
492 


tend=50 * yr #  time to end simulation 
493 
ahallam 
2495 
outputs = 200 # number of time steps required. 
494 


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


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


i=0 #loop counter 
498 
ahallam 
2775 
\end{python} 
499 
ahallam 
2979 
Now that we know our inputs we will build a domain using the 
500 
caltinay 
2982 
\verbRectangle() function from \FINLEY. The four arguments allow us to 
501 


define our domain \verbmodel as: 
502 
ahallam 
2775 
\begin{python} 
503 
ahallam 
2606 
#generate domain using rectangle 
504 
gross 
2905 
blocks = Rectangle(l0=mx,l1=my,n0=ndx, n1=ndy) 
505 
ahallam 
2775 
\end{python} 
506 
caltinay 
2982 
\verbblocks now describes a domain in the manner of Section \ref{ss:domcon}. 
507 
ahallam 
2658 

508 
caltinay 
2982 
With a domain and all the required variables established, it is now possible to 
509 
ahallam 
2979 
set up our PDE so that it can be solved by \esc. The first step is to define the 
510 


type of PDE that we are trying to solve in each time step. In this example it is 
511 
caltinay 
2982 
a single linear PDE\footnote{in contrast to a system of PDEs which we discuss 
512 
ahallam 
2979 
later.}. We also need to state the values of our general form variables. 
513 
ahallam 
2775 
\begin{python} 
514 
gross 
2905 
mypde=LinearPDE(blocks) 
515 
gross 
2878 
A=zeros((2,2))) 
516 


A[0,0]=kappa 
517 
gross 
2905 
mypde.setValue(A=A, D=rhocp/h) 
518 
ahallam 
2775 
\end{python} 
519 
caltinay 
2982 
In many cases it may be possible to decrease the computational time of the 
520 
ahallam 
2979 
solver if the PDE is symmetric. 
521 
gross 
2878 
Symmetry of a PDE is defined by; 
522 
ahallam 
2495 
\begin{equation}\label{eqn:symm} 
523 
jfenwick 
3308 
A_{jl}=A_{lj} 
524 
ahallam 
2495 
\end{equation} 
525 
ahallam 
2979 
Symmetry is only dependent on the $A$ coefficient in the general form and the 
526 


other coefficients $D$ as well as the right hand side $Y$. From the above 
527 
caltinay 
2982 
definition we can see that our PDE is symmetric. The \verbLinearPDE class 
528 
ahallam 
2979 
provides the method \method{checkSymmetry} to check if the given PDE is 
529 


symmetric. As our PDE is symmetrical we enable symmetry via; 
530 
ahallam 
2775 
\begin{python} 
531 
caltinay 
2982 
myPDE.setSymmetryOn() 
532 
ahallam 
2775 
\end{python} 
533 
ahallam 
2979 
Next we need to establish the initial temperature distribution \verbT. We need 
534 


to 
535 


assign the value \verbT1 to all sample points left to the contact interface at 
536 
jfenwick 
3308 
$x_{0}=\frac{mx}{2}$ 
537 
gross 
2905 
and the value \verbT2 right to the contact interface. \esc 
538 
caltinay 
2982 
provides the \verbwhereNegative function to construct this. More 
539 


specifically, \verbwhereNegative returns the value $1$ at those sample points 
540 


where the argument has a negative value. Otherwise zero is returned. 
541 
jfenwick 
3308 
If \verbx are the $x_{0}$ 
542 
gross 
2905 
coordinates of the sample points used to represent the temperature distribution 
543 


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


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


the right of the interface. So with; 
546 
gross 
2878 
\begin{python} 
547 


# ... set initial temperature .... 
548 
gross 
2905 
T= T1*whereNegative(x[0]boundloc)+T2*(1whereNegative(x[0]boundloc)) 
549 
gross 
2878 
\end{python} 
550 
ahallam 
2979 
we get the desired temperature distribution. To get the actual sample points 
551 
caltinay 
2982 
\verbx we use the \verbgetX() method of the function space 
552 


\verbSolution(blocks) which is used to represent the solution of a PDE; 
553 
gross 
2878 
\begin{python} 
554 
gross 
2905 
x=Solution(blocks).getX() 
555 


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


\verbSolution(blocks) 
558 


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


representation. 
560 


Although \esc is trying to be forgiving with the choice of sample points and to 
561 
caltinay 
2982 
convert 
562 
ahallam 
2979 
where necessary the adjustment of the function space is not always possible. So 
563 
caltinay 
2982 
it is advisable to make a careful choice on the function space used. 
564 
gross 
2905 

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


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


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


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


distribution and \verb totT the total heat in the system. 
570 
gross 
2905 
\begin{python} 
571 
gross 
2878 
while t < tend: 
572 


i+=1 #increment the counter 
573 


t+=h #increment the current time 
574 


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


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


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


\end{python} 
578 
ahallam 
2979 
The last statement in this script calculates the total energy in the system as 
579 
jfenwick 
3308 
the volume integral of $\rho c_{p} T$ over the block. 
580 
caltinay 
2982 
As the blocks are insulated no energy should be lost or added. 
581 
gross 
2905 
The total energy should stay constant for the example discussed here. 
582 
ahallam 
2401 

583 
gross 
2905 
\subsection{Running the Script} 
584 
ahallam 
2975 
The script presented so far is available under 
585 
gross 
2949 
\verbexample01a.py. You can edit this file with your favourite text editor. 
586 
caltinay 
2982 
On most operating systems\footnote{The \texttt{runescript} launcher is not 
587 


supported under {\it MS Windows} yet.} you can use the 
588 
ahallam 
2979 
\program{runescript} command 
589 
gross 
2905 
to launch {\it escript} scripts. For the example script use; 
590 


\begin{verbatim} 
591 
gross 
2949 
runescript example01a.py 
592 
gross 
2905 
\end{verbatim} 
593 


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


the python interpreter directly; 
595 


\begin{verbatim} 
596 
gross 
2949 
python example01a.py 
597 
gross 
2905 
\end{verbatim} 
598 
caltinay 
2982 
if the system is configured correctly (please talk to your system 
599 
ahallam 
2979 
administrator). 
600 
gross 
2905 

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

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


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


Two types will be demonstrated in this cookbook; 
607 


\mpl\footnote{\url{http://matplotlib.sourceforge.net/}} and 
608 
caltinay 
2982 
\verbVTK\footnote{\url{http://www.vtk.org/}}. 
609 
ahallam 
2979 
The \mpl package is a component of SciPy\footnote{\url{http://www.scipy.org}} 
610 


and is good for basic graphs and plots. 
611 
caltinay 
2982 
For more complex visualisation tasks, in particular two and three dimensional 
612 
ahallam 
2979 
problems we recommend the use of more advanced tools. For instance, \mayavi 
613 


\footnote{\url{http://code.enthought.com/projects/mayavi/}} 
614 
ahallam 
2975 
which is based upon the \verbVTK toolkit. The usage of \verbVTK based 
615 
caltinay 
2982 
visualisation is discussed in Chapter~\ref{Sec:2DHD} which focuses on a two 
616 
ahallam 
2979 
dimensional PDE. 
617 
gross 
2878 

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


are interested in showing the 
620 


behaviour of the total energy over time and secondly, how the temperature 
621 
caltinay 
2982 
distribution within the block is developing over time. 
622 


Let us start with the first task. 
623 
gross 
2878 

624 
caltinay 
2982 
The idea is to create a record of the time marks and the corresponding total 
625 
ahallam 
2979 
energies observed. 
626 
gross 
2878 
\pyt provides the concept of lists for this. Before 
627 
ahallam 
2979 
the time loop is opened we create empty lists for the time marks \verbt_list 
628 


and the total energies \verbE_list. 
629 


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


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


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


modifications our script looks as follows: 
633 
ahallam 
2775 
\begin{python} 
634 
gross 
2878 
t_list=[] 
635 


E_list=[] 
636 


# ... start iteration: 
637 


while t<tend: 
638 


t+=h 
639 


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


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


totE=integrate(rhocp*T) 
642 


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


E_list.append(totE) # add current total energy to record 
644 
ahallam 
2775 
\end{python} 
645 
ahallam 
2979 
To plot $t$ over $totE$ we use \mpl a module contained within \pylab which needs 
646 
caltinay 
2982 
to be loaded before use; 
647 
ahallam 
2775 
\begin{python} 
648 
gross 
2878 
import pylab as pl # plotting package. 
649 
ahallam 
2775 
\end{python} 
650 
caltinay 
2982 
Here we are not using \verbfrom pylab import * in order to avoid name 
651 


clashes for function names within \esc. 
652 
ahallam 
2401 

653 
ahallam 
2979 
The following statements are added to the script after the time loop has been 
654 


completed; 
655 
ahallam 
2775 
\begin{python} 
656 
gross 
2878 
pl.plot(t_list,E_list) 
657 


pl.title("Total Energy") 
658 
gross 
2905 
pl.axis([0,max(t_list),0,max(E_list)*1.1]) 
659 
gross 
2878 
pl.savefig("totE.png") 
660 
ahallam 
2775 
\end{python} 
661 
ahallam 
2979 
The first statement hands over the time marks and corresponding total energies 
662 
caltinay 
2982 
to the plotter. 
663 


The second statement sets the title for the plot. The third statement 
664 
ahallam 
2979 
sets the axis ranges. In most cases these are set appropriately by the plotter. 
665 



666 
caltinay 
2982 
The last statement generates the plot and writes the result into the file 
667 


\verbtotE.png which can be displayed by (almost) any image viewer. 
668 
ahallam 
2979 
As expected the total energy is constant over time, see 
669 


\reffig{fig:onedheatout1}. 
670 
ahallam 
2401 

671 
ahallam 
3370 
\begin{figure}[ht] 
672 
ahallam 
3029 
\begin{center} 
673 


\includegraphics[width=4in]{figures/ttblockspyplot150} 
674 


\caption{Example 1b: Total Energy in the Blocks over Time (in seconds)} 
675 


\label{fig:onedheatout1} 
676 


\end{center} 
677 


\end{figure} 
678 


\clearpage 
679 



680 
gross 
2878 
\subsection{Plotting the Temperature Distribution} 
681 
gross 
2931 
\label{sec: plot T} 
682 
gross 
2949 
\sslist{example01c.py} 
683 
ahallam 
2979 
For plotting the spatial distribution of the temperature we need to modify the 
684 
caltinay 
2982 
strategy we have used for the total energy. 
685 


Instead of producing a final plot at the end we will generate a 
686 
ahallam 
2979 
picture at each time step which can be browsed as a slide show or composed into 
687 


a movie. 
688 


The first problem we encounter is that if we produce an image at each time step 
689 
caltinay 
2982 
we need to make sure that the images previously generated are not overwritten. 
690 
gross 
2878 

691 
ahallam 
2979 
To develop an incrementing file name we can use the following convention. It is 
692 
caltinay 
2982 
convenient to put all image files showing the same variable  in our case the 
693 


temperature distribution  into a separate directory. 
694 


As part of the \verbos module\footnote{The \texttt{os} module provides 
695 
ahallam 
2979 
a powerful interface to interact with the operating system, see 
696 


\url{http://docs.python.org/library/os.html}.} \pyt 
697 
caltinay 
2982 
provides the \verbos.path.join command to build file and directory names in a 
698 


platform independent way. Assuming that 
699 


\verbsave_path is the name of the directory we want to put the results in the 
700 


command is; 
701 
ahallam 
2775 
\begin{python} 
702 
gross 
2878 
import os 
703 


os.path.join(save_path, "tempT%03d.png"%i ) 
704 
ahallam 
2775 
\end{python} 
705 
gross 
2878 
where \verbi is the time step counter. 
706 
caltinay 
2982 
There are two arguments to the \verbjoin command. The \verbsave_path 
707 


variable is a predefined string pointing to the directory we want to save our 
708 


data, for example a single subfolder called \verbdata would be defined by; 
709 
gross 
2878 
\begin{verbatim} 
710 


save_path = "data" 
711 


\end{verbatim} 
712 
caltinay 
2982 
while a subfolder of \verbdata called \verbexample01 would be defined by; 
713 
gross 
2878 
\begin{verbatim} 
714 
gross 
2949 
save_path = os.path.join("data","example01") 
715 
gross 
2878 
\end{verbatim} 
716 
caltinay 
2982 
The second argument of \verbjoin contains a string which is the file 
717 
ahallam 
2979 
name or subdirectory name. We can use the operator \verb% to use the value of 
718 
caltinay 
2982 
\verbi as part of our filename. The substring \verb%03d indicates that we 
719 


want to substitute a value into the name; 
720 
gross 
2878 
\begin{itemize} 
721 
jfenwick 
2961 
\item \verb 0 means that small numbers should have leading zeroes; 
722 
ahallam 
2979 
\item \verb 3 means that numbers should be written using at least 3 digits; 
723 


and 
724 
caltinay 
2982 
\item \verb d means that the value to substitute will be a decimal integer. 
725 
gross 
2878 
\end{itemize} 
726 
ahallam 
3370 

727 
caltinay 
2982 
To actually substitute the value of \verbi into the name write \verb%i after 
728 
ahallam 
2979 
the string. 
729 
caltinay 
2982 
When done correctly, the output files from this command will be placed in the 
730 
ahallam 
2979 
directory defined by \verb save_path as; 
731 
gross 
2878 
\begin{verbatim} 
732 
jfenwick 
2961 
blockspyplot001.png 
733 


blockspyplot002.png 
734 


blockspyplot003.png 
735 
gross 
2878 
... 
736 


\end{verbatim} 
737 


and so on. 
738 



739 


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


\begin{verbatim} 
741 


mkDir(save_path) 
742 


\end{verbatim} 
743 
caltinay 
2982 
will check for the existence of \verb save_path and if missing, create the 
744 
ahallam 
2979 
required directories. 
745 
gross 
2878 

746 
artak 
2964 
We start by modifying our solution script. 
747 
caltinay 
2982 
Prior to the \verbwhile loop we need to extract our finite solution 
748 
ahallam 
2979 
points to a data object that is compatible with \mpl. First we create the node 
749 


coordinates of the sample points used to represent 
750 


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


plotting function. 
752 


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


\verbSolution(blocks).getX() into a \pyt list 
754 
jfenwick 
3308 
and then to a \numpy array. The $x_{0}$ component is then extracted via 
755 
ahallam 
2979 
an array slice to the variable \verbplx; 
756 
ahallam 
2775 
\begin{python} 
757 
gross 
2878 
import numpy as np # array package. 
758 


#convert solution points for plotting 
759 
gross 
2905 
plx = x.toListOfTuples() 
760 
gross 
2878 
plx = np.array(plx) # convert to tuple to numpy array 
761 


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

764 
ahallam 
2645 
\begin{figure} 
765 


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


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


\includegraphics[width=4in]{figures/blockspyplot200} 
769 
ahallam 
2979 
\caption{Example 1c: Temperature ($T$) distribution in the blocks at time steps 
770 
caltinay 
2982 
$1$, $50$ and $200$} 
771 
ahallam 
2645 
\label{fig:onedheatout} 
772 


\end{center} 
773 


\end{figure} 
774 
ahallam 
3029 
\clearpage 
775 
ahallam 
2645 

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


energy over time. 
778 


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


each to a file. 
780 
caltinay 
2982 
The following is appended to the end of the \verbwhile loop and creates one 
781 
ahallam 
2979 
figure of the temperature distribution. We start by converting the solution to a 
782 
caltinay 
2982 
tuple and then plotting this against our \textit{x coordinates} \verbplx we 
783 
ahallam 
2979 
have generated before. We add a title to the diagram before it is rendered into 
784 


a file. 
785 


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


following iteration. 
787 
ahallam 
2775 
\begin{python} 
788 
gross 
2878 
# ... start iteration: 
789 


while t<tend: 
790 
ahallam 
3373 
i+=1 
791 


t+=h 
792 


mypde.setValue(Y=qH+rhocp/h*T) 
793 


T=mypde.getSolution() 
794 


totE=integrate(rhocp*T) 
795 


print "time step %s at t=%e days completed. total energy = %e."%(i,t/day,totE) 
796 


t_list.append(t) 
797 


E_list.append(totE) 
798 



799 


#establish figure 1 for temperature vs x plots 
800 


tempT = T.toListOfTuples() 
801 


pl.figure(1) #current figure 
802 


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


# add title 
804 


pl.axis([0,mx,T1*.9,T2*1.1]) 
805 


pl.title("Temperature across blocks at time %d days"%(t/day)) 
806 


#save figure to file 
807 


pl.savefig(os.path.join(save_path,"tempT", "blockspyplot%03d.png"%i)) 
808 


pl.clf() #clear figure 
809 


\end{python} 
810 
gross 
2905 
Some results are shown in \reffig{fig:onedheatout}. 
811 
jfenwick 
2657 

812 
caltinay 
2982 
\subsection{Making a Video} 
813 
ahallam 
2979 
Our saved plots from the previous section can be cast into a video using the 
814 


following command appended to the end of the script. The \verb mencoder command 
815 
caltinay 
2982 
is not available on every platform, so some users need to use an alternative 
816 


video encoder. 
817 
ahallam 
2775 
\begin{python} 
818 
caltinay 
2982 
# compile the *.png files to create a *.avi video that shows T change 
819 
gross 
2878 
# with time. This operation uses Linux mencoder. For other operating 
820 
gross 
2905 
# systems it is possible to use your favourite video compiler to 
821 
ahallam 
2606 
# convert image files to videos. 
822 
gross 
2477 

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