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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 573 - (hide annotations)
Thu Mar 2 00:42:53 2006 UTC (13 years, 6 months ago) by lkettle
File MIME type: application/x-tex
File size: 17649 byte(s)
I have made a few changes to the documentation for the online reference
guide for the Poisson and diffusion examples.

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