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

Revision 6928 - (hide annotations)
Mon Dec 9 05:29:21 2019 UTC (2 years, 5 months ago) by acodd
File MIME type: application/x-tex
File size: 18457 byte(s)
more updates to userguide


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

## Properties

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