/[escript]/trunk/doc/user/diffusion.tex
ViewVC logotype

Contents of /trunk/doc/user/diffusion.tex

Parent Directory Parent Directory | Revision Log Revision Log


Revision 2881 - (show annotations)
Thu Jan 28 02:03:15 2010 UTC (9 years, 4 months ago) by jfenwick
File MIME type: application/x-tex
File size: 18907 byte(s)
Don't panic.
Updating copyright stamps

1
2 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
3 %
4 % Copyright (c) 2003-2010 by University of Queensland
5 % Earth Systems Science Computational Center (ESSCC)
6 % http://www.uq.edu.au/esscc
7 %
8 % Primary Business: Queensland, Australia
9 % Licensed under the Open Software License version 3.0
10 % http://www.opensource.org/licenses/osl-3.0.php
11 %
12 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
13
14
15 \section{The Diffusion Problem}
16 \label{DIFFUSION CHAP}
17
18 \begin{figure}
19 \centerline{\includegraphics[width=\figwidth]{figures/DiffusionDomain}}
20 \caption{Temperature Diffusion Problem with Circular Heat Source}
21 \label{DIFFUSION FIG 1}
22 \end{figure}
23
24 \subsection{\label{DIFFUSION OUT SEC}Outline}
25 In this chapter we will discuss how to solve a time-dependent temperature diffusion\index{diffusion equation} PDE for
26 a given block of material. Within the block there is a heat source which drives the temperature diffusion.
27 On the surface, energy can radiate into the surrounding environment.
28 \fig{DIFFUSION FIG 1} shows the configuration.
29
30 In the next \Sec{DIFFUSION TEMP SEC} we will present the relevant model. A
31 time integration scheme is introduced to calculate the temperature at given time nodes $t^{(n)}$.
32 We will see that at each time step a Helmholtz equation \index{Helmholtz equation}
33 must be solved.
34 The implementation of a Helmholtz equation solver will be discussed in \Sec{DIFFUSION HELM SEC}.
35 In Section~\ref{DIFFUSION TRANS SEC} the solver of the Helmholtz equation is used to build a
36 solver for the temperature diffusion problem.
37
38 \subsection{\label{DIFFUSION TEMP SEC}Temperature Diffusion}
39 The unknown temperature $T$ is a function of its location in the domain and time $t>0$. The governing equation
40 in the interior of the domain is given by
41 \begin{equation}
42 \rho c\hackscore p T\hackscore{,t} - (\kappa T\hackscore{,i})\hackscore{,i} = q\hackscore H
43 \label{DIFFUSION TEMP EQ 1}
44 \end{equation}
45 where $\rho c\hackscore p$ and $\kappa$ are given material constants. In case of a composite
46 material the parameters depend on their location in the domain. $q\hackscore H$ is
47 a heat source (or sink) within the domain. We are using the Einstein summation convention \index{summation convention}
48 as introduced in \Chap{FirstSteps}. In our case we assume $q\hackscore H$ to be equal to a constant heat production rate
49 $q^{c}$ on a circle or sphere with center $x^c$ and radius $r$ and $0$ elsewhere:
50 \begin{equation}
51 q\hackscore H(x,t)=
52 \left\{
53 \begin{array}{lcl}
54 q^c & & \|x-x^c\| \le r \\
55 & \mbox{if} \\
56 0 & & \mbox{else} \\
57 \end{array}
58 \right.
59 \label{DIFFUSION TEMP EQ 1b}
60 \end{equation}
61 for all $x$ in the domain and all time $t>0$.
62
63 On the surface of the domain we are
64 specifying a radiation condition
65 which prescribes the normal component of the flux $\kappa T\hackscore{,i}$ to be proportional
66 to the difference of the current temperature to the surrounding temperature $T\hackscore{ref}$:
67 \begin{equation}
68 \kappa T\hackscore{,i} n\hackscore i = \eta (T\hackscore{ref}-T)
69 \label{DIFFUSION TEMP EQ 2}
70 \end{equation}
71 $\eta$ is a given material coefficient depending on the material of the block and the surrounding medium.
72 $n\hackscore i$ is the $i$-th component of the outer normal field \index{outer normal field}
73 at the surface of the domain.
74
75 To solve the time-dependent \eqn{DIFFUSION TEMP EQ 1} the initial temperature at time
76 $t=0$ has to be given. Here we assume that the initial temperature is the surrounding temperature:
77 \begin{equation}
78 T(x,0)=T\hackscore{ref}
79 \label{DIFFUSION TEMP EQ 4}
80 \end{equation}
81 for all $x$ in the domain. It is pointed out that
82 the initial conditions satisfy the
83 boundary condition defined by \eqn{DIFFUSION TEMP EQ 2}.
84
85 The temperature is calculated at discrete time nodes $t^{(n)}$ where
86 $t^{(0)}=0$ and $t^{(n)}=t^{(n-1)}+h$ where $h>0$ is the step size which is assumed to be constant.
87 In the following the upper index ${(n)}$ refers to a value at time $t^{(n)}$. The simplest
88 and most robust scheme to approximate the time derivative of the the temperature is
89 the backward Euler
90 \index{backward Euler} scheme. The backward Euler
91 scheme is based
92 on the Taylor expansion of $T$ at time $t^{(n)}$:
93 \begin{equation}
94 T^{(n)}\approx T^{(n-1)}+T\hackscore{,t}^{(n)}(t^{(n)}-t^{(n-1)})
95 =T^{(n-1)} + h \cdot T\hackscore{,t}^{(n)}
96 \label{DIFFUSION TEMP EQ 6}
97 \end{equation}
98 This is inserted into \eqn{DIFFUSION TEMP EQ 1}. By separating the terms at
99 $t^{(n)}$ and $t^{(n-1)}$ one gets for $n=1,2,3\ldots$
100 \begin{equation}
101 \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^{(n-1)}
102 \label{DIFFUSION TEMP EQ 7}
103 \end{equation}
104 where $T^{(0)}=T\hackscore{ref}$ is taken form the initial condition given by \eqn{DIFFUSION TEMP EQ 4}.
105 Together with the natural boundary condition
106 \begin{equation}
107 \kappa T\hackscore{,i}^{(n)} n\hackscore i = \eta (T\hackscore{ref}-T^{(n)})
108 \label{DIFFUSION TEMP EQ 2222}
109 \end{equation}
110 taken from \eqn{DIFFUSION TEMP EQ 2}
111 this forms a boundary value problem that has to be solved for each time step.
112 As a first step to implement a solver for the temperature diffusion problem we will
113 first implement a solver for the boundary value problem that has to be solved at each time step.
114
115 \subsection{\label{DIFFUSION HELM SEC}Helmholtz Problem}
116 The partial differential equation to be solved for $T^{(n)}$ has the form
117 \begin{equation}
118 \omega T^{(n)} - (\kappa T^{(n)}\hackscore{,i})\hackscore{,i} = f
119 \label{DIFFUSION HELM EQ 1}
120 \end{equation}
121 and we set
122 \begin{equation}
123 \omega=\frac{\rho c\hackscore p}{h} \mbox{ and } f=q\hackscore H +\frac{\rho c\hackscore p}{h}T^{(n-1)} \;.
124 \label{DIFFUSION HELM EQ 1b}
125 \end{equation}
126 With $g=\eta T\hackscore{ref}$ the radiation condition defined by \eqn{DIFFUSION TEMP EQ 2222}
127 takes the form
128 \begin{equation}
129 \kappa T^{(n)}\hackscore{,i} n\hackscore{i} = g - \eta T^{(n)}\mbox{ on } \Gamma
130 \label{DIFFUSION HELM EQ 2}
131 \end{equation}
132 The partial differential \eqn{DIFFUSION HELM EQ 1} together with boundary conditions of \eqn{DIFFUSION HELM EQ 2}
133 is called the Helmholtz equation \index{Helmholtz equation}.
134
135 We want to use the \LinearPDE class provided by \escript to define and solve a general linear,steady, second order PDE such as the
136 Helmholtz equation. For a single PDE the \LinearPDE class supports the following form:
137 \begin{equation}\label{LINEARPDE.SINGLE.1 TUTORIAL}
138 -(A\hackscore{jl} u\hackscore{,l})\hackscore{,j}+D u = Y \; .
139 \end{equation}
140 where we show only the coefficients relevant for the problem discussed here. For the general form of
141 single PDE see \eqn{LINEARPDE.SINGLE.1}.
142 The coefficients $A$, and $Y$ have to be specified through \Data objects in the
143 \Function on the PDE or objects that can be converted into such \Data objects.
144 $A$ is a \RankTwo and $D$ and $Y$ are scalar.
145 The following natural
146 boundary conditions are considered \index{boundary condition!natural} on $\Gamma$:
147 \begin{equation}\label{LINEARPDE.SINGLE.2 TUTORIAL}
148 n\hackscore{j}A\hackscore{jl} u\hackscore{,l}+d u= y \;.
149 \end{equation}
150 Notice that the coefficient $A$ is the same like in the PDE~\eqn{LINEARPDE.SINGLE.1 TUTORIAL}.
151 The coefficients $d$ and $y$ are
152 each a \Scalar in the \FunctionOnBoundary. Constraints \index{constraint} for the solution prescribing the value of the
153 solution at certain locations in the domain. They have the form
154 \begin{equation}\label{LINEARPDE.SINGLE.3 TUTORIAL}
155 u=r \mbox{ where } q>0
156 \end{equation}
157 $r$ and $q$ are each \Scalar where $q$ is the characteristic function
158 \index{characteristic function} defining where the constraint is applied.
159 The constraints defined by \eqn{LINEARPDE.SINGLE.3 TUTORIAL} override any other condition set by
160 \eqn{LINEARPDE.SINGLE.1 TUTORIAL} or \eqn{LINEARPDE.SINGLE.2 TUTORIAL}.
161 The \Poisson class of the \linearPDEs module,
162 which we have already used in \Chap{FirstSteps}, is in fact a subclass of the more general
163 \LinearPDE class. The \linearPDEs module provides a \Helmholtz class but
164 we will make direct use of the general \LinearPDE class.
165
166 By inspecting the Helmholtz equation \index{Helmholtz equation}
167 (\ref{DIFFUSION HELM EQ 1}) and boundary condition (\ref{DIFFUSION HELM EQ 2}) and
168 substituting $u$ for $T^{(n)}$
169 we can easily assign values to the coefficients in the
170 general PDE of the \LinearPDE class:
171 \begin{equation}\label{DIFFUSION HELM EQ 3}
172 \begin{array}{llllll}
173 A\hackscore{ij}=\kappa \delta\hackscore{ij} & D=\omega & Y=f \\
174 d=\eta & y= g & \\
175 \end{array}
176 \end{equation}
177 $\delta\hackscore{ij}$ is the Kronecker symbol \index{Kronecker symbol} defined by $\delta\hackscore{ij}=1$ for
178 $i=j$ and $0$ otherwise. Undefined coefficients are assumed to be not present.\footnote{There is a difference
179 in \escript of being not present and set to zero. As not present coefficients are not processed,
180 it is more efficient to leave a coefficient undefined instead of assigning zero to it.}
181 In this diffusion example we do not need to define a characteristic function $q$ because the
182 boundary conditions we consider in \eqn{DIFFUSION HELM EQ 2} are just the natural boundary
183 conditions which are already defined in the \LinearPDE class (shown in \eqn{LINEARPDE.SINGLE.2 TUTORIAL}).
184
185 The Helmholtz equation can be set up by following way~\footnote{Please, note that this is not a complete code. The complete code can be found in ``helmholtz.py''. } :
186 \begin{python}
187 mypde=LinearPDE(mydomain)
188 mypde.setValue(A=kappa*kronecker(mydomain),D=omega,Y=f,d=eta,y=g)
189 u=mypde.getSolution()
190 \end{python}
191 where we assume that \code{mydomain} is a \Domain object and
192 \code{kappa} \code{omega} \code{eta} and \code{g} are given scalar values
193 typically \code{float} or \Data objects. The \method{setValue} method
194 assigns values to the coefficients of the general PDE. The \method{getSolution} method solves
195 the PDE and returns the solution \code{u} of the PDE. \function{kronecker} is \escript function
196 returning the Kronecker symbol.
197
198 The coefficients can set by several calls of \method{setValue} where the order can be chosen arbitrarily.
199 If a value is assigned to a coefficient several times, the last assigned value is used when
200 the solution is calculated:
201 \begin{python}
202 mypde=LinearPDE(mydomain)
203 mypde.setValue(A=kappa*kronecker(mydomain),d=eta)
204 mypde.setValue(D=omega,Y=f,y=g)
205 mypde.setValue(d=2*eta) # overwrites d=eta
206 u=mypde.getSolution()
207 \end{python}
208 In some cases the solver of the PDE can make use of the fact that the PDE is symmetric\index{symmetric PDE} where the
209 PDE is called symmetric if
210 \begin{equation}\label{LINEARPDE.SINGLE.4 TUTORIAL}
211 A\hackscore{jl}=A\hackscore{lj}\;.
212 \end{equation}
213 Note that $D$ and $d$ may have any value and the right hand sides $Y$, $y$ as well as the constraints
214 are not relevant. The Helmholtz problem is symmetric.
215 The \LinearPDE class provides the method \method{checkSymmetry} to check if the given PDE is symmetric.
216 \begin{python}
217 mypde=LinearPDE(mydomain)
218 mypde.setValue(A=kappa*kronecker(mydomain),d=eta)
219 print mypde.checkSymmetry() # returns True
220 mypde.setValue(B=kronecker(mydomain)[0])
221 print mypde.checkSymmetry() # returns False
222 mypde.setValue(C=kronecker(mydomain)[0])
223 print mypde.checkSymmetry() # returns True
224 \end{python}
225 Unfortunately, a \method{checkSymmetry} is very expensive and is only recommended only for
226 testing and debugging purposes only. The \method{setSymmetryOn} method is used to
227 declare a PDE symmetric:
228 \begin{python}
229 mypde = LinearPDE(mydomain)
230 mypde.setValue(A=kappa*kronecker(mydomain))
231 mypde.setSymmetryOn()
232 \end{python}
233 Now we want to see how we actually solve the Helmholtz equation.
234 on a rectangular domain
235 of length $l\hackscore{0}=5$ and height $l\hackscore{1}=1$. We choose a simple test solution such that we
236 can verify the returned solution against the exact answer. Actually, we
237 take $T=x\hackscore{0}$ (here $q\hackscore H = 0$) and then calculate the right hand side terms $f$ and $g$ such that
238 the test solution becomes the solution of the problem. If we assume $\kappa$ as being constant,
239 an easy calculation shows that we have to choose $f=\omega \cdot x\hackscore{0}$. On the boundary we get
240 $\kappa n\hackscore{i} u\hackscore{,i}=\kappa n\hackscore{0}$.
241 So we have to set $g=\kappa n\hackscore{0}+\eta x\hackscore{0}$. The following script \file{helmholtz.py}
242 \index{scripts!\file{helmholtz.py}} which is available in the \ExampleDirectory
243 implements this test problem using the \finley PDE solver:
244 \begin{python}
245 from esys.escript import *
246 from esys.escript.linearPDEs import LinearPDE
247 from esys.finley import Rectangle
248 #... set some parameters ...
249 kappa=1.
250 omega=0.1
251 eta=10.
252 #... generate domain ...
253 mydomain = Rectangle(l0=5.,l1=1.,n0=50, n1=10)
254 #... open PDE and set coefficients ...
255 mypde=LinearPDE(mydomain)
256 mypde.setSymmetryOn()
257 n=mydomain.getNormal()
258 x=mydomain.getX()
259 mypde.setValue(A=kappa*kronecker(mydomain),D=omega,Y=omega*x[0], \
260 d=eta,y=kappa*n[0]+eta*x[0])
261 #... calculate error of the PDE solution ...
262 u=mypde.getSolution()
263 print "error is ",Lsup(u-x[0])
264 saveVTK("x0.xml",sol=u)
265 \end{python}
266 To visualize the solution `x0.~xml' just use the command
267 \begin{python}
268 mayavi -d u.xml -m SurfaceMap &
269 \end{python}
270 and it is easy to see that the solution $T=x\hackscore{0}$ is calculated.
271
272 The script is similar to the script \file{poisson.py} discussed in \Chap{FirstSteps}.
273 \code{mydomain.getNormal()} returns the outer normal field on the surface of the domain. The function \function{Lsup}
274 imported by the \code{from esys.escript import *} statement and returns the maximum absolute value of its argument.
275 The error shown by the print statement should be in the order of $10^{-7}$. As piecewise bi-linear interpolation is
276 used by \finley approximate the solution and our solution is a linear function of the spatial coordinates one might
277 expect that the error would be zero or in the order of machine precision (typically $\approx 10^{-15}$).
278 However most PDE packages use an iterative solver which is terminated
279 when a given tolerance has been reached. The default tolerance is $10^{-8}$. This value can be altered by using the
280 \method{setTolerance} of the \LinearPDE class.
281
282 \subsection{The Transition Problem}
283 \label{DIFFUSION TRANS SEC}
284 Now we are ready to solve the original time-dependent problem. The main
285 part of the script is the loop over time $t$ which takes the following form:
286 \begin{python}
287 t=0
288 T=Tref
289 mypde=LinearPDE(mydomain)
290 mypde.setValue(A=kappa*kronecker(mydomain),D=rhocp/h,d=eta,y=eta*Tref)
291 while t<t_end:
292 mypde.setValue(Y=q+rhocp/h*T)
293 T=mypde.getSolution()
294 t+=h
295 \end{python}
296 \var{kappa}, \var{rhocp}, \var{eta} and \var{Tref} are input parameters of the model. \var{q} is the heat source
297 in the domain and \var{h} is the time step size.
298 The variable \var{T}
299 holds the current temperature. It is used to calculate the right hand side coefficient \var{f} in the
300 Helmholtz equation in \eqn{DIFFUSION HELM EQ 1}. The statement \code{T=mypde.getSolution()} overwrites \var{T} with the
301 temperature of the new time step $\var{t}+\var{h}$. To get this iterative process going we need to specify the
302 initial temperature distribution, which equal to $T\hackscore{ref}$.
303 The \LinearPDE class object \var{mypde}
304 and coefficients not changing over time are set up before the loop over time is entered. In each time step only the coefficient
305 $Y$ is reset as it depends on the temperature of the previous time step. This allows the PDE solver to reuse information
306 from previous time steps as much as possible.
307
308 The heat source $q\hackscore H$ which is defined in \eqn{DIFFUSION TEMP EQ 1b} is \var{qc}
309 in an area defined as a circle of radius \var{r} and center \var{xc} and zero outside this circle.
310 \var{q0} is a fixed constant. The following script defines $q\hackscore H$ as desired:
311 \begin{python}
312 from esys.escript import length,whereNegative
313 xc=[0.02,0.002]
314 r=0.001
315 x=mydomain.getX()
316 qH=q0*whereNegative(length(x-xc)-r)
317 \end{python}
318 \var{x} is a \Data class object of
319 the \escript module defining locations in the \Domain \var{mydomain}.
320 The \function{length()} imported from the \escript module returns the
321 Euclidean norm:
322 \begin{equation}\label{DIFFUSION HELM EQ 3aba}
323 \|y\|=\sqrt{
324 y\hackscore{i}
325 y\hackscore{i}
326 } = \function{esys.escript.length}(y)
327 \end{equation}
328 So \code{length(x-xc)} calculates the distances
329 of the location \var{x} to the center of the circle \var{xc} where the heat source is acting.
330 Note that the coordinates of \var{xc} are defined as a list of floating point numbers. It is automatically
331 converted into a \Data class object before being subtracted from \var{x}. The function \function{whereNegative}
332 applied to
333 \code{length(x-xc)-r}, returns a \Data object which is equal to one where the object is negative (inside the circle) and
334 zero elsewhere. After multiplication with \var{qc} we get a function with the desired property of having value \var{qc} inside
335 the circle and zero elsewhere.
336
337 Now we can put the components together to create the script \file{diffusion.py} which is available in the \ExampleDirectory:
338 \index{scripts!\file{diffusion.py}}:
339 \begin{python}
340 from esys.escript import *
341 from esys.escript.linearPDEs import LinearPDE
342 from esys.finley import Rectangle
343 #... set some parameters ...
344 xc=[0.02,0.002]
345 r=0.001
346 qc=50.e6
347 Tref=0.
348 rhocp=2.6e6
349 eta=75.
350 kappa=240.
351 tend=5.
352 # ... time, time step size and counter ...
353 t=0
354 h=0.1
355 i=0
356 #... generate domain ...
357 mydomain = Rectangle(l0=0.05,l1=0.01,n0=250, n1=50)
358 #... open PDE ...
359 mypde=LinearPDE(mydomain)
360 mypde.setSymmetryOn()
361 mypde.setValue(A=kappa*kronecker(mydomain),D=rhocp/h,d=eta,y=eta*Tref)
362 # ... set heat source: ....
363 x=mydomain.getX()
364 qH=qc*whereNegative(length(x-xc)-r)
365 # ... set initial temperature ....
366 T=Tref
367 # ... start iteration:
368 while t<tend:
369 i+=1
370 t+=h
371 print "time step :",t
372 mypde.setValue(Y=qH+rhocp/h*T)
373 T=mypde.getSolution()
374 saveVTK("T.%d.xml"%i,temp=T)
375 \end{python}
376 The script will create the files \file{T.1.xml},
377 \file{T.2.xml}, $\ldots$, \file{T.50.xml} in the directory where the script has been started. The files give the
378 temperature distributions at time steps $1$, $2$, $\ldots$, $50$ in the \VTK file format.
379
380 \begin{figure}
381 \centerline{\includegraphics[width=\figwidth]{figures/DiffusionRes1}}
382 \centerline{\includegraphics[width=\figwidth]{figures/DiffusionRes16}}
383 \centerline{\includegraphics[width=\figwidth]{figures/DiffusionRes32}}
384 \centerline{\includegraphics[width=\figwidth]{figures/DiffusionRes48}}
385 \caption{Results of the Temperature Diffusion Problem for Time Steps $1$ $16$, $32$ and $48$.}
386 \label{DIFFUSION FIG 2}
387 \end{figure}
388 \fig{DIFFUSION FIG 2} shows the result for some selected time steps.
389 An easy way to visualize the results is the command
390 \begin{verbatim}
391 mayavi -d T.1.xml -m SurfaceMap &
392 \end{verbatim}
393 Use the \texttt{Configure Data} window in mayavi
394 to move forward and and backwards in time.

Properties

Name Value
svn:eol-style native
svn:keywords Author Date Id Revision

  ViewVC Help
Powered by ViewVC 1.1.26