1 
% $Id$ 
2 
\chapter{How to Solve The Diffusion Equation} 
3 
\label{DIFFUSION CHAP} 
4 

5 
\begin{figure} 
6 
\centerline{\includegraphics[width=\figwidth]{DiffusionDomain}} 
7 
\caption{Temperature Diffusion Problem with Circular Heat Source} 
8 
\label{DIFFUSION FIG 1} 
9 
\end{figure} 
10 

11 
\begin{figure} 
12 
\centerline{\includegraphics[width=\figwidth]{DiffusionRes1}} 
13 
\centerline{\includegraphics[width=\figwidth]{DiffusionRes16}} 
14 
\centerline{\includegraphics[width=\figwidth]{DiffusionRes32}} 
15 
\centerline{\includegraphics[width=\figwidth]{DiffusionRes48}} 
16 
\caption{Results of the Temperture Diffusion Problem for Time Steps $1$ $16$, $32$ and $48$.} 
17 
\label{DIFFUSION FIG 2} 
18 
\end{figure} 
19 

20 

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

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

34 
\section{\label{DIFFUSION TEMP SEC}Temperature Diffusion} 
35 

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

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

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

82 
The temperature is calculated discrete time nodes $t^{(n)}$ where 
83 
$t^{(0)}=0$ and $t^{(n)}=t^{(n1)}+h$ where $h>0$ is the step size which is assumed to be constant. 
84 
In the following the upper index ${(n)}$ is refering to a value at time $t^{(n)}$. The simplest 
85 
and most robust scheme to approximate the time derivative of the the temperature is the 
86 
\index{backward Euler} scheme, see~\cite{XXX} for alternatives. The backward Euler scheme bases 
87 
on the Taylor expansion 
88 
\begin{equation} 
89 
T^{(n1)}\approx T^{(n)}+T\hackscore{,t}^{(n)}(t^{(n1)}t^{(n)}) 
90 
=T^{(n1)} + h \cdot T\hackscore{,t}^{(n)} 
91 
\label{DIFFUSION TEMP EQ 6} 
92 
\end{equation} 
93 
which 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 + \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 from \eqn{DIFFUSION TEMP EQ 2}. 
101 
this forms a boundary value problem that has to be solved for each time step. 
102 
As a first step to implement a solver for the temperature diffusion problem we will 
103 
first implement a solver for the boundary value problem that has to be solved at each time step. 
104 

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

126 
We want to use the \LinearPDE class provided by \escript to define and solve a general linear PDE such as the 
127 
Helmholtz equation. We have used a special case of the \LinearPDE class, namely the 
128 
\Poisson class already in \Chap{FirstSteps}. 
129 
Here we will write our own specialized class of the \LinearPDE to solve the Helmholtz equation. 
130 

131 
The general form of a single PDE that can be handeled by the \LinearPDE class is 
132 
\begin{equation}\label{EQU.FEM.1} 
133 
(A\hackscore{jl} u\hackscore{,l})\hackscore{,j}+D u =Y \; . 
134 
\end{equation} 
135 
The general form and systems is discussed in \Sec{SEC LinearPDE}. 
136 
$A$, $D$ and $Y$ are the known coeffecients of the PDE \index{partial differential equation!coefficients}. 
137 
Notice that $A$ is a matrix or tensor of order 2 and $D$ and $Y$ are scalar. 
138 
They may be constant or may depend on their 
139 
location in the domain but must not depend on the unknown solution $u$. 
140 
The following natural boundary conditions \index{boundary condition!natural} that 
141 
are used in the \LinearPDE class have the form 
142 
\begin{equation}\label{EQU.FEM.2} 
143 
n\hackscore{j}A\hackscore{jl} u\hackscore{,l}+du=y \;. 
144 
\end{equation} 
145 
where, as usual, $n$ denotes the outer normal field on the surface of the domain. Notice that 
146 
the coefficient $A$ is already used in the PDE in \eqn{EQU.FEM.1}. $d$ and $y$ are given scalar coefficients. 
147 

148 
By inpecting the Helmholtz equation \index{Helmholtz equation} 
149 
we can easily assign values to the coefficients in the 
150 
general PDE of the \LinearPDE class: 
151 
\begin{equation}\label{DIFFUSION HELM EQ 3} 
152 
\begin{array}{llllll} 
153 
A\hackscore{ij}=\kappa \delta\hackscore{ij} & D=\omega & Y=f \\ 
154 
d=\eta & y= g & \\ 
155 
\end{array} 
156 
\end{equation} 
157 
$\delta\hackscore{ij}$ is the Kronecker symbol \index{Kronecker symbol} defined by $\delta=\hackscore{ij}=1$ for 
158 
$i=j$ and $0$ otherwise. 
159 

160 
We want to implement a 
161 
new class which we will call \class{Helmholtz} that provides the same methods like the \LinearPDE class but 
162 
is defined through the coefficeints $\kappa$, $\omega$, $f$, $\eta$, 
163 
$g$ rather than the general form given by \eqn{EQU.FEM.1}. 
164 
Python's 
165 
mechanism of inhertence allows doing this in a very easy way. 
166 
The advantage is that our new \class{Helmholtz} can be used in any context 
167 
that works with a \LinearPDE class but with an easier interface to define the PDE. 
168 
This improves reuasablity as well as maintainability of program codes. 
169 

170 
We want to implement a 
171 
new class which we will call \class{Helmholtz} that provides the same methods like the \LinearPDE class but 
172 
is defined through the coefficeints $\kappa$, $\omega$, $f$, $\alpha$, 
173 
$g$ rather than the general form given by \eqn{EQU.FEM.1}. 
174 
Python's mechanism of subclasses allows doing this in a very easy way. 
175 
The \Poisson class of the \linearPDEsPack module, 
176 
which we have already used in \Chap{FirstSteps}, is in fact a subclass of the 
177 
\LinearPDE class. That means that it all methods (such as the \method{getSolution}) 
178 
from the parent class \LinearPDE are defined for any \Poisson object. However, new 
179 
methods can be added and methods of the parent class can be redefined. In fact, 
180 
the \Poisson class redefines the \method{setValue} of the \LinearPDE class which is called to assign 
181 
values to the coefficients of the PDE. This is exactly what we will do when we define 
182 
our new \class{Helmholtz} class: 
183 
\begin{python} 
184 
from esys.linearPDEs import LinearPDE 
185 
import numarray 
186 
class Helmholtz(LinearPDE) 
187 
def setValue(self,kappa=0,omega=1,f=0,eta=0,g=0) 
188 
self._setValue(A=kappa*numarray.identity(self.getDim()),D=omega,Y=f,d=eta,y=g) 
189 
\end{python} 
190 
\code{class Helmholtz(linearPDE)} declares the new \class{Helmholtz} class as a subclass 
191 
of the \LinearPDE which have imported in the first line of the script. 
192 
We add the method \method{setValue} to the class which overwrites the 
193 
\method{setValue} method of the \LinearPDE class. The new methods which has the 
194 
parameters of the Helmholtz \eqn{DIFFUSION HELM EQ 1} as arguments 
195 
maps the parameters of the coefficients of the general PDE defined 
196 
in \eqn{EQU.FEM.1}. The coefficient \var{A} is defined by Kroneckers symbol. we use the 
197 
\numarray function \function{identity} which return a square matrix which has ones on the 
198 
main diagonal and zeros out side the main diagonal. The argument of \function{identity} gives the order of the matrix. 
199 
Here we use 
200 
the \method{getDim} of the \LinearPDE class object \var{self} to get the spatial dimension of the domain of the 
201 
PDE. As we will make use of the \class{Helmholtz} class several times, it is convient to 
202 
put is definition into a file which we name \file{mytools.py} available in the \ExampleDirectory. 
203 
You can use your favourite editor to create and edit the file. 
204 

205 
An object of the \class{Helmholtz} class is created through the statments: 
206 
\begin{python} 
207 
from mytools import * 
208 
mypde=Helmholtz(mydomain) 
209 
mypde.setValue(kappa=10.,omega=0.1,f=12) 
210 
u=mypde.getSolution() 
211 
\end{python} 
212 
In the first statement we import all definition from the \file{mytools.py} \index{scripts!\file{mytools.py}}. Make sure 
213 
that \file{mytools.py} is in directory from where you have started Python. 
214 
\var{mydomain} is the \Domain of the PDE. In the third statment values are 
215 
assigned to the PDE parameters. As no values for arguments \var{eta} and \var{g} are 
216 
specified the default values $0$ are used \footnote{It would be better to use the default value 
217 
\var{escript.Data()} rather then $0$ as then the coefficient would be defined as being not present and 
218 
would not be processed when the PDE is evaluated.}. In the forth statement the solution of the 
219 
PDE is returned. 
220 

221 
We want to test our \class{Helmholtz} class on a rectangular domain 
222 
of length $l$ and height $h$. We do this by choosing a simple test solution, 
223 
here we take $u=x\hackscore{0}$ and then calculate the right hand side terms $f$ and $g$ such that 
224 
the test solution becomes the solution of the problem. If we take $\kappa=1$ 
225 
an easy calculation shows that we have to choose $f=\omega \cdot x\hackscore{0}$. On the boundary we get 
226 
$\kappa n\hackscore{i} u\hackscore{,i}=n\hackscore{0}$. 
227 
So we have to set $g=n\hackscore{0}+\eta x\hackscore{0}$. The following script \file{helmholtztest.py} 
228 
\index{scripts!\file{helmholtztest.py}} which is available in the \ExampleDirectory 
229 
implements this test problem using the \finley PDE solver: 
230 
\begin{python} 
231 
from mytools import * 
232 
from esys.escript import * 
233 
import esys.finley 
234 
#... set some parameters ... 
235 
omega=0.1 
236 
eta=10. 
237 
#... generate domain ... 
238 
mydomain = esys.finley.Rectangle(l0=5.,l1=1.,n0=50, n1=10) 
239 
#... open PDE and set coefficients ... 
240 
mypde=Helmholtz(mydomain) 
241 
n=mydomain.getNormal() 
242 
x=mydomain.getX() 
243 
mypde.setValue(1,omega,omega*x[0],eta,n[0]+eta*x[0]) 
244 
#... calculate error of the PDE solution ... 
245 
u=mypde.getSolution() 
246 
print "error is ",Lsup(ux[0]) 
247 
\end{python} 
248 
The script is similar to the script \file{mypoisson.py} dicussed in \Chap{FirstSteps}. 
249 
\code{mydomain.getNormal()} returns the outer normal field on the surface of the domain. The function \function{Lsup} 
250 
is imported by the \code{from escript import *} statement and returns the maximum absulute value of it argument. To run 
251 
The error shown by the print statement should be in the order of $10^{7}$. As piecewise linear interpolation is 
252 
used to approximate the solution and our solution is a linear function of the spatial coordinates one may 
253 
expect that the error is zero. However, as most PDE packages uses an iterative solver which is terminated 
254 
when a given tolerance has been reached. The default tolerance is $10^{8}$. Thsi value can be altered by using the 
255 
\method{setTolerance} of the \LinearPDE class. 
256 

257 
\section{The Transition Problem} 
258 
\label{DIFFUSION TRANS SEC} 
259 
Now we are ready to solve the original time dependent problem. The main 
260 
part of the script is the loop over the time $t$ which takes the following form: 
261 
\begin{python} 
262 
mypde=Helmholtz(mydomain) 
263 
while t<t_end: 
264 
mypde.setValue(kappa,rhocp/h,q+rhocp/h*T,eta,eta*Tref) 
265 
T=mypde.getSolution() 
266 
t+=h 
267 
\end{python} 
268 
\var{kappa}, \var{rhocp}, \var{eta} and \var{Tref} are input parameters of the model. \var{q} is the heat source 
269 
in the domain and \var{h} is the time step size which has to be chosen. Notice that the \class{Hemholtz} 
270 
is created before the loop over time is entered while in each time step only the coefficients 
271 
are reset in each time step. This way some information about the reperesentation of the PDE can be reused 
272 
\footnote{The efficience can be improved further by setting the coefficients in the operator 
273 
\var{kappa}, \var{omega} and \var{eta} before entering the \code{while}loop and only update the coefficients 
274 
in the right hand side \var{f} and \var{g}. This needs a more careful implementation of the \method{setValue} 
275 
method but gives the advantage that the \LinearPDE class can save rebuilding the PDE operator}. The variable \var{T} 
276 
holds the current temperature. It is used to calculate the right hand side coefficient \var{f} in the 
277 
Helmholtz \eqn{DIFFUSION HELM EQ 1}. Statement \code{T=mypde.getSolution()} overwrites \var{T} with the 
278 
temperature of the new time step $\var{t}+\var{h}$. To get this iterative process going we need to sepcify the 
279 
initial temperature distribution, which equal to $T\hackscore{ref}$. 
280 

281 
The heat source \var{q} which is defined in \eqn{DIFFUSION TEMP EQ 1b} shall be \var{q0} 
282 
at an area defined as a circle of radius \var{r} and center \var{xc} and zero outside this circle. 
283 
\var{q0} is a fixed constant. The following script defines \var{q} as desired: 
284 
\begin{python} 
285 
xc=[0.02,0.002] 
286 
r=0.001 
287 
x=mydomain.getX() 
288 
q=q0*(length(xxc)r).whereNegative() 
289 
\end{python} 
290 
\var{x} is a \Data class object of 
291 
the \escript module defining the locations of points in the \Domain \var{mydomain}. 
292 
\code{length(xxc)} calculates the distances in the Euclidean norm 
293 
of the locations \var{x} to the center of the circle \var{xc} where the heat source is acting. 
294 
Notice that the coordinates of \var{xc} are defined as a list of floating point numbers. It is independently 
295 
converted into a \Data class object before subtracted from \var{x}. The method \method{whereNegative} of 
296 
a \Data class object, in this case the result of the expression 
297 
\code{length(xxc)r}, returns a \Data class which is equal one where the object negative and 
298 
zero elsewhere. After multiplication with \var{q0} we get a function with the deired property. 
299 

300 
Now we can put the components together to the script \file{diffusion.py} which is available in the \ExampleDirectory: 
301 
\index{scripts!\file{diffusion.py}}: 
302 
\begin{python} 
303 
from mytools import * 
304 
from esys.escript import * 
305 
import esys.finley 
306 
#... set some parameters ... 
307 
x_c=[0.02,0.002] 
308 
r=0.001 
309 
q0=50.e6 
310 
Tref=0. 
311 
rhocp=2.6e6 
312 
eta=75. 
313 
kappa=240. 
314 
t_end=5. 
315 
# ...time step size and counter ... 
316 
h=0.1 
317 
i=0 
318 
t=0 
319 
#... generate domain ... 
320 
mydomain = esys.finley.Rectangle(l0=0.05,l1=0.01,n0=250, n1=50) 
321 
#... open PDE ... 
322 
mypde=Helmholtz(mydomain) 
323 
# ... set heat source: .... 
324 
x=mydomain.getX() 
325 
q=q0*(length(xx_c)r).whereNegative() 
326 
# ... set initial temperature .... 
327 
T=Tref 
328 
# ... start iteration: 
329 
while t<t_end: 
330 
i+=1 
331 
t+=h 
332 
print "time step :",t 
333 
mypde.setValue(kappa=kappa,omega=rhocp/h,f=q+rhocp/h*T,eta=eta,g=eta*Tref) 
334 
T=mypde.getSolution() 
335 
T.saveDX("T%d.dx"%i) 
336 
\end{python} 
337 
The script will create the files \file{T.1.dx}, 
338 
\file{T.2.dx}, $\ldots$, \file{T.50.dx} in the directory where the script has been started. The files give the 
339 
temperature distributions at time steps $1$, $2$, $\ldots$, $50$ in the \OpenDX file format. 
340 
An easy way to visualize the results is the command 
341 
\begin{verbatim} 
342 
dx edit diffusion.net 
343 
\end{verbatim} 
344 
where \file{diffusion.net} is an \OpenDX script available in the \ExampleDirectory. 
345 
\fig{DIFFUSION FIG 2} shows the result for some selected time steps. 
346 
