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

Revision 2335 - (hide annotations)
Thu Mar 26 04:33:44 2009 UTC (10 years, 6 months ago) by jfenwick
File MIME type: application/x-tex
File size: 18913 byte(s)
The user guide now builds under pdflatex (if you have converted the figures). Unfortunately, it has a really ugly title.

 1 ksteube 1811 2 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 3 ksteube 1316 % 4 ksteube 1811 % Copyright (c) 2003-2008 by University of Queensland 5 % Earth Systems Science Computational Center (ESSCC) 6 7 gross 625 % 8 ksteube 1811 % Primary Business: Queensland, Australia 9 % Licensed under the Open Software License version 3.0 10 11 gross 625 % 12 ksteube 1811 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 13 gross 625 14 ksteube 1811 15 jgs 121 \section{The Diffusion Problem} 16 jgs 102 \label{DIFFUSION CHAP} 17 18 \begin{figure} 19 jfenwick 2335 \centerline{\includegraphics[width=\figwidth]{figures/DiffusionDomain}} 20 jgs 102 \caption{Temperature Diffusion Problem with Circular Heat Source} 21 \label{DIFFUSION FIG 1} 22 \end{figure} 23 24 jgs 121 \subsection{\label{DIFFUSION OUT SEC}Outline} 25 ksteube 1316 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 jgs 107 On the surface, energy can radiate into the surrounding environment. 28 jgs 102 \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 jgs 107 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 jgs 102 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 jgs 121 \subsection{\label{DIFFUSION TEMP SEC}Temperature Diffusion} 39 jgs 107 The unknown temperature $T$ is a function of its location in the domain and time $t>0$. The governing equation 40 jgs 102 in the interior of the domain is given by 41 \begin{equation} 42 lkettle 575 \rho c\hackscore p T\hackscore{,t} - (\kappa T\hackscore{,i})\hackscore{,i} = q\hackscore H 43 jgs 102 \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 lkettle 575 material the parameters depend on their location in the domain. $q\hackscore H$ is 47 ksteube 1316 a heat source (or sink) within the domain. We are using the Einstein summation convention \index{summation convention} 48 lkettle 575 as introduced in \Chap{FirstSteps}. In our case we assume $q\hackscore H$ to be equal to a constant heat production rate 49 jgs 107 $q^{c}$ on a circle or sphere with center $x^c$ and radius $r$ and $0$ elsewhere: 50 jgs 102 \begin{equation} 51 lkettle 575 q\hackscore H(x,t)= 52 jgs 102 \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 jgs 107 specifying a radiation condition 65 ksteube 1316 which prescribes the normal component of the flux $\kappa T\hackscore{,i}$ to be proportional 66 jgs 102 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 jgs 107 $\eta$ is a given material coefficient depending on the material of the block and the surrounding medium. 72 ksteube 1316 $n\hackscore i$ is the $i$-th component of the outer normal field \index{outer normal field} 73 jgs 102 at the surface of the domain. 74 75 ksteube 1316 To solve the time-dependent \eqn{DIFFUSION TEMP EQ 1} the initial temperature at time 76 jgs 102 $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 jgs 107 the initial conditions satisfy the 83 jgs 102 boundary condition defined by \eqn{DIFFUSION TEMP EQ 2}. 84 85 jgs 107 The temperature is calculated at discrete time nodes $t^{(n)}$ where 86 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. 87 jgs 107 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 gross 660 \index{backward Euler} scheme. The backward Euler 91 jgs 107 scheme is based 92 on the Taylor expansion of $T$ at time $t^{(n)}$: 93 jgs 102 \begin{equation} 94 lkettle 573 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 jgs 102 \label{DIFFUSION TEMP EQ 6} 97 \end{equation} 98 jgs 107 This is inserted into \eqn{DIFFUSION TEMP EQ 1}. By separating the terms at 99 jgs 102 $t^{(n)}$ and $t^{(n-1)}$ one gets for $n=1,2,3\ldots$ 100 \begin{equation} 101 lkettle 575 \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 jgs 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 jgs 107 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 jgs 102 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 jgs 121 \subsection{\label{DIFFUSION HELM SEC}Helmholtz Problem} 116 jgs 102 The partial differential equation to be solved for $T^{(n)}$ has the form 117 \begin{equation} 118 lkettle 575 \omega T^{(n)} - (\kappa T^{(n)}\hackscore{,i})\hackscore{,i} = f 119 jgs 102 \label{DIFFUSION HELM EQ 1} 120 \end{equation} 121 lkettle 575 and we set 122 jgs 102 \begin{equation} 123 lkettle 575 \omega=\frac{\rho c\hackscore p}{h} \mbox{ and } f=q\hackscore H +\frac{\rho c\hackscore p}{h}T^{(n-1)} \;. 124 jgs 102 \label{DIFFUSION HELM EQ 1b} 125 \end{equation} 126 jgs 107 With $g=\eta T\hackscore{ref}$ the radiation condition defined by \eqn{DIFFUSION TEMP EQ 2222} 127 jgs 102 takes the form 128 \begin{equation} 129 lkettle 575 \kappa T^{(n)}\hackscore{,i} n\hackscore{i} = g - \eta T^{(n)}\mbox{ on } \Gamma 130 jgs 102 \label{DIFFUSION HELM EQ 2} 131 \end{equation} 132 gross 568 The partial differential \eqn{DIFFUSION HELM EQ 1} together with boundary conditions of \eqn{DIFFUSION HELM EQ 2} 133 jgs 107 is called the Helmholtz equation \index{Helmholtz equation}. 134 jgs 102 135 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 136 Helmholtz equation. For a single PDE the \LinearPDE class supports the following form: 137 gross 625 \begin{equation}\label{LINEARPDE.SINGLE.1 TUTORIAL} 138 -(A\hackscore{jl} u\hackscore{,l})\hackscore{,j}+D u = Y \; . 139 jgs 102 \end{equation} 140 gross 625 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 gross 568 \Function on the PDE or objects that can be converted into such \Data objects. 144 gross 625 $A$ is a \RankTwo and $D$ and $Y$ are scalar. 145 gross 568 The following natural 146 boundary conditions are considered \index{boundary condition!natural} on $\Gamma$: 147 gross 625 \begin{equation}\label{LINEARPDE.SINGLE.2 TUTORIAL} 148 n\hackscore{j}A\hackscore{jl} u\hackscore{,l}+d u= y \;. 149 jgs 102 \end{equation} 150 gross 625 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 gross 568 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 gross 625 \begin{equation}\label{LINEARPDE.SINGLE.3 TUTORIAL} 155 gross 568 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 gross 625 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 gross 568 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 jgs 102 166 jgs 107 By inspecting the Helmholtz equation \index{Helmholtz equation} 167 lkettle 575 (\ref{DIFFUSION HELM EQ 1}) and boundary condition (\ref{DIFFUSION HELM EQ 2}) and 168 substituting $u$ for $T^{(n)}$ 169 jgs 102 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 jgs 107 $\delta\hackscore{ij}$ is the Kronecker symbol \index{Kronecker symbol} defined by $\delta\hackscore{ij}=1$ for 178 lkettle 575 $i=j$ and $0$ otherwise. Undefined coefficients are assumed to be not present.\footnote{There is a difference 179 gross 568 in \escript of being not present and set to zero. As not present coefficients are not processed, 180 lkettle 575 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 gross 625 conditions which are already defined in the \LinearPDE class (shown in \eqn{LINEARPDE.SINGLE.2 TUTORIAL}). 184 jgs 102 185 artak 1971 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 gross 568 \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 gross 569 the PDE and returns the solution \code{u} of the PDE. \function{kronecker} is \escript function 196 gross 568 returning the Kronecker symbol. 197 198 The coefficients can set by several calls of \method{setValue} where the order can be chosen arbitrarily. 199 lkettle 573 If a value is assigned to a coefficient several times, the last assigned value is used when 200 gross 568 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 gross 625 \begin{equation}\label{LINEARPDE.SINGLE.4 TUTORIAL} 211 A\hackscore{jl}=A\hackscore{lj}\;. 212 gross 568 \end{equation} 213 gross 625 Note that $D$ and $d$ may have any value and the right hand sides $Y$, $y$ as well as the constraints 214 gross 568 are not relevant. The Helmholtz problem is symmetric. 215 gross 569 The \LinearPDE class provides the method \method{checkSymmetry} method to check if the given PDE is symmetric. 216 gross 568 \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 gross 569 Unfortunately, a \method{checkSymmetry} is very expensive and is recommendable to use for 226 testing and debugging purposes only. The \method{setSymmetryOn} method is used to 227 gross 568 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 gross 569 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 lkettle 575 take $T=x\hackscore{0}$ (here $q\hackscore H = 0$) and then calculate the right hand side terms $f$ and $g$ such that 238 jgs 107 the test solution becomes the solution of the problem. If we assume $\kappa$ as being constant, 239 jgs 102 an easy calculation shows that we have to choose $f=\omega \cdot x\hackscore{0}$. On the boundary we get 240 jgs 107 $\kappa n\hackscore{i} u\hackscore{,i}=\kappa n\hackscore{0}$. 241 gross 569 So we have to set $g=\kappa n\hackscore{0}+\eta x\hackscore{0}$. The following script \file{helmholtz.py} 242 gross 568 \index{scripts!\file{helmholtz.py}} which is available in the \ExampleDirectory 243 jgs 102 implements this test problem using the \finley PDE solver: 244 \begin{python} 245 gross 568 from esys.escript import * 246 gross 569 from esys.escript.linearPDEs import LinearPDE 247 jgs 107 from esys.finley import Rectangle 248 jgs 102 #... set some parameters ... 249 jgs 107 kappa=1. 250 jgs 102 omega=0.1 251 eta=10. 252 #... generate domain ... 253 jgs 107 mydomain = Rectangle(l0=5.,l1=1.,n0=50, n1=10) 254 jgs 102 #... open PDE and set coefficients ... 255 gross 568 mypde=LinearPDE(mydomain) 256 mypde.setSymmetryOn() 257 jgs 102 n=mydomain.getNormal() 258 x=mydomain.getX() 259 gross 569 mypde.setValue(A=kappa*kronecker(mydomain),D=omega,Y=omega*x[0], \ 260 d=eta,y=kappa*n[0]+eta*x[0]) 261 jgs 102 #... calculate error of the PDE solution ... 262 u=mypde.getSolution() 263 print "error is ",Lsup(u-x[0]) 264 lkettle 575 saveVTK("x0.xml",sol=u) 265 jgs 102 \end{python} 266 lgraham 1700 To visualize the solution `x0.~xml' just use the command 267 lkettle 575 \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 lgraham 1700 The script is similar to the script \file{poisson.py} discussed in \Chap{FirstSteps}. 273 jgs 102 \code{mydomain.getNormal()} returns the outer normal field on the surface of the domain. The function \function{Lsup} 274 lkettle 573 imported by the \code{from esys.escript import *} statement and returns the maximum absolute value of its argument. 275 jgs 107 The error shown by the print statement should be in the order of $10^{-7}$. As piecewise bi-linear interpolation is 276 gross 569 used by \finley approximate the solution and our solution is a linear function of the spatial coordinates one might 277 jgs 107 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 jgs 102 \method{setTolerance} of the \LinearPDE class. 281 282 jgs 121 \subsection{The Transition Problem} 283 jgs 102 \label{DIFFUSION TRANS SEC} 284 ksteube 1316 Now we are ready to solve the original time-dependent problem. The main 285 jgs 107 part of the script is the loop over time $t$ which takes the following form: 286 jgs 102 \begin{python} 287 jgs 107 t=0 288 T=Tref 289 gross 569 mypde=LinearPDE(mydomain) 290 mypde.setValue(A=kappa*kronecker(mydomain),D=rhocp/h,d=eta,y=eta*Tref) 291 jgs 102 while t

## Properties

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