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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 107 - (hide annotations)
Thu Jan 27 06:21:48 2005 UTC (14 years, 5 months ago) by jgs
Original Path: trunk/esys2/doc/user/diffusion.tex
File MIME type: application/x-tex
File size: 18183 byte(s)
commit of branch dev-01 back to main trunk on 2005-01-27

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 dependent-temperature 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 & & \|x-x^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^{(n-1)}+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^{(n-1)}\approx T^{(n)}+T\hackscore{,t}^{(n)}(t^{(n-1)}-t^{(n)})
83 jgs 107 =T^{(n-1)} - 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^{(n-1)}$ 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^{(n-1)}
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^{(n-1)} \;.
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 sub-class 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(u-x[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 bi-linear 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(x-xc)-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(x-xc)} 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(x-xc)-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(x-xc)-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.

Properties

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

  ViewVC Help
Powered by ViewVC 1.1.26