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

Revision 102 - (show annotations)
Wed Dec 15 07:08:39 2004 UTC (17 years, 6 months ago) by jgs
File MIME type: application/x-tex
File size: 17701 byte(s)
*** empty log message ***


 1 % $Id$ 2 \chapter{How to Solve The Diffusion Equation} 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 \begin{figure} 12 \centerline{\includegraphics[width=\figwidth]{DiffusionRes1}} 13 \centerline{\includegraphics[width=\figwidth]{DiffusionRes16}} 14 \centerline{\includegraphics[width=\figwidth]{DiffusionRes32}} 15 \centerline{\includegraphics[width=\figwidth]{DiffusionRes48}} 16 \caption{Results of the Temperture Diffusion Problem for Time Steps $1$ $16$, $32$ and $48$.} 17 \label{DIFFUSION FIG 2} 18 \end{figure} 19 20 21 \section{\label{DIFFUSION OUT SEC}Outline} 22 In this chapter we will discuss how to solve the time depeneded-temperature diffusion\index{diffusion equation} within 23 a block of material. Within the block there is a heat source which drives the temperature diffusion. 24 On the surface energy can radiate into the surrounding environment. 25 \fig{DIFFUSION FIG 1} shows the configuration. 26 27 In the next \Sec{DIFFUSION TEMP SEC} we will present the relevant model. A 28 time integration scheme is introduced to calculate the temperature at given time nodes $t^{(n)}$. 29 We will see that at time step a so-called Helmholtz equation \index{Helmholtz equation} has to be solved. 30 The implementation iof a Helmholtz equation solver will be discussed in \Sec{DIFFUSION HELM SEC}. 31 In Section~\ref{DIFFUSION TRANS SEC} the solver of the Helmholtz equation is used to build a 32 solver for the temperature diffusion problem. 33 34 \section{\label{DIFFUSION TEMP SEC}Temperature Diffusion} 35 36 The temperature $T$ is a function of its location in the domain and time $t>0$. The governing equation 37 in the interior of the domain is given by 38 \begin{equation} 39 \rho c\hackscore p T\hackscore{,t} - (\kappa T\hackscore{,i})\hackscore{,i} = q 40 \label{DIFFUSION TEMP EQ 1} 41 \end{equation} 42 where $\rho c\hackscore p$ and $\kappa$ are given material constants. In case of a composite 43 material the parameters are depending on the their location in the domain. $q$ is 44 a heat source (or sink) within the domain. We are using Einstein summation convention \index{summation convention} 45 as introduced in \Chap{FirstSteps}. In our case we assume $q$ to be equal to a constant $q^{c}$ 46 on a circle or sphere with center $x^c$ and radius $r$ and $0$ elsewhere: 47 \begin{equation} 48 q(x,t)= 49 \left\{ 50 \begin{array}{lcl} 51 q^c & & \|x-x^c\| \le r \\ 52 & \mbox{if} \\ 53 0 & & \mbox{else} \\ 54 \end{array} 55 \right. 56 \label{DIFFUSION TEMP EQ 1b} 57 \end{equation} 58 for all $x$ in the domain and all time $t>0$. 59 60 On the surface of the domain we are 61 are specifying a radiation condition 62 which precribes the normal component of the flux $J\hackscore i= \kappa T\hackscore{,i}$ to be proportional 63 to the difference of the current temperature to the surrounding temperature $T\hackscore{ref}$: 64 \begin{equation} 65 \kappa T\hackscore{,i} n\hackscore i = \eta (T\hackscore{ref}-T) 66 \label{DIFFUSION TEMP EQ 2} 67 \end{equation} 68 $\eta$ is a given material coefficient depending on the material and the surrounding medium. 69 As usual $n_i$ is the $i$-th component of the outer normal field \index{outer normal field} 70 at the surface of the domain. 71 72 To solve the the time depended \eqn{DIFFUSION TEMP EQ 1} the initial temperature at time 73 $t=0$ has to be given. Here we assume that the initial temperature is the surrounding temperature: 74 \begin{equation} 75 T(x,0)=T\hackscore{ref} 76 \label{DIFFUSION TEMP EQ 4} 77 \end{equation} 78 for all $x$ in the domain. It is pointed out that 79 the initial conditions is fullfilling the 80 boundary condition defined by \eqn{DIFFUSION TEMP EQ 2}. 81 82 The temperature is calculated discrete time nodes $t^{(n)}$ where 83 $t^{(0)}=0$ and $t^{(n)}=t^{(n-1)}+h$ where $h>0$ is the step size which is assumed to be constant. 84 In the following the upper index ${(n)}$ is refering to a value at time $t^{(n)}$. The simplest 85 and most robust scheme to approximate the time derivative of the the temperature is the 86 \index{backward Euler} scheme, see~\cite{XXX} for alternatives. The backward Euler scheme bases 87 on the Taylor expansion 88 \begin{equation} 89 T^{(n-1)}\approx T^{(n)}+T\hackscore{,t}^{(n)}(t^{(n-1)}-t^{(n)}) 90 =T^{(n-1)} + h \cdot T\hackscore{,t}^{(n)} 91 \label{DIFFUSION TEMP EQ 6} 92 \end{equation} 93 which is inserted into \eqn{DIFFUSION TEMP EQ 1}. By separating the terms at 94 $t^{(n)}$ and $t^{(n-1)}$ one gets for $n=1,2,3\ldots$ 95 \begin{equation} 96 \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)} 97 \label{DIFFUSION TEMP EQ 7} 98 \end{equation} 99 where $T^{(0)}=T\hackscore{ref}$ is taken form the initial condition given by \eqn{DIFFUSION TEMP EQ 4}. 100 Together with the natural boundary condition from \eqn{DIFFUSION TEMP EQ 2}. 101 this forms a boundary value problem that has to be solved for each time step. 102 As a first step to implement a solver for the temperature diffusion problem we will 103 first implement a solver for the boundary value problem that has to be solved at each time step. 104 105 \section{\label{DIFFUSION HELM SEC}Helmholtz Problem} 106 The partial differential equation to be solved for $T^{(n)}$ has the form 107 \begin{equation} 108 \omega u - (\kappa u\hackscore{,i})\hackscore{,i} = f 109 \label{DIFFUSION HELM EQ 1} 110 \end{equation} 111 where $u$ plays the role of $T^{(n)}$ and we set 112 \begin{equation} 113 \omega=\frac{\rho c\hackscore p}{h} \mbox{ and } f=q+\frac{\rho c\hackscore p}{h}T^{(n-1)} \;. 114 \label{DIFFUSION HELM EQ 1b} 115 \end{equation} 116 With $g=\eta T\hackscore{ref}$ the radiation condition defined by \eqn{DIFFUSION TEMP EQ 2} 117 takes the form 118 \begin{equation} 119 \kappa u\hackscore{,i} n\hackscore{i} = g - \eta u\mbox{ on } \Gamma 120 \label{DIFFUSION HELM EQ 2} 121 \end{equation} 122 The partial differential 123 \eqn{DIFFUSION HELM EQ 1} together with boundary conditions of \eqn{DIFFUSION HELM EQ 2} 124 is called a Helmholtz equation \index{Helmholtz equation}. 125 126 We want to use the \LinearPDE class provided by \escript to define and solve a general linear PDE such as the 127 Helmholtz equation. We have used a special case of the \LinearPDE class, namely the 128 \Poisson class already in \Chap{FirstSteps}. 129 Here we will write our own specialized class of the \LinearPDE to solve the Helmholtz equation. 130 131 The general form of a single PDE that can be handeled by the \LinearPDE class is 132 \begin{equation}\label{EQU.FEM.1} 133 -(A\hackscore{jl} u\hackscore{,l})\hackscore{,j}+D u =Y \; . 134 \end{equation} 135 The general form and systems is discussed in \Sec{SEC LinearPDE}. 136 $A$, $D$ and $Y$ are the known coeffecients of the PDE \index{partial differential equation!coefficients}. 137 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 By inpecting the Helmholtz equation \index{Helmholtz equation} 149 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 $\delta\hackscore{ij}$ is the Kronecker symbol \index{Kronecker symbol} defined by $\delta=\hackscore{ij}=1$ for 158 $i=j$ and $0$ otherwise. 159 160 We want to implement a 161 new class which we will call \class{Helmholtz} that provides the same methods like the \LinearPDE class but 162 is defined through the coefficeints $\kappa$, $\omega$, $f$, $\eta$, 163 $g$ rather than the general form given by \eqn{EQU.FEM.1}. 164 Python's 165 mechanism of inhertence allows doing this in a very easy way. 166 The advantage is that our new \class{Helmholtz} can be used in any context 167 that works with a \LinearPDE class but with an easier interface to define the PDE. 168 This improves reuasablity as well as maintainability of program codes. 169 170 We want to implement a 171 new class which we will call \class{Helmholtz} that provides the same methods like the \LinearPDE class but 172 is defined through the coefficeints $\kappa$, $\omega$, $f$, $\alpha$, 173 $g$ rather than the general form given by \eqn{EQU.FEM.1}. 174 Python's mechanism of subclasses allows doing this in a very easy way. 175 The \Poisson class of the \linearPDEsPack module, 176 which we have already used in \Chap{FirstSteps}, is in fact a subclass of the 177 \LinearPDE class. That means that it all methods (such as the \method{getSolution}) 178 from the parent class \LinearPDE are defined for any \Poisson object. However, new 179 methods can be added and methods of the parent class can be redefined. In fact, 180 the \Poisson class redefines the \method{setValue} of the \LinearPDE class which is called to assign 181 values to the coefficients of the PDE. This is exactly what we will do when we define 182 our new \class{Helmholtz} class: 183 \begin{python} 184 from esys.linearPDEs import LinearPDE 185 import numarray 186 class Helmholtz(LinearPDE) 187 def setValue(self,kappa=0,omega=1,f=0,eta=0,g=0) 188 self._setValue(A=kappa*numarray.identity(self.getDim()),D=omega,Y=f,d=eta,y=g) 189 \end{python} 190 \code{class Helmholtz(linearPDE)} declares the new \class{Helmholtz} class as a subclass 191 of the \LinearPDE which have imported in the first line of the script. 192 We add the method \method{setValue} to the class which overwrites the 193 \method{setValue} method of the \LinearPDE class. The new methods which has the 194 parameters of the Helmholtz \eqn{DIFFUSION HELM EQ 1} as arguments 195 maps the parameters of the coefficients of the general PDE defined 196 in \eqn{EQU.FEM.1}. The coefficient \var{A} is defined by Kroneckers symbol. we use the 197 \numarray function \function{identity} which return a square matrix which has ones on the 198 main diagonal and zeros out side the main diagonal. The argument of \function{identity} gives the order of the matrix. 199 Here we use 200 the \method{getDim} of the \LinearPDE class object \var{self} to get the spatial dimension of the domain of the 201 PDE. As we will make use of the \class{Helmholtz} class several times, it is convient to 202 put is definition into a file which we name \file{mytools.py} available in the \ExampleDirectory. 203 You can use your favourite editor to create and edit the file. 204 205 An object of the \class{Helmholtz} class is created through the statments: 206 \begin{python} 207 from mytools import * 208 mypde=Helmholtz(mydomain) 209 mypde.setValue(kappa=10.,omega=0.1,f=12) 210 u=mypde.getSolution() 211 \end{python} 212 In the first statement we import all definition from the \file{mytools.py} \index{scripts!\file{mytools.py}}. Make sure 213 that \file{mytools.py} is in directory from where you have started Python. 214 \var{mydomain} is the \Domain of the PDE. In the third statment values are 215 assigned to the PDE parameters. As no values for arguments \var{eta} and \var{g} are 216 specified the default values $0$ are used \footnote{It would be better to use the default value 217 \var{escript.Data()} rather then $0$ as then the coefficient would be defined as being not present and 218 would not be processed when the PDE is evaluated.}. In the forth statement the solution of the 219 PDE is returned. 220 221 We want to test our \class{Helmholtz} class on a rectangular domain 222 of length $l$ and height $h$. We do this by choosing a simple test solution, 223 here we take $u=x\hackscore{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 take $\kappa=1$ 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}=n\hackscore{0}$. 227 So we have to set $g=n\hackscore{0}+\eta x\hackscore{0}$. The following script \file{helmholtztest.py} 228 \index{scripts!\file{helmholtztest.py}} which is available in the \ExampleDirectory 229 implements this test problem using the \finley PDE solver: 230 \begin{python} 231 from mytools import * 232 from esys.escript import * 233 import esys.finley 234 #... set some parameters ... 235 omega=0.1 236 eta=10. 237 #... generate domain ... 238 mydomain = esys.finley.Rectangle(l0=5.,l1=1.,n0=50, n1=10) 239 #... open PDE and set coefficients ... 240 mypde=Helmholtz(mydomain) 241 n=mydomain.getNormal() 242 x=mydomain.getX() 243 mypde.setValue(1,omega,omega*x[0],eta,n[0]+eta*x[0]) 244 #... calculate error of the PDE solution ... 245 u=mypde.getSolution() 246 print "error is ",Lsup(u-x[0]) 247 \end{python} 248 The script is similar to the script \file{mypoisson.py} dicussed in \Chap{FirstSteps}. 249 \code{mydomain.getNormal()} returns the outer normal field on the surface of the domain. The function \function{Lsup} 250 is imported by the \code{from escript import *} statement and returns the maximum absulute value of it argument. To run 251 The error shown by the print statement should be in the order of $10^{-7}$. As piecewise linear interpolation is 252 used to approximate the solution and our solution is a linear function of the spatial coordinates one may 253 expect that the error is zero. However, as most PDE packages uses an iterative solver which is terminated 254 when a given tolerance has been reached. The default tolerance is $10^{-8}$. Thsi value can be altered by using the 255 \method{setTolerance} of the \LinearPDE class. 256 257 \section{The Transition Problem} 258 \label{DIFFUSION TRANS SEC} 259 Now we are ready to solve the original time dependent problem. The main 260 part of the script is the loop over the time $t$ which takes the following form: 261 \begin{python} 262 mypde=Helmholtz(mydomain) 263 while t

## Properties

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