1 
% $Id$ 
2 
\section{The Diffusion Problem} 
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 
\subsection{\label{DIFFUSION OUT SEC}Outline} 
12 
In this chapter we will discuss how to solve the time dependenttemperature diffusion\index{diffusion equation} for 
13 
a block of material. Within the block there is a heat source which drives the temperature diffusion. 
14 
On the surface, energy can radiate into the surrounding environment. 
15 
\fig{DIFFUSION FIG 1} shows the configuration. 
16 

17 
In the next \Sec{DIFFUSION TEMP SEC} we will present the relevant model. A 
18 
time integration scheme is introduced to calculate the temperature at given time nodes $t^{(n)}$. 
19 
We will see that at each time step a Helmholtz equation \index{Helmholtz equation} 
20 
must be solved. 
21 
The implementation of a Helmholtz equation solver will be discussed in \Sec{DIFFUSION HELM SEC}. 
22 
In Section~\ref{DIFFUSION TRANS SEC} the solver of the Helmholtz equation is used to build a 
23 
solver for the temperature diffusion problem. 
24 

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

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

62 
To solve the time dependent \eqn{DIFFUSION TEMP EQ 1} the initial temperature at time 
63 
$t=0$ has to be given. Here we assume that the initial temperature is the surrounding temperature: 
64 
\begin{equation} 
65 
T(x,0)=T\hackscore{ref} 
66 
\label{DIFFUSION TEMP EQ 4} 
67 
\end{equation} 
68 
for all $x$ in the domain. It is pointed out that 
69 
the initial conditions satisfy the 
70 
boundary condition defined by \eqn{DIFFUSION TEMP EQ 2}. 
71 

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

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

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

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

169 
Defining and solving the Helmholtz equation is very easy now: 
170 
\begin{python} 
171 
from esys.escript import * 
172 
from linearPDEs import LinearPDE 
173 
mypde=LinearPDE(mydomain) 
174 
mypde.setValue(A=kappa*kronecker(mydomain),D=omega,Y=f,d=eta,y=g) 
175 
u=mypde.getSolution() 
176 
\end{python} 
177 
where we assume that \code{mydomain} is a \Domain object and 
178 
\code{kappa} \code{omega} \code{eta} and \code{g} are given scalar values 
179 
typically \code{float} or \Data objects. The \method{setValue} method 
180 
assigns values to the coefficients of the general PDE. The \method{getSolution} method solves 
181 
the PDE and returns the solution \code{u} of the PDE. \function{kronecker} is \escript function 
182 
returning the Kronecker symbol. 
183 

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

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

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

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

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

365 
\begin{figure} 
366 
\centerline{\includegraphics[width=\figwidth]{DiffusionRes1}} 
367 
\centerline{\includegraphics[width=\figwidth]{DiffusionRes16}} 
368 
\centerline{\includegraphics[width=\figwidth]{DiffusionRes32}} 
369 
\centerline{\includegraphics[width=\figwidth]{DiffusionRes48}} 
370 
\caption{Results of the Temperature Diffusion Problem for Time Steps $1$ $16$, $32$ and $48$.} 
371 
\label{DIFFUSION FIG 2} 
372 
\end{figure} 
373 
An easy way to visualize the results is the command 
374 
\begin{verbatim} 
375 
mayavi d T.1.xml m SurfaceMap & 
376 
\end{verbatim} 
377 
Use the \texttt{Configure Data} 
378 
to move forward and and backwards in time. 
379 
\fig{DIFFUSION FIG 2} shows the result for some selected time steps. 