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

Revision 575 - (show annotations)
Fri Mar 3 03:33:07 2006 UTC (14 years, 4 months ago) by lkettle
File MIME type: application/x-tex
File size: 18371 byte(s)
I have changed some of the documentation and added more explanations for
the online reference guide for esys13. I have modified two of the
example source codes to write out the results for Helmholtz problem and
changed one variable name in the diffusion.py code to avoid confusion.


 1 % $Id$ 2 \section{The Diffusion Problem} 3 \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 \subsection{\label{DIFFUSION OUT SEC}Outline} 12 In this chapter we will discuss how to solve the time dependent-temperature diffusion\index{diffusion equation} for 13 a block of material. Within the block there is a heat source which drives the temperature diffusion. 14 On the surface, energy can radiate into the surrounding environment. 15 \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 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 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 \subsection{\label{DIFFUSION TEMP SEC}Temperature Diffusion} 26 The unknown temperature $T$ is a function of its location in the domain and time $t>0$. The governing equation 27 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\hackscore H 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 material the parameters depend on their location in the domain. $q\hackscore H$ is 34 a heat source (or sink) within the domain. We are using Einstein summation convention \index{summation convention} 35 as introduced in \Chap{FirstSteps}. In our case we assume $q\hackscore H$ 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 \begin{equation} 38 q\hackscore H(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 specifying a radiation condition 52 which precribes the normal component of the flux $\kappa T\hackscore{,i}$ to be proportional 53 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 $\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 at the surface of the domain. 61 62 To solve the time dependent \eqn{DIFFUSION TEMP EQ 1} the initial temperature at time 63 $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 the initial conditions satisfy the 70 boundary condition defined by \eqn{DIFFUSION TEMP EQ 2}. 71 72 The temperature is calculated at discrete time nodes $t^{(n)}$ where 73 $t^{(0)}=0$ and $t^{(n)}=t^{(n-1)}+h$ where $h>0$ is the step size which is assumed to be constant. 74 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 \begin{equation} 81 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 \label{DIFFUSION TEMP EQ 6} 84 \end{equation} 85 This is inserted into \eqn{DIFFUSION TEMP EQ 1}. By separating the terms at 86 $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\hackscore H + \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 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 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 \subsection{\label{DIFFUSION HELM SEC}Helmholtz Problem} 103 The partial differential equation to be solved for $T^{(n)}$ has the form 104 \begin{equation} 105 \omega T^{(n)} - (\kappa T^{(n)}\hackscore{,i})\hackscore{,i} = f 106 \label{DIFFUSION HELM EQ 1} 107 \end{equation} 108 and we set 109 \begin{equation} 110 \omega=\frac{\rho c\hackscore p}{h} \mbox{ and } f=q\hackscore H +\frac{\rho c\hackscore p}{h}T^{(n-1)} \;. 111 \label{DIFFUSION HELM EQ 1b} 112 \end{equation} 113 With $g=\eta T\hackscore{ref}$ the radiation condition defined by \eqn{DIFFUSION TEMP EQ 2222} 114 takes the form 115 \begin{equation} 116 \kappa T^{(n)}\hackscore{,i} n\hackscore{i} = g - \eta T^{(n)}\mbox{ on } \Gamma 117 \label{DIFFUSION HELM EQ 2} 118 \end{equation} 119 The partial differential \eqn{DIFFUSION HELM EQ 1} together with boundary conditions of \eqn{DIFFUSION HELM EQ 2} 120 is called the Helmholtz equation \index{Helmholtz equation}. 121 122 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 \end{equation} 127 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 \end{equation} 135 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 150 By inspecting the Helmholtz equation \index{Helmholtz equation} 151 (\ref{DIFFUSION HELM EQ 1}) and boundary condition (\ref{DIFFUSION HELM EQ 2}) and 152 substituting $u$ for $T^{(n)}$ 153 we can easily assign values to the coefficients in the 154 general PDE of the \LinearPDE class: 155 \begin{equation}\label{DIFFUSION HELM EQ 3} 156 \begin{array}{llllll} 157 A\hackscore{ij}=\kappa \delta\hackscore{ij} & D=\omega & Y=f \\ 158 d=\eta & y= g & \\ 159 \end{array} 160 \end{equation} 161 $\delta\hackscore{ij}$ is the Kronecker symbol \index{Kronecker symbol} defined by $\delta\hackscore{ij}=1$ for 162 $i=j$ and $0$ otherwise. Undefined coefficients are assumed to be not present.\footnote{There is a difference 163 in \escript of being not present and set to zero. As not present coefficients are not processed, 164 it is more efficient to leave a coefficient undefined instead of assigning zero to it.} 165 In this diffusion example we do not need to define a characteristic function $q$ because the 166 boundary conditions we consider in \eqn{DIFFUSION HELM EQ 2} are just the natural boundary 167 conditions which are already defined in the \LinearPDE class (shown in \eqn{LINEARPDE.SINGLE.2}). 168 169 Defining and solving the Helmholtz equation is very easy now: 170 \begin{python} 171 from esys.escript import * 172 from linearPDEs import LinearPDE 173 mypde=LinearPDE(mydomain) 174 mypde.setValue(A=kappa*kronecker(mydomain),D=omega,Y=f,d=eta,y=g) 175 u=mypde.getSolution() 176 \end{python} 177 where we assume that \code{mydomain} is a \Domain object and 178 \code{kappa} \code{omega} \code{eta} and \code{g} are given scalar values 179 typically \code{float} or \Data objects. The \method{setValue} method 180 assigns values to the coefficients of the general PDE. The \method{getSolution} method solves 181 the PDE and returns the solution \code{u} of the PDE. \function{kronecker} is \escript function 182 returning the Kronecker symbol. 183 184 The coefficients can set by several calls of \method{setValue} where the order can be chosen arbitrarily. 185 If a value is assigned to a coefficient several times, the last assigned value is used when 186 the solution is calculated: 187 \begin{python} 188 mypde=LinearPDE(mydomain) 189 mypde.setValue(A=kappa*kronecker(mydomain),d=eta) 190 mypde.setValue(D=omega,Y=f,y=g) 191 mypde.setValue(d=2*eta) # overwrites d=eta 192 u=mypde.getSolution() 193 \end{python} 194 In some cases the solver of the PDE can make use of the fact that the PDE is symmetric\index{symmetric PDE} where the 195 PDE is called symmetric if 196 \begin{equation}\label{LINEARPDE.SINGLE.4} 197 A\hackscore{jl}=A\hackscore{lj} \mbox{ and } B\hackscore{j}=C\hackscore{j} \;. 198 \end{equation} 199 Note that $D$ and $d$ may have any value and the right hand side $X$, $Y$, $y$ as well as the constraints 200 are not relevant. The Helmholtz problem is symmetric. 201 The \LinearPDE class provides the method \method{checkSymmetry} method to check if the given PDE is symmetric. 202 \begin{python} 203 mypde=LinearPDE(mydomain) 204 mypde.setValue(A=kappa*kronecker(mydomain),d=eta) 205 print mypde.checkSymmetry() # returns True 206 mypde.setValue(B=kronecker(mydomain)) 207 print mypde.checkSymmetry() # returns False 208 mypde.setValue(C=kronecker(mydomain)) 209 print mypde.checkSymmetry() # returns True 210 \end{python} 211 Unfortunately, a \method{checkSymmetry} is very expensive and is recommendable to use for 212 testing and debugging purposes only. The \method{setSymmetryOn} method is used to 213 declare a PDE symmetric: 214 \begin{python} 215 mypde = LinearPDE(mydomain) 216 mypde.setValue(A=kappa*kronecker(mydomain)) 217 mypde.setSymmetryOn() 218 \end{python} 219 Now we want to see how we actually solve the Helmholtz equation. 220 on a rectangular domain 221 of length $l\hackscore{0}=5$ and height $l\hackscore{1}=1$. We choose a simple test solution such that we 222 can verify the returned solution against the exact answer. Actually, we 223 take $T=x\hackscore{0}$ (here $q\hackscore H = 0$) and then calculate the right hand side terms $f$ and $g$ such that 224 the test solution becomes the solution of the problem. If we assume $\kappa$ as being constant, 225 an easy calculation shows that we have to choose $f=\omega \cdot x\hackscore{0}$. On the boundary we get 226 $\kappa n\hackscore{i} u\hackscore{,i}=\kappa n\hackscore{0}$. 227 So we have to set $g=\kappa n\hackscore{0}+\eta x\hackscore{0}$. The following script \file{helmholtz.py} 228 \index{scripts!\file{helmholtz.py}} which is available in the \ExampleDirectory 229 implements this test problem using the \finley PDE solver: 230 \begin{python} 231 from esys.escript import * 232 from esys.escript.linearPDEs import LinearPDE 233 from esys.finley import Rectangle 234 #... set some parameters ... 235 kappa=1. 236 omega=0.1 237 eta=10. 238 #... generate domain ... 239 mydomain = Rectangle(l0=5.,l1=1.,n0=50, n1=10) 240 #... open PDE and set coefficients ... 241 mypde=LinearPDE(mydomain) 242 mypde.setSymmetryOn() 243 n=mydomain.getNormal() 244 x=mydomain.getX() 245 mypde.setValue(A=kappa*kronecker(mydomain),D=omega,Y=omega*x, \ 246 d=eta,y=kappa*n+eta*x) 247 #... calculate error of the PDE solution ... 248 u=mypde.getSolution() 249 print "error is ",Lsup(u-x) 250 saveVTK("x0.xml",sol=u) 251 \end{python} 252 To visualise the solution `x0.~xml' just use the command 253 \begin{python} 254 mayavi -d u.xml -m SurfaceMap & 255 \end{python} 256 and it is easy to see that the solution $T=x\hackscore{0}$ is calculated. 257 258 The script is similar to the script \file{poisson.py} dicussed in \Chap{FirstSteps}. 259 \code{mydomain.getNormal()} returns the outer normal field on the surface of the domain. The function \function{Lsup} 260 imported by the \code{from esys.escript import *} statement and returns the maximum absolute value of its argument. 261 The error shown by the print statement should be in the order of $10^{-7}$. As piecewise bi-linear interpolation is 262 used by \finley approximate the solution and our solution is a linear function of the spatial coordinates one might 263 expect that the error would be zero or in the order of machine precision (typically $\approx 10^{-15}$). 264 However most PDE packages use an iterative solver which is terminated 265 when a given tolerance has been reached. The default tolerance is $10^{-8}$. This value can be altered by using the 266 \method{setTolerance} of the \LinearPDE class. 267 268 \subsection{The Transition Problem} 269 \label{DIFFUSION TRANS SEC} 270 Now we are ready to solve the original time dependent problem. The main 271 part of the script is the loop over time $t$ which takes the following form: 272 \begin{python} 273 t=0 274 T=Tref 275 mypde=LinearPDE(mydomain) 276 mypde.setValue(A=kappa*kronecker(mydomain),D=rhocp/h,d=eta,y=eta*Tref) 277 while t

## Properties

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