1 
jgs 
102 
% $Id$ 
2 
gross 
625 
% 
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 
jgs 
121 
\section{The Diffusion Problem} 
11 
jgs 
102 
\label{DIFFUSION CHAP} 
12 



13 


\begin{figure} 
14 
gross 
599 
\centerline{\includegraphics[width=\figwidth]{figures/DiffusionDomain.eps}} 
15 
jgs 
102 
\caption{Temperature Diffusion Problem with Circular Heat Source} 
16 


\label{DIFFUSION FIG 1} 
17 


\end{figure} 
18 



19 
jgs 
121 
\subsection{\label{DIFFUSION OUT SEC}Outline} 
20 
jgs 
107 
In this chapter we will discuss how to solve the time dependenttemperature diffusion\index{diffusion equation} for 
21 
jgs 
102 
a block of material. Within the block there is a heat source which drives the temperature diffusion. 
22 
jgs 
107 
On the surface, energy can radiate into the surrounding environment. 
23 
jgs 
102 
\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 
jgs 
107 
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 
jgs 
102 
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 
jgs 
121 
\subsection{\label{DIFFUSION TEMP SEC}Temperature Diffusion} 
34 
jgs 
107 
The unknown temperature $T$ is a function of its location in the domain and time $t>0$. The governing equation 
35 
jgs 
102 
in the interior of the domain is given by 
36 


\begin{equation} 
37 
lkettle 
575 
\rho c\hackscore p T\hackscore{,t}  (\kappa T\hackscore{,i})\hackscore{,i} = q\hackscore H 
38 
jgs 
102 
\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 
lkettle 
575 
material the parameters depend on their location in the domain. $q\hackscore H$ is 
42 
jgs 
102 
a heat source (or sink) within the domain. We are using Einstein summation convention \index{summation convention} 
43 
lkettle 
575 
as introduced in \Chap{FirstSteps}. In our case we assume $q\hackscore H$ to be equal to a constant heat production rate 
44 
jgs 
107 
$q^{c}$ on a circle or sphere with center $x^c$ and radius $r$ and $0$ elsewhere: 
45 
jgs 
102 
\begin{equation} 
46 
lkettle 
575 
q\hackscore H(x,t)= 
47 
jgs 
102 
\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 
jgs 
107 
specifying a radiation condition 
60 


which precribes the normal component of the flux $\kappa T\hackscore{,i}$ to be proportional 
61 
jgs 
102 
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 
jgs 
107 
$\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 
jgs 
102 
at the surface of the domain. 
69 



70 
jgs 
107 
To solve the time dependent \eqn{DIFFUSION TEMP EQ 1} the initial temperature at time 
71 
jgs 
102 
$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 
jgs 
107 
the initial conditions satisfy the 
78 
jgs 
102 
boundary condition defined by \eqn{DIFFUSION TEMP EQ 2}. 
79 



80 
jgs 
107 
The temperature is calculated at discrete time nodes $t^{(n)}$ where 
81 
jgs 
102 
$t^{(0)}=0$ and $t^{(n)}=t^{(n1)}+h$ where $h>0$ is the step size which is assumed to be constant. 
82 
jgs 
107 
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 
gross 
660 
\index{backward Euler} scheme. The backward Euler 
86 
jgs 
107 
scheme is based 
87 


on the Taylor expansion of $T$ at time $t^{(n)}$: 
88 
jgs 
102 
\begin{equation} 
89 
lkettle 
573 
T^{(n)}\approx T^{(n1)}+T\hackscore{,t}^{(n)}(t^{(n)}t^{(n1)}) 
90 


=T^{(n1)} + h \cdot T\hackscore{,t}^{(n)} 
91 
jgs 
102 
\label{DIFFUSION TEMP EQ 6} 
92 


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


\begin{equation} 
96 
lkettle 
575 
\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 
jgs 
102 
\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 
jgs 
107 
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 
jgs 
102 
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 
jgs 
121 
\subsection{\label{DIFFUSION HELM SEC}Helmholtz Problem} 
111 
jgs 
102 
The partial differential equation to be solved for $T^{(n)}$ has the form 
112 


\begin{equation} 
113 
lkettle 
575 
\omega T^{(n)}  (\kappa T^{(n)}\hackscore{,i})\hackscore{,i} = f 
114 
jgs 
102 
\label{DIFFUSION HELM EQ 1} 
115 


\end{equation} 
116 
lkettle 
575 
and we set 
117 
jgs 
102 
\begin{equation} 
118 
lkettle 
575 
\omega=\frac{\rho c\hackscore p}{h} \mbox{ and } f=q\hackscore H +\frac{\rho c\hackscore p}{h}T^{(n1)} \;. 
119 
jgs 
102 
\label{DIFFUSION HELM EQ 1b} 
120 


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


\begin{equation} 
124 
lkettle 
575 
\kappa T^{(n)}\hackscore{,i} n\hackscore{i} = g  \eta T^{(n)}\mbox{ on } \Gamma 
125 
jgs 
102 
\label{DIFFUSION HELM EQ 2} 
126 


\end{equation} 
127 
gross 
568 
The partial differential \eqn{DIFFUSION HELM EQ 1} together with boundary conditions of \eqn{DIFFUSION HELM EQ 2} 
128 
jgs 
107 
is called the Helmholtz equation \index{Helmholtz equation}. 
129 
jgs 
102 

130 
gross 
568 
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 
gross 
625 
\begin{equation}\label{LINEARPDE.SINGLE.1 TUTORIAL} 
133 


(A\hackscore{jl} u\hackscore{,l})\hackscore{,j}+D u = Y \; . 
134 
jgs 
102 
\end{equation} 
135 
gross 
625 
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 
gross 
568 
\Function on the PDE or objects that can be converted into such \Data objects. 
139 
gross 
625 
$A$ is a \RankTwo and $D$ and $Y$ are scalar. 
140 
gross 
568 
The following natural 
141 


boundary conditions are considered \index{boundary condition!natural} on $\Gamma$: 
142 
gross 
625 
\begin{equation}\label{LINEARPDE.SINGLE.2 TUTORIAL} 
143 


n\hackscore{j}A\hackscore{jl} u\hackscore{,l}+d u= y \;. 
144 
jgs 
102 
\end{equation} 
145 
gross 
625 
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 
gross 
568 
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 
gross 
625 
\begin{equation}\label{LINEARPDE.SINGLE.3 TUTORIAL} 
150 
gross 
568 
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 
gross 
625 
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 
gross 
568 
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 
jgs 
102 

161 
jgs 
107 
By inspecting the Helmholtz equation \index{Helmholtz equation} 
162 
lkettle 
575 
(\ref{DIFFUSION HELM EQ 1}) and boundary condition (\ref{DIFFUSION HELM EQ 2}) and 
163 


substituting $u$ for $T^{(n)}$ 
164 
jgs 
102 
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 
jgs 
107 
$\delta\hackscore{ij}$ is the Kronecker symbol \index{Kronecker symbol} defined by $\delta\hackscore{ij}=1$ for 
173 
lkettle 
575 
$i=j$ and $0$ otherwise. Undefined coefficients are assumed to be not present.\footnote{There is a difference 
174 
gross 
568 
in \escript of being not present and set to zero. As not present coefficients are not processed, 
175 
lkettle 
575 
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 
gross 
625 
conditions which are already defined in the \LinearPDE class (shown in \eqn{LINEARPDE.SINGLE.2 TUTORIAL}). 
179 
jgs 
102 

180 
gross 
568 
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 
gross 
569 
the PDE and returns the solution \code{u} of the PDE. \function{kronecker} is \escript function 
193 
gross 
568 
returning the Kronecker symbol. 
194 



195 


The coefficients can set by several calls of \method{setValue} where the order can be chosen arbitrarily. 
196 
lkettle 
573 
If a value is assigned to a coefficient several times, the last assigned value is used when 
197 
gross 
568 
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 
gross 
625 
\begin{equation}\label{LINEARPDE.SINGLE.4 TUTORIAL} 
208 


A\hackscore{jl}=A\hackscore{lj}\;. 
209 
gross 
568 
\end{equation} 
210 
gross 
625 
Note that $D$ and $d$ may have any value and the right hand sides $Y$, $y$ as well as the constraints 
211 
gross 
568 
are not relevant. The Helmholtz problem is symmetric. 
212 
gross 
569 
The \LinearPDE class provides the method \method{checkSymmetry} method to check if the given PDE is symmetric. 
213 
gross 
568 
\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 
gross 
569 
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 
gross 
568 
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 
gross 
569 
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 
lkettle 
575 
take $T=x\hackscore{0}$ (here $q\hackscore H = 0$) and then calculate the right hand side terms $f$ and $g$ such that 
235 
jgs 
107 
the test solution becomes the solution of the problem. If we assume $\kappa$ as being constant, 
236 
jgs 
102 
an easy calculation shows that we have to choose $f=\omega \cdot x\hackscore{0}$. On the boundary we get 
237 
jgs 
107 
$\kappa n\hackscore{i} u\hackscore{,i}=\kappa n\hackscore{0}$. 
238 
gross 
569 
So we have to set $g=\kappa n\hackscore{0}+\eta x\hackscore{0}$. The following script \file{helmholtz.py} 
239 
gross 
568 
\index{scripts!\file{helmholtz.py}} which is available in the \ExampleDirectory 
240 
jgs 
102 
implements this test problem using the \finley PDE solver: 
241 


\begin{python} 
242 
gross 
568 
from esys.escript import * 
243 
gross 
569 
from esys.escript.linearPDEs import LinearPDE 
244 
jgs 
107 
from esys.finley import Rectangle 
245 
jgs 
102 
#... set some parameters ... 
246 
jgs 
107 
kappa=1. 
247 
jgs 
102 
omega=0.1 
248 


eta=10. 
249 


#... generate domain ... 
250 
jgs 
107 
mydomain = Rectangle(l0=5.,l1=1.,n0=50, n1=10) 
251 
jgs 
102 
#... open PDE and set coefficients ... 
252 
gross 
568 
mypde=LinearPDE(mydomain) 
253 


mypde.setSymmetryOn() 
254 
jgs 
102 
n=mydomain.getNormal() 
255 


x=mydomain.getX() 
256 
gross 
569 
mypde.setValue(A=kappa*kronecker(mydomain),D=omega,Y=omega*x[0], \ 
257 


d=eta,y=kappa*n[0]+eta*x[0]) 
258 
jgs 
102 
#... calculate error of the PDE solution ... 
259 


u=mypde.getSolution() 
260 


print "error is ",Lsup(ux[0]) 
261 
lkettle 
575 
saveVTK("x0.xml",sol=u) 
262 
jgs 
102 
\end{python} 
263 
lkettle 
575 
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 
gross 
569 
The script is similar to the script \file{poisson.py} dicussed in \Chap{FirstSteps}. 
270 
jgs 
102 
\code{mydomain.getNormal()} returns the outer normal field on the surface of the domain. The function \function{Lsup} 
271 
lkettle 
573 
imported by the \code{from esys.escript import *} statement and returns the maximum absolute value of its argument. 
272 
jgs 
107 
The error shown by the print statement should be in the order of $10^{7}$. As piecewise bilinear interpolation is 
273 
gross 
569 
used by \finley approximate the solution and our solution is a linear function of the spatial coordinates one might 
274 
jgs 
107 
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 
jgs 
102 
\method{setTolerance} of the \LinearPDE class. 
278 



279 
jgs 
121 
\subsection{The Transition Problem} 
280 
jgs 
102 
\label{DIFFUSION TRANS SEC} 
281 


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


T=Tref 
286 
gross 
569 
mypde=LinearPDE(mydomain) 
287 


mypde.setValue(A=kappa*kronecker(mydomain),D=rhocp/h,d=eta,y=eta*Tref) 
288 
jgs 
102 
while t<t_end: 
289 
gross 
569 
mypde.setValue(Y=q+rhocp/h*T) 
290 
jgs 
102 
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 
gross 
569 
in the domain and \var{h} is the time step size. 
295 


The variable \var{T} 
296 
jgs 
102 
holds the current temperature. It is used to calculate the right hand side coefficient \var{f} in the 
297 
jgs 
107 
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 
jgs 
102 
initial temperature distribution, which equal to $T\hackscore{ref}$. 
300 
gross 
569 
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 
jgs 
102 

305 
lkettle 
575 
The heat source $q\hackscore H$ which is defined in \eqn{DIFFUSION TEMP EQ 1b} is \var{qc} 
306 
jgs 
107 
in an area defined as a circle of radius \var{r} and center \var{xc} and zero outside this circle. 
307 
lkettle 
575 
\var{q0} is a fixed constant. The following script defines $q\hackscore H$ as desired: 
308 
jgs 
102 
\begin{python} 
309 
lkettle 
573 
from esys.escript import length,whereNegative 
310 
jgs 
102 
xc=[0.02,0.002] 
311 


r=0.001 
312 


x=mydomain.getX() 
313 
lkettle 
575 
qH=q0*whereNegative(length(xxc)r) 
314 
jgs 
102 
\end{python} 
315 


\var{x} is a \Data class object of 
316 
jgs 
107 
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 
gross 
569 
converted into a \Data class object before being subtracted from \var{x}. The function \function{whereNegative} 
329 


applied to 
330 
jgs 
107 
\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 
jgs 
102 

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


\begin{python} 
336 
gross 
569 
from esys.escript import * 
337 


from esys.escript.linearPDEs import LinearPDE 
338 
jgs 
107 
from esys.finley import Rectangle 
339 
jgs 
102 
#... set some parameters ... 
340 
jgs 
107 
xc=[0.02,0.002] 
341 
jgs 
102 
r=0.001 
342 
jgs 
107 
qc=50.e6 
343 
jgs 
102 
Tref=0. 
344 


rhocp=2.6e6 
345 


eta=75. 
346 


kappa=240. 
347 
jgs 
107 
tend=5. 
348 


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


t=0 
350 
jgs 
102 
h=0.1 
351 


i=0 
352 


#... generate domain ... 
353 
jgs 
107 
mydomain = Rectangle(l0=0.05,l1=0.01,n0=250, n1=50) 
354 
jgs 
102 
#... open PDE ... 
355 
gross 
569 
mypde=LinearPDE(mydomain) 
356 


mypde.setSymmetryOn() 
357 


mypde.setValue(A=kappa*kronecker(mydomain),D=rhocp/h,d=eta,y=eta*Tref) 
358 
jgs 
102 
# ... set heat source: .... 
359 


x=mydomain.getX() 
360 
lkettle 
575 
qH=qc*whereNegative(length(xxc)r) 
361 
jgs 
102 
# ... set initial temperature .... 
362 


T=Tref 
363 


# ... start iteration: 
364 
jgs 
107 
while t<tend: 
365 
jgs 
102 
i+=1 
366 


t+=h 
367 


print "time step :",t 
368 
lkettle 
575 
mypde.setValue(Y=qH+rhocp/h*T) 
369 
jgs 
102 
T=mypde.getSolution() 
370 
gross 
569 
saveVTK("T.%d.xml"%i,temp=T) 
371 
jgs 
102 
\end{python} 
372 
gross 
569 
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 
jgs 
107 

376 


\begin{figure} 
377 
gross 
599 
\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 
lkettle 
581 
\caption{Results of the Temperature Diffusion Problem for Time Steps $1$ $16$, $32$ and $48$.} 
382 
jgs 
107 
\label{DIFFUSION FIG 2} 
383 


\end{figure} 
384 
jgs 
102 
An easy way to visualize the results is the command 
385 


\begin{verbatim} 
386 
gross 
569 
mayavi d T.1.xml m SurfaceMap & 
387 
jgs 
102 
\end{verbatim} 
388 
gross 
569 
Use the \texttt{Configure Data} 
389 


to move forward and and backwards in time. 
390 
jgs 
121 
\fig{DIFFUSION FIG 2} shows the result for some selected time steps. 