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

Revision 569 - (show annotations)
Tue Feb 28 05:34:37 2006 UTC (14 years, 7 months ago) by gross
File MIME type: application/x-tex
File size: 17630 byte(s)
modified presentation of LinearPDE

 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 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$ 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$ 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(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-1)}\approx T^{(n)}+T\hackscore{,t}^{(n)}(t^{(n-1)}-t^{(n)}) 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 + \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 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 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 u\hackscore{,i} n\hackscore{i} = g - \eta u\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 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 $\delta\hackscore{ij}$ is the Kronecker symbol \index{Kronecker symbol} defined by $\delta\hackscore{ij}=1$ for 160 $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 it is more efficient to leave a coefficient undefined insted assigning zero to it.} 163 164 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 the PDE and returns the solution \code{u} of the PDE. \function{kronecker} is \escript function 177 returning the Kronecker symbol. 178 179 The coefficients can set by several calls of \method{setValue} where the order can be chosen arbitrarily. 180 If a value is assigned to a coefficint several times, the last assigned value is used when 181 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 The \LinearPDE class provides the method \method{checkSymmetry} method to check if the given PDE is symmetric. 197 \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 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 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 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 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 an easy calculation shows that we have to choose $f=\omega \cdot x\hackscore{0}$. On the boundary we get 221 $\kappa n\hackscore{i} u\hackscore{,i}=\kappa n\hackscore{0}$. 222 So we have to set $g=\kappa n\hackscore{0}+\eta x\hackscore{0}$. The following script \file{helmholtz.py} 223 \index{scripts!\file{helmholtz.py}} which is available in the \ExampleDirectory 224 implements this test problem using the \finley PDE solver: 225 \begin{python} 226 from esys.escript import * 227 from esys.escript.linearPDEs import LinearPDE 228 from esys.finley import Rectangle 229 #... set some parameters ... 230 kappa=1. 231 omega=0.1 232 eta=10. 233 #... generate domain ... 234 mydomain = Rectangle(l0=5.,l1=1.,n0=50, n1=10) 235 #... open PDE and set coefficients ... 236 mypde=LinearPDE(mydomain) 237 mypde.setSymmetryOn() 238 n=mydomain.getNormal() 239 x=mydomain.getX() 240 mypde.setValue(A=kappa*kronecker(mydomain),D=omega,Y=omega*x[0], \ 241 d=eta,y=kappa*n[0]+eta*x[0]) 242 #... calculate error of the PDE solution ... 243 u=mypde.getSolution() 244 print "error is ",Lsup(u-x[0]) 245 \end{python} 246 The script is similar to the script \file{poisson.py} dicussed in \Chap{FirstSteps}. 247 \code{mydomain.getNormal()} returns the outer normal field on the surface of the domain. The function \function{Lsup} 248 imported by the \code{from esys.escript import *} statement and returns the maximum absulute value of its argument. 249 The error shown by the print statement should be in the order of $10^{-7}$. As piecewise bi-linear interpolation is 250 used by \finley approximate the solution and our solution is a linear function of the spatial coordinates one might 251 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 \method{setTolerance} of the \LinearPDE class. 255 256 \subsection{The Transition Problem} 257 \label{DIFFUSION TRANS SEC} 258 Now we are ready to solve the original time dependent problem. The main 259 part of the script is the loop over time $t$ which takes the following form: 260 \begin{python} 261 t=0 262 T=Tref 263 mypde=LinearPDE(mydomain) 264 mypde.setValue(A=kappa*kronecker(mydomain),D=rhocp/h,d=eta,y=eta*Tref) 265 while t

## Properties

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