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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 6678 - (show annotations)
Mon Jun 11 04:05:30 2018 UTC (16 months, 1 week ago) by jfenwick
File MIME type: application/x-tex
File size: 18546 byte(s)
Fixing overfull hboxes

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

Properties

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

  ViewVC Help
Powered by ViewVC 1.1.26