1 
jgs 
102 
% $Id$ 
2 
jgs 
107 
\chapter{How to Solve A Diffusion Problem} 
3 
jgs 
102 
\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 


\section{\label{DIFFUSION OUT SEC}Outline} 
12 
jgs 
107 
In this chapter we will discuss how to solve the time dependenttemperature diffusion\index{diffusion equation} for 
13 
jgs 
102 
a block of material. Within the block there is a heat source which drives the temperature diffusion. 
14 
jgs 
107 
On the surface, energy can radiate into the surrounding environment. 
15 
jgs 
102 
\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 
jgs 
107 
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 
jgs 
102 
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 


\section{\label{DIFFUSION TEMP SEC}Temperature Diffusion} 
26 



27 
jgs 
107 
The unknown temperature $T$ is a function of its location in the domain and time $t>0$. The governing equation 
28 
jgs 
102 
in the interior of the domain is given by 
29 


\begin{equation} 
30 


\rho c\hackscore p T\hackscore{,t}  (\kappa T\hackscore{,i})\hackscore{,i} = q 
31 


\label{DIFFUSION TEMP EQ 1} 
32 


\end{equation} 
33 


where $\rho c\hackscore p$ and $\kappa$ are given material constants. In case of a composite 
34 
jgs 
107 
material the parameters depend on their location in the domain. $q$ is 
35 
jgs 
102 
a heat source (or sink) within the domain. We are using Einstein summation convention \index{summation convention} 
36 
jgs 
107 
as introduced in \Chap{FirstSteps}. In our case we assume $q$ to be equal to a constant heat production rate 
37 


$q^{c}$ on a circle or sphere with center $x^c$ and radius $r$ and $0$ elsewhere: 
38 
jgs 
102 
\begin{equation} 
39 


q(x,t)= 
40 


\left\{ 
41 


\begin{array}{lcl} 
42 


q^c & & \xx^c\ \le r \\ 
43 


& \mbox{if} \\ 
44 


0 & & \mbox{else} \\ 
45 


\end{array} 
46 


\right. 
47 


\label{DIFFUSION TEMP EQ 1b} 
48 


\end{equation} 
49 


for all $x$ in the domain and all time $t>0$. 
50 



51 


On the surface of the domain we are 
52 
jgs 
107 
specifying a radiation condition 
53 


which precribes the normal component of the flux $\kappa T\hackscore{,i}$ to be proportional 
54 
jgs 
102 
to the difference of the current temperature to the surrounding temperature $T\hackscore{ref}$: 
55 


\begin{equation} 
56 


\kappa T\hackscore{,i} n\hackscore i = \eta (T\hackscore{ref}T) 
57 


\label{DIFFUSION TEMP EQ 2} 
58 


\end{equation} 
59 
jgs 
107 
$\eta$ is a given material coefficient depending on the material of the block and the surrounding medium. 
60 


As usual $n\hackscore i$ is the $i$th component of the outer normal field \index{outer normal field} 
61 
jgs 
102 
at the surface of the domain. 
62 



63 
jgs 
107 
To solve the time dependent \eqn{DIFFUSION TEMP EQ 1} the initial temperature at time 
64 
jgs 
102 
$t=0$ has to be given. Here we assume that the initial temperature is the surrounding temperature: 
65 


\begin{equation} 
66 


T(x,0)=T\hackscore{ref} 
67 


\label{DIFFUSION TEMP EQ 4} 
68 


\end{equation} 
69 


for all $x$ in the domain. It is pointed out that 
70 
jgs 
107 
the initial conditions satisfy the 
71 
jgs 
102 
boundary condition defined by \eqn{DIFFUSION TEMP EQ 2}. 
72 



73 
jgs 
107 
The temperature is calculated at discrete time nodes $t^{(n)}$ where 
74 
jgs 
102 
$t^{(0)}=0$ and $t^{(n)}=t^{(n1)}+h$ where $h>0$ is the step size which is assumed to be constant. 
75 
jgs 
107 
In the following the upper index ${(n)}$ refers to a value at time $t^{(n)}$. The simplest 
76 


and most robust scheme to approximate the time derivative of the the temperature is 
77 


the backward Euler 
78 


\index{backward Euler} scheme, see~\cite{XXX} for alternatives. The backward Euler 
79 


scheme is based 
80 


on the Taylor expansion of $T$ at time $t^{(n)}$: 
81 
jgs 
102 
\begin{equation} 
82 


T^{(n1)}\approx T^{(n)}+T\hackscore{,t}^{(n)}(t^{(n1)}t^{(n)}) 
83 
jgs 
107 
=T^{(n1)}  h \cdot T\hackscore{,t}^{(n)} 
84 
jgs 
102 
\label{DIFFUSION TEMP EQ 6} 
85 


\end{equation} 
86 
jgs 
107 
This is inserted into \eqn{DIFFUSION TEMP EQ 1}. By separating the terms at 
87 
jgs 
102 
$t^{(n)}$ and $t^{(n1)}$ one gets for $n=1,2,3\ldots$ 
88 


\begin{equation} 
89 


\frac{\rho c\hackscore p}{h} T^{(n)}  (\kappa T^{(n)}\hackscore{,i})\hackscore{,i} = q + \frac{\rho c\hackscore p}{h} T^{(n1)} 
90 


\label{DIFFUSION TEMP EQ 7} 
91 


\end{equation} 
92 


where $T^{(0)}=T\hackscore{ref}$ is taken form the initial condition given by \eqn{DIFFUSION TEMP EQ 4}. 
93 
jgs 
107 
Together with the natural boundary condition 
94 


\begin{equation} 
95 


\kappa T\hackscore{,i}^{(n)} n\hackscore i = \eta (T\hackscore{ref}T^{(n)}) 
96 


\label{DIFFUSION TEMP EQ 2222} 
97 


\end{equation} 
98 


taken from \eqn{DIFFUSION TEMP EQ 2} 
99 
jgs 
102 
this forms a boundary value problem that has to be solved for each time step. 
100 


As a first step to implement a solver for the temperature diffusion problem we will 
101 


first implement a solver for the boundary value problem that has to be solved at each time step. 
102 



103 


\section{\label{DIFFUSION HELM SEC}Helmholtz Problem} 
104 


The partial differential equation to be solved for $T^{(n)}$ has the form 
105 


\begin{equation} 
106 


\omega u  (\kappa u\hackscore{,i})\hackscore{,i} = f 
107 


\label{DIFFUSION HELM EQ 1} 
108 


\end{equation} 
109 


where $u$ plays the role of $T^{(n)}$ and we set 
110 


\begin{equation} 
111 


\omega=\frac{\rho c\hackscore p}{h} \mbox{ and } f=q+\frac{\rho c\hackscore p}{h}T^{(n1)} \;. 
112 


\label{DIFFUSION HELM EQ 1b} 
113 


\end{equation} 
114 
jgs 
107 
With $g=\eta T\hackscore{ref}$ the radiation condition defined by \eqn{DIFFUSION TEMP EQ 2222} 
115 
jgs 
102 
takes the form 
116 


\begin{equation} 
117 


\kappa u\hackscore{,i} n\hackscore{i} = g  \eta u\mbox{ on } \Gamma 
118 


\label{DIFFUSION HELM EQ 2} 
119 


\end{equation} 
120 


The partial differential 
121 


\eqn{DIFFUSION HELM EQ 1} together with boundary conditions of \eqn{DIFFUSION HELM EQ 2} 
122 
jgs 
107 
is called the Helmholtz equation \index{Helmholtz equation}. 
123 
jgs 
102 

124 


We want to use the \LinearPDE class provided by \escript to define and solve a general linear PDE such as the 
125 


Helmholtz equation. We have used a special case of the \LinearPDE class, namely the 
126 


\Poisson class already in \Chap{FirstSteps}. 
127 
jgs 
107 
Here we will write our own specialized subclass of the \LinearPDE to define the Helmholtz equation 
128 


and use the \method{getSolution} method of parent class to actually solve the problem. 
129 
jgs 
102 

130 
jgs 
107 
The form of a single PDE that can be handled by the \LinearPDE class is 
131 
jgs 
102 
\begin{equation}\label{EQU.FEM.1} 
132 


(A\hackscore{jl} u\hackscore{,l})\hackscore{,j}+D u =Y \; . 
133 


\end{equation} 
134 
jgs 
107 
We show here the terms which are relevant for the Helmholtz problem. 
135 
jgs 
102 
The general form and systems is discussed in \Sec{SEC LinearPDE}. 
136 
jgs 
107 
$A$, $D$ and $Y$ are the known coeffecients of the PDE. \index{partial differential equation!coefficients} 
137 
jgs 
102 
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 
jgs 
107 
By inspecting the Helmholtz equation \index{Helmholtz equation} 
149 
jgs 
102 
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 
jgs 
107 
$\delta\hackscore{ij}$ is the Kronecker symbol \index{Kronecker symbol} defined by $\delta\hackscore{ij}=1$ for 
158 
jgs 
102 
$i=j$ and $0$ otherwise. 
159 



160 


We want to implement a 
161 
jgs 
107 
new class which we will call \class{Helmholtz} that provides the same methods as the \LinearPDE class but 
162 


is described using the coefficients $\kappa$, $\omega$, $f$, $\eta$, 
163 
jgs 
102 
$g$ rather than the general form given by \eqn{EQU.FEM.1}. 
164 
jgs 
107 
Python's mechanism of subclasses allows us to do this in a very simple way. 
165 
jgs 
102 
The \Poisson class of the \linearPDEsPack module, 
166 
jgs 
107 
which we have already used in \Chap{FirstSteps}, is in fact a subclass of the more general 
167 


\LinearPDE class. That means that all methods (such as the \method{getSolution}) 
168 


from the parent class \LinearPDE are available for any \Poisson object. However, new 
169 
jgs 
102 
methods can be added and methods of the parent class can be redefined. In fact, 
170 


the \Poisson class redefines the \method{setValue} of the \LinearPDE class which is called to assign 
171 


values to the coefficients of the PDE. This is exactly what we will do when we define 
172 


our new \class{Helmholtz} class: 
173 


\begin{python} 
174 


from esys.linearPDEs import LinearPDE 
175 


import numarray 
176 
jgs 
107 
class Helmholtz(LinearPDE): 
177 
jgs 
102 
def setValue(self,kappa=0,omega=1,f=0,eta=0,g=0) 
178 
jgs 
107 
ndim=self.getDim() 
179 


kronecker=numarray.identity(ndim) 
180 


self._setValue(A=kappa*kronecker,D=omega,Y=f,d=eta,y=g) 
181 
jgs 
102 
\end{python} 
182 


\code{class Helmholtz(linearPDE)} declares the new \class{Helmholtz} class as a subclass 
183 
jgs 
107 
of \LinearPDE which we have imported in the first line of the script. 
184 
jgs 
102 
We add the method \method{setValue} to the class which overwrites the 
185 
jgs 
107 
\method{setValue} method of the \LinearPDE class. The new method which has the 
186 
jgs 
102 
parameters of the Helmholtz \eqn{DIFFUSION HELM EQ 1} as arguments 
187 


maps the parameters of the coefficients of the general PDE defined 
188 
jgs 
107 
in \eqn{EQU.FEM.1}. We are actually using the \method{_setValue} of 
189 


the \LinearPDE class (notice the leeding underscoure). 
190 


The coefficient \var{A} involves the Kronecker symbol. We use the 
191 


\numarray function \function{identity} which returns a square matrix with ones on the 
192 


main diagonal and zeros off the main diagonal. The argument of \function{identity} gives the order of the matrix. 
193 


The \method{getDim} of the \LinearPDE class object \var{self} to get the spatial dimensions of the domain of the 
194 


PDE. As we will make use of the \class{Helmholtz} class several times, it is convenient to 
195 


put its definition into a file which we name \file{mytools.py} available in the \ExampleDirectory. 
196 
jgs 
102 
You can use your favourite editor to create and edit the file. 
197 



198 


An object of the \class{Helmholtz} class is created through the statments: 
199 


\begin{python} 
200 
jgs 
107 
from mytools import Helmholtz 
201 


mypde = Helmholtz(mydomain) 
202 
jgs 
102 
mypde.setValue(kappa=10.,omega=0.1,f=12) 
203 
jgs 
107 
u = mypde.getSolution() 
204 
jgs 
102 
\end{python} 
205 
jgs 
107 
In the first statement we import all definition from the \file{mytools.py} \index{scripts!\file{mytools.py}}. Make sure 
206 


that \file{mytools.py} is in the directory from where you started Python. 
207 


\var{mydomain} is the \Domain of the PDE which we have undefined here. In the third statment values are 
208 
jgs 
102 
assigned to the PDE parameters. As no values for arguments \var{eta} and \var{g} are 
209 
jgs 
107 
specified, the default values $0$ are used. \footnote{It would be better to use the default value 
210 
jgs 
102 
\var{escript.Data()} rather then $0$ as then the coefficient would be defined as being not present and 
211 
jgs 
107 
would not be processed when the PDE is evaluated}. In the fourth statement the solution of the 
212 
jgs 
102 
PDE is returned. 
213 



214 
jgs 
107 
To test our \class{Helmholtz} class on a rectangular domain 
215 


of length $l\hackscore{0}=5$ and height $l\hackscore{1}=1$, we choose a simple test solution. Actually, we 
216 


we take $u=x\hackscore{0}$ and then calculate the right hand side terms $f$ and $g$ such that 
217 


the test solution becomes the solution of the problem. If we assume $\kappa$ as being constant, 
218 
jgs 
102 
an easy calculation shows that we have to choose $f=\omega \cdot x\hackscore{0}$. On the boundary we get 
219 
jgs 
107 
$\kappa n\hackscore{i} u\hackscore{,i}=\kappa n\hackscore{0}$. 
220 


So we have to set $g=\kappa n\hackscore{0}+\eta x\hackscore{0}$. The following script \file{helmholtztest.py} 
221 
jgs 
102 
\index{scripts!\file{helmholtztest.py}} which is available in the \ExampleDirectory 
222 


implements this test problem using the \finley PDE solver: 
223 


\begin{python} 
224 
jgs 
107 
from mytools import Helmholtz 
225 


from esys.escript import Lsup 
226 


from esys.finley import Rectangle 
227 
jgs 
102 
#... set some parameters ... 
228 
jgs 
107 
kappa=1. 
229 
jgs 
102 
omega=0.1 
230 


eta=10. 
231 


#... generate domain ... 
232 
jgs 
107 
mydomain = Rectangle(l0=5.,l1=1.,n0=50, n1=10) 
233 
jgs 
102 
#... open PDE and set coefficients ... 
234 


mypde=Helmholtz(mydomain) 
235 


n=mydomain.getNormal() 
236 


x=mydomain.getX() 
237 
jgs 
107 
mypde.setValue(kappa,omega,omega*x[0],eta,kappa*n[0]+eta*x[0]) 
238 
jgs 
102 
#... calculate error of the PDE solution ... 
239 


u=mypde.getSolution() 
240 


print "error is ",Lsup(ux[0]) 
241 


\end{python} 
242 


The script is similar to the script \file{mypoisson.py} dicussed in \Chap{FirstSteps}. 
243 


\code{mydomain.getNormal()} returns the outer normal field on the surface of the domain. The function \function{Lsup} 
244 
jgs 
107 
is imported by the \code{from esys.escript import Lsup} statement and returns the maximum absulute value of its argument. 
245 


The error shown by the print statement should be in the order of $10^{7}$. As piecewise bilinear interpolation is 
246 


used to approximate the solution and our solution is a linear function of the spatial coordinates one might 
247 


expect that the error would be zero or in the order of machine precision (typically $\approx 10^{15}$). 
248 


However most PDE packages use an iterative solver which is terminated 
249 


when a given tolerance has been reached. The default tolerance is $10^{8}$. This value can be altered by using the 
250 
jgs 
102 
\method{setTolerance} of the \LinearPDE class. 
251 



252 


\section{The Transition Problem} 
253 


\label{DIFFUSION TRANS SEC} 
254 


Now we are ready to solve the original time dependent problem. The main 
255 
jgs 
107 
part of the script is the loop over time $t$ which takes the following form: 
256 
jgs 
102 
\begin{python} 
257 
jgs 
107 
t=0 
258 


T=Tref 
259 
jgs 
102 
mypde=Helmholtz(mydomain) 
260 


while t<t_end: 
261 


mypde.setValue(kappa,rhocp/h,q+rhocp/h*T,eta,eta*Tref) 
262 


T=mypde.getSolution() 
263 


t+=h 
264 


\end{python} 
265 


\var{kappa}, \var{rhocp}, \var{eta} and \var{Tref} are input parameters of the model. \var{q} is the heat source 
266 
jgs 
107 
in the domain and \var{h} is the time step size. Notice that the \class{Hemholtz} class object \var{mypde} 
267 
jgs 
102 
is created before the loop over time is entered while in each time step only the coefficients 
268 
jgs 
107 
are reset in each time step. This way some information about the representation of the PDE can be reused 
269 
jgs 
102 
\footnote{The efficience can be improved further by setting the coefficients in the operator 
270 
jgs 
107 
\var{kappa}, \var{omega} and \var{eta} before entering the \code{while}loop and only updating the coefficients 
271 
jgs 
102 
in the right hand side \var{f} and \var{g}. This needs a more careful implementation of the \method{setValue} 
272 
jgs 
107 
method but gives the advantage that the \LinearPDE class can save rebuilding the PDE operator.}. The variable \var{T} 
273 
jgs 
102 
holds the current temperature. It is used to calculate the right hand side coefficient \var{f} in the 
274 
jgs 
107 
Helmholtz equation in \eqn{DIFFUSION HELM EQ 1}. The statement \code{T=mypde.getSolution()} overwrites \var{T} with the 
275 


temperature of the new time step $\var{t}+\var{h}$. To get this iterative process going we need to specify the 
276 
jgs 
102 
initial temperature distribution, which equal to $T\hackscore{ref}$. 
277 



278 
jgs 
107 
The heat source \var{q} which is defined in \eqn{DIFFUSION TEMP EQ 1b} is \var{qc} 
279 


in an area defined as a circle of radius \var{r} and center \var{xc} and zero outside this circle. 
280 
jgs 
102 
\var{q0} is a fixed constant. The following script defines \var{q} as desired: 
281 


\begin{python} 
282 
jgs 
107 
from esys.escript import length 
283 
jgs 
102 
xc=[0.02,0.002] 
284 


r=0.001 
285 


x=mydomain.getX() 
286 


q=q0*(length(xxc)r).whereNegative() 
287 


\end{python} 
288 


\var{x} is a \Data class object of 
289 
jgs 
107 
the \escript module defining locations in the \Domain \var{mydomain}. 
290 


The \function{length()} imported from the \escript module returns the 
291 


Euclidean norm: 
292 


\begin{equation}\label{DIFFUSION HELM EQ 3aba} 
293 


\y\=\sqrt{ 
294 


y\hackscore{i} 
295 


y\hackscore{i} 
296 


} = \function{esys.escript.length}(y) 
297 


\end{equation} 
298 


So \code{length(xxc)} calculates the distances 
299 


of the location \var{x} to the center of the circle \var{xc} where the heat source is acting. 
300 


Note that the coordinates of \var{xc} are defined as a list of floating point numbers. It is independently 
301 


converted into a \Data class object before being subtracted from \var{x}. The method \method{whereNegative} of 
302 
jgs 
102 
a \Data class object, in this case the result of the expression 
303 
jgs 
107 
\code{length(xxc)r}, returns a \Data class which is equal to one where the object is negative and 
304 


zero elsewhere. After multiplication with \var{qc} we get a function with the desired property. 
305 
jgs 
102 

306 
jgs 
107 
Now we can put the components together to create the script \file{diffusion.py} which is available in the \ExampleDirectory: 
307 
jgs 
102 
\index{scripts!\file{diffusion.py}}: 
308 


\begin{python} 
309 
jgs 
107 
from mytools import Helmholtz 
310 


from esys.escript import Lsup 
311 


from esys.finley import Rectangle 
312 
jgs 
102 
#... set some parameters ... 
313 
jgs 
107 
xc=[0.02,0.002] 
314 
jgs 
102 
r=0.001 
315 
jgs 
107 
qc=50.e6 
316 
jgs 
102 
Tref=0. 
317 


rhocp=2.6e6 
318 


eta=75. 
319 


kappa=240. 
320 
jgs 
107 
tend=5. 
321 


# ... time, time step size and counter ... 
322 


t=0 
323 
jgs 
102 
h=0.1 
324 


i=0 
325 


#... generate domain ... 
326 
jgs 
107 
mydomain = Rectangle(l0=0.05,l1=0.01,n0=250, n1=50) 
327 
jgs 
102 
#... open PDE ... 
328 


mypde=Helmholtz(mydomain) 
329 


# ... set heat source: .... 
330 


x=mydomain.getX() 
331 
jgs 
107 
q=qc*(length(xxc)r).whereNegative() 
332 
jgs 
102 
# ... set initial temperature .... 
333 


T=Tref 
334 


# ... start iteration: 
335 
jgs 
107 
while t<tend: 
336 
jgs 
102 
i+=1 
337 


t+=h 
338 


print "time step :",t 
339 


mypde.setValue(kappa=kappa,omega=rhocp/h,f=q+rhocp/h*T,eta=eta,g=eta*Tref) 
340 


T=mypde.getSolution() 
341 


T.saveDX("T%d.dx"%i) 
342 


\end{python} 
343 


The script will create the files \file{T.1.dx}, 
344 


\file{T.2.dx}, $\ldots$, \file{T.50.dx} in the directory where the script has been started. The files give the 
345 


temperature distributions at time steps $1$, $2$, $\ldots$, $50$ in the \OpenDX file format. 
346 
jgs 
107 

347 


\begin{figure} 
348 


\centerline{\includegraphics[width=\figwidth]{DiffusionRes1}} 
349 


\centerline{\includegraphics[width=\figwidth]{DiffusionRes16}} 
350 


\centerline{\includegraphics[width=\figwidth]{DiffusionRes32}} 
351 


\centerline{\includegraphics[width=\figwidth]{DiffusionRes48}} 
352 


\caption{Results of the Temperture Diffusion Problem for Time Steps $1$ $16$, $32$ and $48$.} 
353 


\label{DIFFUSION FIG 2} 
354 


\end{figure} 
355 



356 
jgs 
102 
An easy way to visualize the results is the command 
357 


\begin{verbatim} 
358 
jgs 
107 
dx edit diffusion.net & 
359 
jgs 
102 
\end{verbatim} 
360 
jgs 
107 
where \file{diffusion.net} is an \OpenDX script available in the \ExampleDirectory. 
361 


Use the \texttt{Sequencer} to move forward and and backwards in time. 
362 


\fig{DIFFUSION FIG 2} shows the result for some selected time steps. 