26 
The unknown temperature $T$ is a function of its location in the domain and time $t>0$. The governing equation 
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 
in the interior of the domain is given by 
28 
\begin{equation} 
\begin{equation} 
29 
\rho c\hackscore p T\hackscore{,t}  (\kappa T\hackscore{,i})\hackscore{,i} = q 
\rho c\hackscore p T\hackscore{,t}  (\kappa T\hackscore{,i})\hackscore{,i} = q\hackscore H 
30 
\label{DIFFUSION TEMP EQ 1} 
\label{DIFFUSION TEMP EQ 1} 
31 
\end{equation} 
\end{equation} 
32 
where $\rho c\hackscore p$ and $\kappa$ are given material constants. In case of a composite 
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$ is 
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} 
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$ to be equal to a constant heat production rate 
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: 
$q^{c}$ on a circle or sphere with center $x^c$ and radius $r$ and $0$ elsewhere: 
37 
\begin{equation} 
\begin{equation} 
38 
q(x,t)= 
q\hackscore H(x,t)= 
39 
\left\{ 
\left\{ 
40 
\begin{array}{lcl} 
\begin{array}{lcl} 
41 
q^c & & \xx^c\ \le r \\ 
q^c & & \xx^c\ \le r \\ 
85 
This is inserted into \eqn{DIFFUSION TEMP EQ 1}. By separating the terms at 
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$ 
$t^{(n)}$ and $t^{(n1)}$ one gets for $n=1,2,3\ldots$ 
87 
\begin{equation} 
\begin{equation} 
88 
\frac{\rho c\hackscore p}{h} T^{(n)}  (\kappa T^{(n)}\hackscore{,i})\hackscore{,i} = q + \frac{\rho c\hackscore p}{h} T^{(n1)} 
\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} 
\label{DIFFUSION TEMP EQ 7} 
90 
\end{equation} 
\end{equation} 
91 
where $T^{(0)}=T\hackscore{ref}$ is taken form the initial condition given by \eqn{DIFFUSION TEMP EQ 4}. 
where $T^{(0)}=T\hackscore{ref}$ is taken form the initial condition given by \eqn{DIFFUSION TEMP EQ 4}. 
102 
\subsection{\label{DIFFUSION HELM SEC}Helmholtz Problem} 
\subsection{\label{DIFFUSION HELM SEC}Helmholtz Problem} 
103 
The partial differential equation to be solved for $T^{(n)}$ has the form 
The partial differential equation to be solved for $T^{(n)}$ has the form 
104 
\begin{equation} 
\begin{equation} 
105 
\omega u  (\kappa u\hackscore{,i})\hackscore{,i} = f 
\omega T^{(n)}  (\kappa T^{(n)}\hackscore{,i})\hackscore{,i} = f 
106 
\label{DIFFUSION HELM EQ 1} 
\label{DIFFUSION HELM EQ 1} 
107 
\end{equation} 
\end{equation} 
108 
where $u$ plays the role of $T^{(n)}$ and we set 
and we set 
109 
\begin{equation} 
\begin{equation} 
110 
\omega=\frac{\rho c\hackscore p}{h} \mbox{ and } f=q+\frac{\rho c\hackscore p}{h}T^{(n1)} \;. 
\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} 
\label{DIFFUSION HELM EQ 1b} 
112 
\end{equation} 
\end{equation} 
113 
With $g=\eta T\hackscore{ref}$ the radiation condition defined by \eqn{DIFFUSION TEMP EQ 2222} 
With $g=\eta T\hackscore{ref}$ the radiation condition defined by \eqn{DIFFUSION TEMP EQ 2222} 
114 
takes the form 
takes the form 
115 
\begin{equation} 
\begin{equation} 
116 
\kappa u\hackscore{,i} n\hackscore{i} = g  \eta u\mbox{ on } \Gamma 
\kappa T^{(n)}\hackscore{,i} n\hackscore{i} = g  \eta T^{(n)}\mbox{ on } \Gamma 
117 
\label{DIFFUSION HELM EQ 2} 
\label{DIFFUSION HELM EQ 2} 
118 
\end{equation} 
\end{equation} 
119 
The partial differential \eqn{DIFFUSION HELM EQ 1} together with boundary conditions of \eqn{DIFFUSION HELM EQ 2} 
The partial differential \eqn{DIFFUSION HELM EQ 1} together with boundary conditions of \eqn{DIFFUSION HELM EQ 2} 
148 
we will make direct use of the general \LinearPDE class. 
we will make direct use of the general \LinearPDE class. 
149 


150 
By inspecting the Helmholtz equation \index{Helmholtz equation} 
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 
we can easily assign values to the coefficients in the 
154 
general PDE of the \LinearPDE class: 
general PDE of the \LinearPDE class: 
155 
\begin{equation}\label{DIFFUSION HELM EQ 3} 
\begin{equation}\label{DIFFUSION HELM EQ 3} 
159 
\end{array} 
\end{array} 
160 
\end{equation} 
\end{equation} 
161 
$\delta\hackscore{ij}$ is the Kronecker symbol \index{Kronecker symbol} defined by $\delta\hackscore{ij}=1$ for 
$\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 
$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, 
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.} 
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: 
Defining and solving the Helmholtz equation is very easy now: 
170 
\begin{python} 
\begin{python} 
220 
on a rectangular domain 
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 
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 
can verify the returned solution against the exact answer. Actually, we 
223 
we take $u=x\hackscore{0}$ and then calculate the right hand side terms $f$ and $g$ such that 
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, 
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 
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}$. 
$\kappa n\hackscore{i} u\hackscore{,i}=\kappa n\hackscore{0}$. 
247 
#... calculate error of the PDE solution ... 
#... calculate error of the PDE solution ... 
248 
u=mypde.getSolution() 
u=mypde.getSolution() 
249 
print "error is ",Lsup(ux[0]) 
print "error is ",Lsup(ux[0]) 
250 

saveVTK("x0.xml",sol=u) 
251 
\end{python} 
\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}. 
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} 
\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. 
imported by the \code{from esys.escript import *} statement and returns the maximum absolute value of its argument. 
291 
$Y$ is reset as it depends on the temperature of the previous time step. This allows the PDE solver to reuse information 
$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. 
from previous time steps as much as possible. 
293 


294 
The heat source \var{q} which is defined in \eqn{DIFFUSION TEMP EQ 1b} is \var{qc} 
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. 
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 \var{q} as desired: 
\var{q0} is a fixed constant. The following script defines $q\hackscore H$ as desired: 
297 
\begin{python} 
\begin{python} 
298 
from esys.escript import length,whereNegative 
from esys.escript import length,whereNegative 
299 
xc=[0.02,0.002] 
xc=[0.02,0.002] 
300 
r=0.001 
r=0.001 
301 
x=mydomain.getX() 
x=mydomain.getX() 
302 
q=q0*whereNegative(length(xxc)r) 
qH=q0*whereNegative(length(xxc)r) 
303 
\end{python} 
\end{python} 
304 
\var{x} is a \Data class object of 
\var{x} is a \Data class object of 
305 
the \escript module defining locations in the \Domain \var{mydomain}. 
the \escript module defining locations in the \Domain \var{mydomain}. 
346 
mypde.setValue(A=kappa*kronecker(mydomain),D=rhocp/h,d=eta,y=eta*Tref) 
mypde.setValue(A=kappa*kronecker(mydomain),D=rhocp/h,d=eta,y=eta*Tref) 
347 
# ... set heat source: .... 
# ... set heat source: .... 
348 
x=mydomain.getX() 
x=mydomain.getX() 
349 
q=qc*whereNegative(length(xxc)r) 
qH=qc*whereNegative(length(xxc)r) 
350 
# ... set initial temperature .... 
# ... set initial temperature .... 
351 
T=Tref 
T=Tref 
352 
# ... start iteration: 
# ... start iteration: 
354 
i+=1 
i+=1 
355 
t+=h 
t+=h 
356 
print "time step :",t 
print "time step :",t 
357 
mypde.setValue(Y=q+rhocp/h*T) 
mypde.setValue(Y=qH+rhocp/h*T) 
358 
T=mypde.getSolution() 
T=mypde.getSolution() 
359 
saveVTK("T.%d.xml"%i,temp=T) 
saveVTK("T.%d.xml"%i,temp=T) 
360 
\end{python} 
\end{python} 