1 
% $Id$ 
2 
% 
3 
% Copyright © 2006 by ACcESS MNRF 
4 
% \url{http://www.access.edu.au 
5 
% Primary Business: Queensland, Australia. 
6 
% Licensed under the Open Software License version 3.0 
7 
% http://www.opensource.org/licenses/osl3.0.php 
8 
% 
9 

10 
\section{The Diffusion Problem} 
11 
\label{DIFFUSION CHAP} 
12 

13 
\begin{figure} 
14 
\centerline{\includegraphics[width=\figwidth]{figures/DiffusionDomain.eps}} 
15 
\caption{Temperature Diffusion Problem with Circular Heat Source} 
16 
\label{DIFFUSION FIG 1} 
17 
\end{figure} 
18 

19 
\subsection{\label{DIFFUSION OUT SEC}Outline} 
20 
In this chapter we will discuss how to solve the time dependenttemperature diffusion\index{diffusion equation} for 
21 
a block of material. Within the block there is a heat source which drives the temperature diffusion. 
22 
On the surface, energy can radiate into the surrounding environment. 
23 
\fig{DIFFUSION FIG 1} shows the configuration. 
24 

25 
In the next \Sec{DIFFUSION TEMP SEC} we will present the relevant model. A 
26 
time integration scheme is introduced to calculate the temperature at given time nodes $t^{(n)}$. 
27 
We will see that at each time step a Helmholtz equation \index{Helmholtz equation} 
28 
must be solved. 
29 
The implementation of a Helmholtz equation solver will be discussed in \Sec{DIFFUSION HELM SEC}. 
30 
In Section~\ref{DIFFUSION TRANS SEC} the solver of the Helmholtz equation is used to build a 
31 
solver for the temperature diffusion problem. 
32 

33 
\subsection{\label{DIFFUSION TEMP SEC}Temperature Diffusion} 
34 
The unknown temperature $T$ is a function of its location in the domain and time $t>0$. The governing equation 
35 
in the interior of the domain is given by 
36 
\begin{equation} 
37 
\rho c\hackscore p T\hackscore{,t}  (\kappa T\hackscore{,i})\hackscore{,i} = q\hackscore H 
38 
\label{DIFFUSION TEMP EQ 1} 
39 
\end{equation} 
40 
where $\rho c\hackscore p$ and $\kappa$ are given material constants. In case of a composite 
41 
material the parameters depend on their location in the domain. $q\hackscore H$ is 
42 
a heat source (or sink) within the domain. We are using Einstein summation convention \index{summation convention} 
43 
as introduced in \Chap{FirstSteps}. In our case we assume $q\hackscore H$ to be equal to a constant heat production rate 
44 
$q^{c}$ on a circle or sphere with center $x^c$ and radius $r$ and $0$ elsewhere: 
45 
\begin{equation} 
46 
q\hackscore H(x,t)= 
47 
\left\{ 
48 
\begin{array}{lcl} 
49 
q^c & & \xx^c\ \le r \\ 
50 
& \mbox{if} \\ 
51 
0 & & \mbox{else} \\ 
52 
\end{array} 
53 
\right. 
54 
\label{DIFFUSION TEMP EQ 1b} 
55 
\end{equation} 
56 
for all $x$ in the domain and all time $t>0$. 
57 

58 
On the surface of the domain we are 
59 
specifying a radiation condition 
60 
which precribes the normal component of the flux $\kappa T\hackscore{,i}$ to be proportional 
61 
to the difference of the current temperature to the surrounding temperature $T\hackscore{ref}$: 
62 
\begin{equation} 
63 
\kappa T\hackscore{,i} n\hackscore i = \eta (T\hackscore{ref}T) 
64 
\label{DIFFUSION TEMP EQ 2} 
65 
\end{equation} 
66 
$\eta$ is a given material coefficient depending on the material of the block and the surrounding medium. 
67 
As usual $n\hackscore i$ is the $i$th component of the outer normal field \index{outer normal field} 
68 
at the surface of the domain. 
69 

70 
To solve the time dependent \eqn{DIFFUSION TEMP EQ 1} the initial temperature at time 
71 
$t=0$ has to be given. Here we assume that the initial temperature is the surrounding temperature: 
72 
\begin{equation} 
73 
T(x,0)=T\hackscore{ref} 
74 
\label{DIFFUSION TEMP EQ 4} 
75 
\end{equation} 
76 
for all $x$ in the domain. It is pointed out that 
77 
the initial conditions satisfy the 
78 
boundary condition defined by \eqn{DIFFUSION TEMP EQ 2}. 
79 

80 
The temperature is calculated at discrete time nodes $t^{(n)}$ where 
81 
$t^{(0)}=0$ and $t^{(n)}=t^{(n1)}+h$ where $h>0$ is the step size which is assumed to be constant. 
82 
In the following the upper index ${(n)}$ refers to a value at time $t^{(n)}$. The simplest 
83 
and most robust scheme to approximate the time derivative of the the temperature is 
84 
the backward Euler 
85 
\index{backward Euler} scheme. The backward Euler 
86 
scheme is based 
87 
on the Taylor expansion of $T$ at time $t^{(n)}$: 
88 
\begin{equation} 
89 
T^{(n)}\approx T^{(n1)}+T\hackscore{,t}^{(n)}(t^{(n)}t^{(n1)}) 
90 
=T^{(n1)} + h \cdot T\hackscore{,t}^{(n)} 
91 
\label{DIFFUSION TEMP EQ 6} 
92 
\end{equation} 
93 
This is inserted into \eqn{DIFFUSION TEMP EQ 1}. By separating the terms at 
94 
$t^{(n)}$ and $t^{(n1)}$ one gets for $n=1,2,3\ldots$ 
95 
\begin{equation} 
96 
\frac{\rho c\hackscore p}{h} T^{(n)}  (\kappa T^{(n)}\hackscore{,i})\hackscore{,i} = q\hackscore H + \frac{\rho c\hackscore p}{h} T^{(n1)} 
97 
\label{DIFFUSION TEMP EQ 7} 
98 
\end{equation} 
99 
where $T^{(0)}=T\hackscore{ref}$ is taken form the initial condition given by \eqn{DIFFUSION TEMP EQ 4}. 
100 
Together with the natural boundary condition 
101 
\begin{equation} 
102 
\kappa T\hackscore{,i}^{(n)} n\hackscore i = \eta (T\hackscore{ref}T^{(n)}) 
103 
\label{DIFFUSION TEMP EQ 2222} 
104 
\end{equation} 
105 
taken from \eqn{DIFFUSION TEMP EQ 2} 
106 
this forms a boundary value problem that has to be solved for each time step. 
107 
As a first step to implement a solver for the temperature diffusion problem we will 
108 
first implement a solver for the boundary value problem that has to be solved at each time step. 
109 

110 
\subsection{\label{DIFFUSION HELM SEC}Helmholtz Problem} 
111 
The partial differential equation to be solved for $T^{(n)}$ has the form 
112 
\begin{equation} 
113 
\omega T^{(n)}  (\kappa T^{(n)}\hackscore{,i})\hackscore{,i} = f 
114 
\label{DIFFUSION HELM EQ 1} 
115 
\end{equation} 
116 
and we set 
117 
\begin{equation} 
118 
\omega=\frac{\rho c\hackscore p}{h} \mbox{ and } f=q\hackscore H +\frac{\rho c\hackscore p}{h}T^{(n1)} \;. 
119 
\label{DIFFUSION HELM EQ 1b} 
120 
\end{equation} 
121 
With $g=\eta T\hackscore{ref}$ the radiation condition defined by \eqn{DIFFUSION TEMP EQ 2222} 
122 
takes the form 
123 
\begin{equation} 
124 
\kappa T^{(n)}\hackscore{,i} n\hackscore{i} = g  \eta T^{(n)}\mbox{ on } \Gamma 
125 
\label{DIFFUSION HELM EQ 2} 
126 
\end{equation} 
127 
The partial differential \eqn{DIFFUSION HELM EQ 1} together with boundary conditions of \eqn{DIFFUSION HELM EQ 2} 
128 
is called the Helmholtz equation \index{Helmholtz equation}. 
129 

130 
We want to use the \LinearPDE class provided by \escript to define and solve a general linear,steady, second order PDE such as the 
131 
Helmholtz equation. For a single PDE the \LinearPDE class supports the following form: 
132 
\begin{equation}\label{LINEARPDE.SINGLE.1 TUTORIAL} 
133 
(A\hackscore{jl} u\hackscore{,l})\hackscore{,j}+D u = Y \; . 
134 
\end{equation} 
135 
where we show only the coefficients relevant for the problem discussed here. For the general form of 
136 
single PDE see \eqn{LINEARPDE.SINGLE.1}. 
137 
The coefficients $A$, and $Y$ have to be specified through \Data objects in the 
138 
\Function on the PDE or objects that can be converted into such \Data objects. 
139 
$A$ is a \RankTwo and $D$ and $Y$ are scalar. 
140 
The following natural 
141 
boundary conditions are considered \index{boundary condition!natural} on $\Gamma$: 
142 
\begin{equation}\label{LINEARPDE.SINGLE.2 TUTORIAL} 
143 
n\hackscore{j}A\hackscore{jl} u\hackscore{,l}+d u= y \;. 
144 
\end{equation} 
145 
Notice that the coefficient $A$ is the same like in the PDE~\eqn{LINEARPDE.SINGLE.1 TUTORIAL}. 
146 
The coefficients $d$ and $y$ are 
147 
each a \Scalar in the \FunctionOnBoundary. Constraints \index{constraint} for the solution prescribing the value of the 
148 
solution at certain locations in the domain. They have the form 
149 
\begin{equation}\label{LINEARPDE.SINGLE.3 TUTORIAL} 
150 
u=r \mbox{ where } q>0 
151 
\end{equation} 
152 
$r$ and $q$ are each \Scalar where $q$ is the characteristic function 
153 
\index{characteristic function} defining where the constraint is applied. 
154 
The constraints defined by \eqn{LINEARPDE.SINGLE.3 TUTORIAL} override any other condition set by 
155 
\eqn{LINEARPDE.SINGLE.1 TUTORIAL} or \eqn{LINEARPDE.SINGLE.2 TUTORIAL}. 
156 
The \Poisson class of the \linearPDEs module, 
157 
which we have already used in \Chap{FirstSteps}, is in fact a subclass of the more general 
158 
\LinearPDE class. The \linearPDEs module provides a \Helmholtz class but 
159 
we will make direct use of the general \LinearPDE class. 
160 

161 
By inspecting the Helmholtz equation \index{Helmholtz equation} 
162 
(\ref{DIFFUSION HELM EQ 1}) and boundary condition (\ref{DIFFUSION HELM EQ 2}) and 
163 
substituting $u$ for $T^{(n)}$ 
164 
we can easily assign values to the coefficients in the 
165 
general PDE of the \LinearPDE class: 
166 
\begin{equation}\label{DIFFUSION HELM EQ 3} 
167 
\begin{array}{llllll} 
168 
A\hackscore{ij}=\kappa \delta\hackscore{ij} & D=\omega & Y=f \\ 
169 
d=\eta & y= g & \\ 
170 
\end{array} 
171 
\end{equation} 
172 
$\delta\hackscore{ij}$ is the Kronecker symbol \index{Kronecker symbol} defined by $\delta\hackscore{ij}=1$ for 
173 
$i=j$ and $0$ otherwise. Undefined coefficients are assumed to be not present.\footnote{There is a difference 
174 
in \escript of being not present and set to zero. As not present coefficients are not processed, 
175 
it is more efficient to leave a coefficient undefined instead of assigning zero to it.} 
176 
In this diffusion example we do not need to define a characteristic function $q$ because the 
177 
boundary conditions we consider in \eqn{DIFFUSION HELM EQ 2} are just the natural boundary 
178 
conditions which are already defined in the \LinearPDE class (shown in \eqn{LINEARPDE.SINGLE.2 TUTORIAL}). 
179 

180 
Defining and solving the Helmholtz equation is very easy now: 
181 
\begin{python} 
182 
from esys.escript import * 
183 
from linearPDEs import LinearPDE 
184 
mypde=LinearPDE(mydomain) 
185 
mypde.setValue(A=kappa*kronecker(mydomain),D=omega,Y=f,d=eta,y=g) 
186 
u=mypde.getSolution() 
187 
\end{python} 
188 
where we assume that \code{mydomain} is a \Domain object and 
189 
\code{kappa} \code{omega} \code{eta} and \code{g} are given scalar values 
190 
typically \code{float} or \Data objects. The \method{setValue} method 
191 
assigns values to the coefficients of the general PDE. The \method{getSolution} method solves 
192 
the PDE and returns the solution \code{u} of the PDE. \function{kronecker} is \escript function 
193 
returning the Kronecker symbol. 
194 

195 
The coefficients can set by several calls of \method{setValue} where the order can be chosen arbitrarily. 
196 
If a value is assigned to a coefficient several times, the last assigned value is used when 
197 
the solution is calculated: 
198 
\begin{python} 
199 
mypde=LinearPDE(mydomain) 
200 
mypde.setValue(A=kappa*kronecker(mydomain),d=eta) 
201 
mypde.setValue(D=omega,Y=f,y=g) 
202 
mypde.setValue(d=2*eta) # overwrites d=eta 
203 
u=mypde.getSolution() 
204 
\end{python} 
205 
In some cases the solver of the PDE can make use of the fact that the PDE is symmetric\index{symmetric PDE} where the 
206 
PDE is called symmetric if 
207 
\begin{equation}\label{LINEARPDE.SINGLE.4 TUTORIAL} 
208 
A\hackscore{jl}=A\hackscore{lj}\;. 
209 
\end{equation} 
210 
Note that $D$ and $d$ may have any value and the right hand sides $Y$, $y$ as well as the constraints 
211 
are not relevant. The Helmholtz problem is symmetric. 
212 
The \LinearPDE class provides the method \method{checkSymmetry} method to check if the given PDE is symmetric. 
213 
\begin{python} 
214 
mypde=LinearPDE(mydomain) 
215 
mypde.setValue(A=kappa*kronecker(mydomain),d=eta) 
216 
print mypde.checkSymmetry() # returns True 
217 
mypde.setValue(B=kronecker(mydomain)[0]) 
218 
print mypde.checkSymmetry() # returns False 
219 
mypde.setValue(C=kronecker(mydomain)[0]) 
220 
print mypde.checkSymmetry() # returns True 
221 
\end{python} 
222 
Unfortunately, a \method{checkSymmetry} is very expensive and is recommendable to use for 
223 
testing and debugging purposes only. The \method{setSymmetryOn} method is used to 
224 
declare a PDE symmetric: 
225 
\begin{python} 
226 
mypde = LinearPDE(mydomain) 
227 
mypde.setValue(A=kappa*kronecker(mydomain)) 
228 
mypde.setSymmetryOn() 
229 
\end{python} 
230 
Now we want to see how we actually solve the Helmholtz equation. 
231 
on a rectangular domain 
232 
of length $l\hackscore{0}=5$ and height $l\hackscore{1}=1$. We choose a simple test solution such that we 
233 
can verify the returned solution against the exact answer. Actually, we 
234 
take $T=x\hackscore{0}$ (here $q\hackscore H = 0$) and then calculate the right hand side terms $f$ and $g$ such that 
235 
the test solution becomes the solution of the problem. If we assume $\kappa$ as being constant, 
236 
an easy calculation shows that we have to choose $f=\omega \cdot x\hackscore{0}$. On the boundary we get 
237 
$\kappa n\hackscore{i} u\hackscore{,i}=\kappa n\hackscore{0}$. 
238 
So we have to set $g=\kappa n\hackscore{0}+\eta x\hackscore{0}$. The following script \file{helmholtz.py} 
239 
\index{scripts!\file{helmholtz.py}} which is available in the \ExampleDirectory 
240 
implements this test problem using the \finley PDE solver: 
241 
\begin{python} 
242 
from esys.escript import * 
243 
from esys.escript.linearPDEs import LinearPDE 
244 
from esys.finley import Rectangle 
245 
#... set some parameters ... 
246 
kappa=1. 
247 
omega=0.1 
248 
eta=10. 
249 
#... generate domain ... 
250 
mydomain = Rectangle(l0=5.,l1=1.,n0=50, n1=10) 
251 
#... open PDE and set coefficients ... 
252 
mypde=LinearPDE(mydomain) 
253 
mypde.setSymmetryOn() 
254 
n=mydomain.getNormal() 
255 
x=mydomain.getX() 
256 
mypde.setValue(A=kappa*kronecker(mydomain),D=omega,Y=omega*x[0], \ 
257 
d=eta,y=kappa*n[0]+eta*x[0]) 
258 
#... calculate error of the PDE solution ... 
259 
u=mypde.getSolution() 
260 
print "error is ",Lsup(ux[0]) 
261 
saveVTK("x0.xml",sol=u) 
262 
\end{python} 
263 
To visualise the solution `x0.~xml' just use the command 
264 
\begin{python} 
265 
mayavi d u.xml m SurfaceMap & 
266 
\end{python} 
267 
and it is easy to see that the solution $T=x\hackscore{0}$ is calculated. 
268 

269 
The script is similar to the script \file{poisson.py} dicussed in \Chap{FirstSteps}. 
270 
\code{mydomain.getNormal()} returns the outer normal field on the surface of the domain. The function \function{Lsup} 
271 
imported by the \code{from esys.escript import *} statement and returns the maximum absolute value of its argument. 
272 
The error shown by the print statement should be in the order of $10^{7}$. As piecewise bilinear interpolation is 
273 
used by \finley approximate the solution and our solution is a linear function of the spatial coordinates one might 
274 
expect that the error would be zero or in the order of machine precision (typically $\approx 10^{15}$). 
275 
However most PDE packages use an iterative solver which is terminated 
276 
when a given tolerance has been reached. The default tolerance is $10^{8}$. This value can be altered by using the 
277 
\method{setTolerance} of the \LinearPDE class. 
278 

279 
\subsection{The Transition Problem} 
280 
\label{DIFFUSION TRANS SEC} 
281 
Now we are ready to solve the original time dependent problem. The main 
282 
part of the script is the loop over time $t$ which takes the following form: 
283 
\begin{python} 
284 
t=0 
285 
T=Tref 
286 
mypde=LinearPDE(mydomain) 
287 
mypde.setValue(A=kappa*kronecker(mydomain),D=rhocp/h,d=eta,y=eta*Tref) 
288 
while t<t_end: 
289 
mypde.setValue(Y=q+rhocp/h*T) 
290 
T=mypde.getSolution() 
291 
t+=h 
292 
\end{python} 
293 
\var{kappa}, \var{rhocp}, \var{eta} and \var{Tref} are input parameters of the model. \var{q} is the heat source 
294 
in the domain and \var{h} is the time step size. 
295 
The variable \var{T} 
296 
holds the current temperature. It is used to calculate the right hand side coefficient \var{f} in the 
297 
Helmholtz equation in \eqn{DIFFUSION HELM EQ 1}. The statement \code{T=mypde.getSolution()} overwrites \var{T} with the 
298 
temperature of the new time step $\var{t}+\var{h}$. To get this iterative process going we need to specify the 
299 
initial temperature distribution, which equal to $T\hackscore{ref}$. 
300 
The \LinearPDE class object \var{mypde} 
301 
and coefficients not changing over time are set up before the loop over time is entered. In each time step only the coefficient 
302 
$Y$ is reset as it depends on the temperature of the previous time step. This allows the PDE solver to reuse information 
303 
from previous time steps as much as possible. 
304 

305 
The heat source $q\hackscore H$ which is defined in \eqn{DIFFUSION TEMP EQ 1b} is \var{qc} 
306 
in an area defined as a circle of radius \var{r} and center \var{xc} and zero outside this circle. 
307 
\var{q0} is a fixed constant. The following script defines $q\hackscore H$ as desired: 
308 
\begin{python} 
309 
from esys.escript import length,whereNegative 
310 
xc=[0.02,0.002] 
311 
r=0.001 
312 
x=mydomain.getX() 
313 
qH=q0*whereNegative(length(xxc)r) 
314 
\end{python} 
315 
\var{x} is a \Data class object of 
316 
the \escript module defining locations in the \Domain \var{mydomain}. 
317 
The \function{length()} imported from the \escript module returns the 
318 
Euclidean norm: 
319 
\begin{equation}\label{DIFFUSION HELM EQ 3aba} 
320 
\y\=\sqrt{ 
321 
y\hackscore{i} 
322 
y\hackscore{i} 
323 
} = \function{esys.escript.length}(y) 
324 
\end{equation} 
325 
So \code{length(xxc)} calculates the distances 
326 
of the location \var{x} to the center of the circle \var{xc} where the heat source is acting. 
327 
Note that the coordinates of \var{xc} are defined as a list of floating point numbers. It is independently 
328 
converted into a \Data class object before being subtracted from \var{x}. The function \function{whereNegative} 
329 
applied to 
330 
\code{length(xxc)r}, returns a \Data class which is equal to one where the object is negative and 
331 
zero elsewhere. After multiplication with \var{qc} we get a function with the desired property. 
332 

333 
Now we can put the components together to create the script \file{diffusion.py} which is available in the \ExampleDirectory: 
334 
\index{scripts!\file{diffusion.py}}: 
335 
\begin{python} 
336 
from esys.escript import * 
337 
from esys.escript.linearPDEs import LinearPDE 
338 
from esys.finley import Rectangle 
339 
#... set some parameters ... 
340 
xc=[0.02,0.002] 
341 
r=0.001 
342 
qc=50.e6 
343 
Tref=0. 
344 
rhocp=2.6e6 
345 
eta=75. 
346 
kappa=240. 
347 
tend=5. 
348 
# ... time, time step size and counter ... 
349 
t=0 
350 
h=0.1 
351 
i=0 
352 
#... generate domain ... 
353 
mydomain = Rectangle(l0=0.05,l1=0.01,n0=250, n1=50) 
354 
#... open PDE ... 
355 
mypde=LinearPDE(mydomain) 
356 
mypde.setSymmetryOn() 
357 
mypde.setValue(A=kappa*kronecker(mydomain),D=rhocp/h,d=eta,y=eta*Tref) 
358 
# ... set heat source: .... 
359 
x=mydomain.getX() 
360 
qH=qc*whereNegative(length(xxc)r) 
361 
# ... set initial temperature .... 
362 
T=Tref 
363 
# ... start iteration: 
364 
while t<tend: 
365 
i+=1 
366 
t+=h 
367 
print "time step :",t 
368 
mypde.setValue(Y=qH+rhocp/h*T) 
369 
T=mypde.getSolution() 
370 
saveVTK("T.%d.xml"%i,temp=T) 
371 
\end{python} 
372 
The script will create the files \file{T.1.xml}, 
373 
\file{T.2.xml}, $\ldots$, \file{T.50.xml} in the directory where the script has been started. The files give the 
374 
temperature distributions at time steps $1$, $2$, $\ldots$, $50$ in the \VTK file format. 
375 

376 
\begin{figure} 
377 
\centerline{\includegraphics[width=\figwidth]{figures/DiffusionRes1.eps}} 
378 
\centerline{\includegraphics[width=\figwidth]{figures/DiffusionRes16.eps}} 
379 
\centerline{\includegraphics[width=\figwidth]{figures/DiffusionRes32.eps}} 
380 
\centerline{\includegraphics[width=\figwidth]{figures/DiffusionRes48.eps}} 
381 
\caption{Results of the Temperature Diffusion Problem for Time Steps $1$ $16$, $32$ and $48$.} 
382 
\label{DIFFUSION FIG 2} 
383 
\end{figure} 
384 
An easy way to visualize the results is the command 
385 
\begin{verbatim} 
386 
mayavi d T.1.xml m SurfaceMap & 
387 
\end{verbatim} 
388 
Use the \texttt{Configure Data} 
389 
to move forward and and backwards in time. 
390 
\fig{DIFFUSION FIG 2} shows the result for some selected time steps. 