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

Revision 3306 - (hide annotations)
Mon Oct 25 05:09:13 2010 UTC (10 years, 11 months ago) by caltinay
File MIME type: application/x-tex
File size: 18295 byte(s)
Commented declaremodule and modulesynopsis.


 1 ksteube 1811 2 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 3 ksteube 1316 % 4 jfenwick 2881 % Copyright (c) 2003-2010 by University of Queensland 5 ksteube 1811 % 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 caltinay 3279 \centerline{\includegraphics[width=\figwidth]{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 caltinay 3279 In this chapter we will discuss how to solve a time-dependent temperature 26 diffusion\index{diffusion equation} PDE for a given block of material. 27 Within the block there is a heat source which drives the temperature diffusion. 28 jgs 107 On the surface, energy can radiate into the surrounding environment. 29 jgs 102 \fig{DIFFUSION FIG 1} shows the configuration. 30 31 caltinay 3279 In the next \Sec{DIFFUSION TEMP SEC} we will present the relevant model. 32 A time integration scheme is introduced to calculate the temperature at given time nodes $t^{(n)}$. 33 We will see that at each time step a Helmholtz equation\index{Helmholtz equation} must be solved. 34 jgs 107 The implementation of a Helmholtz equation solver will be discussed in \Sec{DIFFUSION HELM SEC}. 35 caltinay 3279 In Section~\ref{DIFFUSION TRANS SEC} this solver is used to build a solver for 36 the temperature diffusion problem. 37 jgs 102 38 jgs 121 \subsection{\label{DIFFUSION TEMP SEC}Temperature Diffusion} 39 caltinay 3279 The unknown temperature $T$ is a function of its location in the domain and time $t>0$. 40 The governing equation in the interior of the domain is given by 41 jgs 102 \begin{equation} 42 jfenwick 3295 \rho c_p T_{,t} - (\kappa T_{,i})_{,i} = q_H 43 jgs 102 \label{DIFFUSION TEMP EQ 1} 44 \end{equation} 45 jfenwick 3295 where $\rho c_p$ and $\kappa$ are given material constants. 46 caltinay 3279 In case of a composite material the parameters depend on their location in the domain. 47 jfenwick 3295 $q_H$ is a heat source (or sink) within the domain. 48 caltinay 3279 We are using the Einstein summation convention\index{summation convention} as introduced in \Chap{FirstSteps}. 49 jfenwick 3295 In our case we assume $q_H$ to be equal to a constant heat 50 caltinay 3279 production rate $q^{c}$ on a circle or sphere with center $x^c$ and radius $r$, and $0$ elsewhere: 51 jgs 102 \begin{equation} 52 jfenwick 3295 q_H(x,t)=\left\{ 53 caltinay 3279 \begin{array}{l l} 54 q^c & \quad\text{if $\|x-x^c\| \le r$}\\ 55 0 & \quad\text{else}\\ 56 \end{array} \right. 57 jgs 102 \label{DIFFUSION TEMP EQ 1b} 58 \end{equation} 59 caltinay 3279 for all $x$ in the domain and time $t>0$. 60 jgs 102 61 caltinay 3279 On the surface of the domain we specify a radiation condition which prescribes 62 jfenwick 3295 the normal component of the flux $\kappa T_{,i}$ to be proportional 63 to the difference of the current temperature to the surrounding temperature $T_{ref}$: 64 jgs 102 \begin{equation} 65 jfenwick 3295 \kappa T_{,i} n_i = \eta (T_{ref}-T) 66 jgs 102 \label{DIFFUSION TEMP EQ 2} 67 \end{equation} 68 caltinay 3279 $\eta$ is a given material coefficient depending on the material of the block and the surrounding medium. 69 jfenwick 3295 $n_i$ is the $i$-th component of the outer normal field\index{outer normal field} at the surface of the domain. 70 jgs 102 71 caltinay 3279 To solve the time-dependent \eqn{DIFFUSION TEMP EQ 1} the initial temperature at time $t=0$ has to be given. 72 Here we assume that the initial temperature is the surrounding temperature: 73 jgs 102 \begin{equation} 74 jfenwick 3295 T(x,0)=T_{ref} 75 jgs 102 \label{DIFFUSION TEMP EQ 4} 76 \end{equation} 77 caltinay 3279 for all $x$ in the domain. Note that the initial conditions satisfy the 78 jgs 102 boundary condition defined by \eqn{DIFFUSION TEMP EQ 2}. 79 80 jgs 107 The temperature is calculated at discrete time nodes $t^{(n)}$ where 81 caltinay 3279 $t^{(0)}=0$ and $t^{(n)}=t^{(n-1)}+h$, where $h>0$ is the step size which is assumed to be constant. 82 In the following, the upper index ${(n)}$ refers to a value at time $t^{(n)}$. 83 The simplest and most robust scheme to approximate the time derivative of the 84 temperature is the backward Euler\index{backward Euler} scheme. 85 The backward Euler scheme is based on the Taylor expansion of $T$ at time $t^{(n)}$: 86 jgs 102 \begin{equation} 87 jfenwick 3295 T^{(n)}\approx T^{(n-1)}+T_{,t}^{(n)}(t^{(n)}-t^{(n-1)}) 88 =T^{(n-1)} + h \cdot T_{,t}^{(n)} 89 jgs 102 \label{DIFFUSION TEMP EQ 6} 90 \end{equation} 91 caltinay 3279 This is inserted into \eqn{DIFFUSION TEMP EQ 1}. By separating the terms at 92 $t^{(n)}$ and $t^{(n-1)}$ one gets for $n=1,2,3\ldots$ 93 jgs 102 \begin{equation} 94 jfenwick 3295 \frac{\rho c_p}{h} T^{(n)} - (\kappa T^{(n)}_{,i})_{,i} = q_H + \frac{\rho c_p}{h} T^{(n-1)} 95 jgs 102 \label{DIFFUSION TEMP EQ 7} 96 \end{equation} 97 jfenwick 3295 where $T^{(0)}=T_{ref}$ is taken form the initial condition given by \eqn{DIFFUSION TEMP EQ 4}. 98 jgs 107 Together with the natural boundary condition 99 \begin{equation} 100 jfenwick 3295 \kappa T_{,i}^{(n)} n_i = \eta (T_{ref}-T^{(n)}) 101 jgs 107 \label{DIFFUSION TEMP EQ 2222} 102 \end{equation} 103 caltinay 3279 taken from \eqn{DIFFUSION TEMP EQ 2} this forms a boundary value problem that 104 has to be solved for each time step. 105 As a first step to implement a solver for the temperature diffusion problem we 106 will implement a solver for the boundary value problem that has to be solved 107 at each time step. 108 jgs 102 109 jgs 121 \subsection{\label{DIFFUSION HELM SEC}Helmholtz Problem} 110 jgs 102 The partial differential equation to be solved for $T^{(n)}$ has the form 111 \begin{equation} 112 jfenwick 3295 \omega T^{(n)} - (\kappa T^{(n)}_{,i})_{,i} = f 113 jgs 102 \label{DIFFUSION HELM EQ 1} 114 \end{equation} 115 lkettle 575 and we set 116 jgs 102 \begin{equation} 117 jfenwick 3295 \omega=\frac{\rho c_p}{h} \mbox{ and } f=q_H +\frac{\rho c_p}{h}T^{(n-1)} \;. 118 jgs 102 \label{DIFFUSION HELM EQ 1b} 119 \end{equation} 120 jfenwick 3295 With $g=\eta T_{ref}$ the radiation condition defined by \eqn{DIFFUSION TEMP EQ 2222} takes the form 121 jgs 102 \begin{equation} 122 jfenwick 3295 \kappa T^{(n)}_{,i} n_{i} = g - \eta T^{(n)}\text{ on $\Gamma$} 123 jgs 102 \label{DIFFUSION HELM EQ 2} 124 \end{equation} 125 caltinay 3279 The partial differential \eqn{DIFFUSION HELM EQ 1} together with boundary 126 conditions of \eqn{DIFFUSION HELM EQ 2} is called the Helmholtz equation\index{Helmholtz equation}. 127 jgs 102 128 caltinay 3279 We want to use the \LinearPDE class provided by \escript to define and solve a 129 general linear, steady, second order PDE such as the Helmholtz equation. 130 For a single PDE the \LinearPDE class supports the following form: 131 gross 625 \begin{equation}\label{LINEARPDE.SINGLE.1 TUTORIAL} 132 jfenwick 3295 -(A_{jl} u_{,l})_{,j}+D u = Y 133 jgs 102 \end{equation} 134 caltinay 3279 where we show only the coefficients relevant for the problem discussed here. 135 For the general form of a single PDE see \eqn{LINEARPDE.SINGLE.1}. 136 The coefficients $A$ and $Y$ have to be specified through \Data objects in 137 the \Function on the PDE or objects that can be converted into such \Data objects. 138 $A$ is a \RankTwo and $D$ and $Y$ are scalar. 139 The following natural boundary conditions\index{boundary condition!natural} 140 are considered on $\Gamma$: 141 gross 625 \begin{equation}\label{LINEARPDE.SINGLE.2 TUTORIAL} 142 jfenwick 3295 n_{j}A_{jl} u_{,l}+d u= y \;. 143 jgs 102 \end{equation} 144 caltinay 3279 Notice that the coefficient $A$ is the same like in the PDE~\eqn{LINEARPDE.SINGLE.1 TUTORIAL}. 145 The coefficients $d$ and $y$ are each a \Scalar in the \FunctionOnBoundary. 146 Constraints\index{constraint} for the solution prescribe the value of the 147 gross 568 solution at certain locations in the domain. They have the form 148 gross 625 \begin{equation}\label{LINEARPDE.SINGLE.3 TUTORIAL} 149 gross 568 u=r \mbox{ where } q>0 150 \end{equation} 151 caltinay 3279 $r$ and $q$ are each \Scalar where $q$ is the characteristic function\index{characteristic function} defining where the constraint is applied. 152 The constraints defined by \eqn{LINEARPDE.SINGLE.3 TUTORIAL} override any 153 other condition set by \eqn{LINEARPDE.SINGLE.1 TUTORIAL} or \eqn{LINEARPDE.SINGLE.2 TUTORIAL}. 154 The \Poisson class of the \linearPDEs module, which we have already used in 155 \Chap{FirstSteps}, is in fact a subclass of the more general \LinearPDE class. 156 The \linearPDEs module provides a \Helmholtz class but we will make direct use 157 of the general \LinearPDE class. 158 jgs 102 159 caltinay 3279 By inspecting the Helmholtz equation\index{Helmholtz equation} 160 (\ref{DIFFUSION HELM EQ 1}) and boundary condition (\ref{DIFFUSION HELM EQ 2}), 161 and substituting $u$ for $T^{(n)}$, we can easily assign values to the 162 coefficients in the general PDE of the \LinearPDE class: 163 jgs 102 \begin{equation}\label{DIFFUSION HELM EQ 3} 164 \begin{array}{llllll} 165 jfenwick 3295 A_{ij}=\kappa \delta_{ij} & D=\omega & Y=f \\ 166 jgs 102 d=\eta & y= g & \\ 167 \end{array} 168 \end{equation} 169 jfenwick 3295 $\delta_{ij}$ is the Kronecker symbol\index{Kronecker symbol} defined 170 by $\delta_{ij}=1$ for $i=j$ and $0$ otherwise. 171 caltinay 3279 Undefined coefficients are assumed to be not present.\footnote{There is a 172 difference in \escript for a coefficient to be not present and set to zero. 173 Since in the former case the coefficient is not processed, it is more efficient 174 to leave it undefined instead of assigning zero to it.} 175 In this diffusion example we do not need to define a characteristic function 176 $q$ because the boundary conditions we consider in \eqn{DIFFUSION HELM EQ 2} 177 are just the natural boundary conditions which are already defined in the 178 \LinearPDE class (shown in \eqn{LINEARPDE.SINGLE.2 TUTORIAL}). 179 jgs 102 180 caltinay 3279 The Helmholtz equation can be set up the following way\footnote{Note that this 181 is not a complete code. The full source code can be found in helmholtz.py''.}: 182 gross 568 \begin{python} 183 caltinay 3279 mypde=LinearPDE(mydomain) 184 mypde.setValue(A=kappa*kronecker(mydomain),D=omega,Y=f,d=eta,y=g) 185 u=mypde.getSolution() 186 gross 568 \end{python} 187 where we assume that \code{mydomain} is a \Domain object and 188 caltinay 3279 \code{kappa}, \code{omega}, \code{eta}, and \code{g} are given scalar values 189 typically \code{float} or \Data objects. The \method{setValue} method 190 assigns values to the coefficients of the general PDE. 191 The \method{getSolution} method solves the PDE and returns the solution 192 \code{u} of the PDE. \function{kronecker} is an \escript function returning 193 the Kronecker symbol. 194 gross 568 195 caltinay 3279 The coefficients can set by several calls to \method{setValue} where the order 196 can be chosen arbitrarily. 197 If a value is assigned to a coefficient several times, the last assigned value 198 is used when the solution is calculated: 199 gross 568 \begin{python} 200 caltinay 3279 mypde = LinearPDE(mydomain) 201 mypde.setValue(A=kappa*kronecker(mydomain), d=eta) 202 mypde.setValue(D=omega, Y=f, y=g) 203 mypde.setValue(d=2*eta) # overwrites d=eta 204 u=mypde.getSolution() 205 gross 568 \end{python} 206 caltinay 3279 In some cases the solver of the PDE can make use of the fact that the PDE is 207 symmetric\index{symmetric PDE}. A PDE is called symmetric if 208 gross 625 \begin{equation}\label{LINEARPDE.SINGLE.4 TUTORIAL} 209 jfenwick 3295 A_{jl}=A_{lj}\;. 210 gross 568 \end{equation} 211 caltinay 3279 Note that $D$ and $d$ may have any value and the right hand sides $Y$, $y$ as 212 well as the constraints are not relevant. The Helmholtz problem is symmetric. 213 The \LinearPDE class provides the method \method{checkSymmetry} to check if 214 the given PDE is symmetric. 215 gross 568 \begin{python} 216 caltinay 3279 mypde = LinearPDE(mydomain) 217 mypde.setValue(A=kappa*kronecker(mydomain), d=eta) 218 print mypde.checkSymmetry() # returns True 219 mypde.setValue(B=kronecker(mydomain)) 220 print mypde.checkSymmetry() # returns False 221 mypde.setValue(C=kronecker(mydomain)) 222 print mypde.checkSymmetry() # returns True 223 gross 568 \end{python} 224 caltinay 3279 Unfortunately, calling \method{checkSymmetry} is very expensive and is only 225 recommended for testing and debugging purposes. 226 The \method{setSymmetryOn} method is used to declare a PDE symmetric: 227 gross 568 \begin{python} 228 caltinay 3279 mypde = LinearPDE(mydomain) 229 mypde.setValue(A=kappa*kronecker(mydomain)) 230 mypde.setSymmetryOn() 231 gross 568 \end{python} 232 caltinay 3279 Now we want to see how we actually solve the Helmholtz equation on a 233 jfenwick 3295 rectangular domain of length $l_{0}=5$ and height $l_{1}=1$. 234 caltinay 3279 We choose a simple test solution such that we can verify the returned solution 235 jfenwick 3295 against the exact answer. Actually, we take $T=x_{0}$ (here 236 $q_H = 0$) and then calculate the right hand side terms $f$ and $g$ 237 caltinay 3279 such that the test solution becomes the solution of the problem. 238 If we assume $\kappa$ as being constant, an easy calculation shows that we 239 jfenwick 3295 have to choose $f=\omega \cdot x_{0}$. 240 On the boundary we get $\kappa n_{i} u_{,i}=\kappa n_{0}$. 241 Thus we have to set $g=\kappa n_{0}+\eta x_{0}$. 242 caltinay 3279 The following script \file{helmholtz.py}\index{scripts!\file{helmholtz.py}} 243 which is available in the \ExampleDirectory implements this test problem using 244 the \finley PDE solver: 245 jgs 102 \begin{python} 246 caltinay 3279 from esys.escript import * 247 from esys.escript.linearPDEs import LinearPDE 248 from esys.finley import Rectangle 249 # set some parameters 250 kappa=1. 251 omega=0.1 252 eta=10. 253 # generate domain 254 mydomain = Rectangle(l0=5., l1=1., n0=50, n1=10) 255 # open PDE and set coefficients 256 mypde=LinearPDE(mydomain) 257 mypde.setSymmetryOn() 258 n=mydomain.getNormal() 259 x=mydomain.getX() 260 mypde.setValue(A=kappa*kronecker(mydomain), D=omega,Y=omega*x, \ 261 d=eta, y=kappa*n+eta*x) 262 # calculate error of the PDE solution 263 u=mypde.getSolution() 264 print("error is ",Lsup(u-x)) 265 saveVTK("x0.vtu", sol=u) 266 jgs 102 \end{python} 267 caltinay 3279 To visualize the solution `x0.vtu' you can use the command 268 \begin{verbatim} 269 mayavi -d x0.vtu -m SurfaceMap 270 \end{verbatim} 271 jfenwick 3295 and it is easy to see that the solution $T=x_{0}$ is calculated. 272 lkettle 575 273 lgraham 1700 The script is similar to the script \file{poisson.py} discussed in \Chap{FirstSteps}. 274 caltinay 3279 \code{mydomain.getNormal()} returns the outer normal field on the surface of the domain. 275 The function \function{Lsup} is imported by the \code{from esys.escript import *} 276 statement and returns the maximum absolute value of its argument. 277 The error shown by the print statement should be in the order of $10^{-7}$. 278 As piecewise bi-linear interpolation is used by \finley to approximate the 279 solution, and our solution is a linear function of the spatial coordinates, 280 one might expect that the error would be zero or in the order of machine 281 precision (typically $\approx 10^{-15}$). 282 However most PDE packages use an iterative solver which is terminated when a 283 given tolerance has been reached. The default tolerance is $10^{-8}$. 284 This value can be altered by using the \method{setTolerance} of the \LinearPDE class. 285 jgs 102 286 jgs 121 \subsection{The Transition Problem} 287 jgs 102 \label{DIFFUSION TRANS SEC} 288 ksteube 1316 Now we are ready to solve the original time-dependent problem. The main 289 jgs 107 part of the script is the loop over time $t$ which takes the following form: 290 jgs 102 \begin{python} 291 caltinay 3279 t=0 292 T=Tref 293 mypde=LinearPDE(mydomain) 294 mypde.setValue(A=kappa*kronecker(mydomain), D=rhocp/h, d=eta, y=eta*Tref) 295 while t

## Properties

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