# Diff of /trunk/doc/cookbook/onedheatdiff001.tex

revision 2860 by ahallam, Thu Dec 3 01:45:48 2009 UTC revision 2861 by gross, Wed Jan 20 02:47:42 2010 UTC
# Line 29  which is defined as: Line 29  which is defined as:
29  \rho c\hackscore p \frac{\partial T}{\partial t} - \kappa \frac{\partial^{2} T}{\partial x^{2}} = q\hackscore H  \rho c\hackscore p \frac{\partial T}{\partial t} - \kappa \frac{\partial^{2} T}{\partial x^{2}} = q\hackscore H
30  \label{eqn:hd}  \label{eqn:hd}
31
32  where $\rho$ is the material density, $c\hackscore p$ is the specific heat and $\kappa$ is the thermal conductivity constant for a given material\footnote{A list of some common thermal conductivities is available from Wikipedia \url{http://en.wikipedia.org/wiki/List_of_thermal_conductivities}}.  where $\rho$ is the material density, $c\hackscore p$ is the specific heat and $\kappa$ is the thermal
33  The heat source is defined by the right hand side of \ref{eqn:hd} as $q\hackscore{H}$; this can take the form of a constant or a function of time and space. For example $q\hackscore{H} = Te^{-\gamma t}$ where we have the output of our heat source decaying with time. There are also two partial derivatives in \ref{eqn:hd}; $\frac{\partial T}{\partial t}$ describes the change in temperature with time while $\frac{\partial ^2 T}{\partial x^2}$ is the spatial change of temperature. As there is only a single spatial dimension to our problem, our temperature solution $T$ is only dependent on the time $t$ and our position along the iron bar $x$.  conductivity\footnote{A list of some common thermal conductivities is available from Wikipedia \url{http://en.wikipedia.org/wiki/List_of_thermal_conductivities}}. Here we assume that these material
34    parameters are \textbf{constant}.
35    The heat source is defined by the right hand side of \refEq{eqn:hd} as $q\hackscore{H}$; this can take the form of a constant or a function of time and space. For example $q\hackscore{H} = q\hackscore{0}e^{-\gamma t}$ where we have the output of our heat source decaying with time. There are also two partial derivatives in \refEq{eqn:hd}; $\frac{\partial T}{\partial t}$ describes the change in temperature with time while $\frac{\partial ^2 T}{\partial x^2}$ is the spatial change of temperature. As there is only a single spatial dimension to our problem, our temperature solution $T$ is only dependent on the time $t$ and our position along the iron bar $x$.
36
37    \subsection{PDEs and the General Form}
38    Potentially, it is now possible to solve PDE \refEq{eqn:hd} analytically and this would produce an exact solution to our problem. However, it is not always possible or practical to solve a problem this way. Alternatively, computers can be used to solve these kinds of problems. To do this, a numerical approach is required to discretised
39    the PDE \refEq{eqn:hd} in time and space so finally we are left with a finite number of equations for a finite number of spatial and time steps in the model. While discretization introduces approximations and a degree of error, we find that a sufficiently sampled model is generally accurate enough for the requirements of the modeler.
40
41    Firstly, we will discretise the PDE \refEq{eqn:hd} in the time direction which will
42    leave as with a steady linear PDE which is involving spatial derivatives only and needs to be solved in each time
43    step to progress in time - \esc can help us here.
44
45  \subsection{\esc, PDEs and The General Form}  For the discretization in time we will use is the Backwards Euler approximation scheme\footnote{see \url{http://en.wikipedia.org/wiki/Euler_method}}. It bases on the
46  Potentially, it is now possible to solve \ref{eqn:hd} analytically and this would produce an exact solution to our problem. However, it is not always possible or practical to solve a problem this way. Alternatively, computers can be used to solve these kinds of problems when a large number of sums or a more complex visualisation is required. To do this, a numerical approach is required - \esc can help us here -  and it becomes necessary to discretize the equation so that we are left with a finite number of equations for a finite number of spatial and time steps in the model. While discretization introduces approximations and a degree of error, we find that a sufficiently sampled model is generally accurate enough for the requirements of the modeller.  approximation
47
48    \frac{\partial T(t)}{\partial t} \approx \frac{T(t)-T(t-h)}{h}
49    \label{eqn:beuler}
50
51    for  $\frac{\partial T}{\partial t}$  at time $t$
52    where $h$ is the time step size. This can also be written as;
53
54    \frac{\partial T}{\partial t}(t^{(n)}) \approx \frac{T^{(n)} - T^{(n-1)}}{h}
55    \label{eqn:Tbeuler}
56
57    where the upper index $n$ denotes the n\textsuperscript{th} time step. So one has
58
59    \begin{array}{rcl}
60    t^{(n)} & = & t^{(n-1)}+h \\
61    T^{(n)} & = & T(t^{(n-1)}) \\
62    \end{array}
63    \label{eqn:Neuler}
64
65    Substituting \refEq{eqn:Tbeuler} into \refEq{eqn:hd} we get;
66
67    \frac{\rho c\hackscore p}{h} (T^{(n)} - T^{(n-1)}) - \kappa \frac{\partial^{2} T^{(n)}}{\partial x^{2}} = q\hackscore H
68    \label{eqn:hddisc}
69
70    Notice that we evaluate the spatial derivative term at current time $t^{(n)}$ - therefore the name \textbf{backward Euler} scheme. Alternatively, one can use evaluate the spatial derivative term at the previous time $t^{(n-1)}$. This
71    approach is called the \textbf{forward Euler} scheme. This scheme can provide some computational advantages which
72    we are not discussed here but has the major disadvantage that depending on the
73    material parameter as well as the discretiztion of the spatial derivative term the time step size $h$ needs to be chosen sufficiently small to achieve a stable temperature when progressing in time. The term \textit{stable} means
74    that the approximation of the temperature will not grow beyond its initial bounds and becomes unphysical.
75    The backward Euler which we use here is unconditionally stable meaning that under the assumption of
76    physically correct problem set-up the temperature approximation remains physical for all times.
77    The user needs to keep in mind that the discretization error introduced by \refEq{eqn:beuler}
78    is sufficiently small so a good approximation of the true temperature is calculated. It is
79    therefore crucial that the user remains critical about his/her results and for instance compares
80    the results for different time and spatial step sizes.
81
82    To get the temperature $T^{(n)}$ at time $t^{(n)}$ we need to solve the linear
83    differential equation \refEq{eqn:hddisc} which is only including spatial derivatives. To solve this problem
84    we want to to use \esc.
85
86  \esc interfaces with any given PDE via a general form. In this example we will illustrate a simpler version of the full linear PDE general form which is available in the \esc user's guide. A simplified form that suits our heat diffusion problem\footnote{In the form of the \esc users guide which using the Einstein convention is written as  \esc interfaces with any given PDE via a general form. For the purpose of this introduction we will illustrate a simpler version of the full linear PDE general form which is available in the \esc user's guide. A simplified form that suits our heat diffusion problem\footnote{In the form of the \esc users guide which using the Einstein convention is written as
87  $-(A\hackscore{jl} u\hackscore{,l})\hackscore{,j}+D u =Y$}  $-(A\hackscore{jl} u\hackscore{,l})\hackscore{,j}+D u =Y$}
88  is described by;  is described by;
89  \label{eqn:commonform nabla}  \label{eqn:commonform nabla}
90  -\nabla\cdot(A\cdot\nabla u) + Du = f  -\nabla\cdot(A\cdot\nabla u) + Du = f
91
92  where $A$, $D$ and $f$ are known values. The symbol $\nabla$ which is called the \textit{Nabla operator} or \textit{del operator} represents  where $A$, $D$ and $f$ are known values and $u$ is the unknown solution. The symbol $\nabla$ which is called the \textit{Nabla operator} or \textit{del operator} represents
93  the spatial derivative of its subject - in this case $u$. Lets assume for a moment that we deal with a one-dimensional problem then ;  the spatial derivative of its subject - in this case $u$. Lets assume for a moment that we deal with a one-dimensional problem then ;
94
95  \nabla = \frac{\partial}{\partial x}  \nabla = \frac{\partial}{\partial x}
96
97  and we can write \ref{eqn:commonform nabla} as;  and we can write \refEq{eqn:commonform nabla} as;
98  \label{eqn:commonform}  \label{eqn:commonform}
99  -A\frac{\partial^{2}u}{\partial x^{2}} + Du = f  -A\frac{\partial^{2}u}{\partial x^{2}} + Du = f
100
101  if $A$ is constant then  \ref{eqn:commonform} is consistent with our heat diffusion problem in \ref{eqn:hd} with the exception of $u$. Thus when comparing equations \ref{eqn:hd} and \ref{eqn:commonform} we see that;  if $A$ is constant. To match this simplified general form to our problem \refEq{eqn:hddisc}
102    we rearrange \refEq{eqn:hddisc};
A = \kappa; D = \rho c \hackscore{p}; f = q \hackscore{H}

We can write the partial $\frac{\partial T}{\partial t}$ in terms of $u$ by discretising the time of our solution. The method we will use is the Backwards Euler approximation, which states;

\frac{\partial f(x)}{\partial x} \approx \frac{f(x+h)-f(x)}{h}
\label{eqn:beuler}

where h is the the discrete step size $\Delta x$.
Now if the temperature $T(t)$ is substituted in by letting $f(x) = T(t)$ then from \ref{eqn:beuler} we see that;

T'(t) \approx \frac{T(t+h) - T(t)}{h}

which can also be written as;

\frac{\partial T^{(n)}}{\partial t} \approx \frac{T^{(n)} - T^{(n-1)}}{h}
\label{eqn:Tbeuler}

where $n$ denotes the n\textsuperscript{th} time step. Substituting \ref{eqn:Tbeuler} into \ref{eqn:hd} we get;

\frac{\rho c\hackscore p}{h} (T^{(n)} - T^{(n-1)}) - \kappa \frac{\partial^{2} T}{\partial x^{2}} = q\hackscore H
\label{eqn:hddisc}

To fit our simplified general form we can rearrange \ref{eqn:hddisc};
103
104  \frac{\rho c\hackscore p}{h} T^{(n)} - \kappa \frac{\partial^2 T^{(n)}}{\partial x^2} = q\hackscore H +  \frac{\rho c\hackscore p}{h} T^{(n-1)}  \frac{\rho c\hackscore p}{h} T^{(n)} - \kappa \frac{\partial^2 T^{(n)}}{\partial x^2} = q\hackscore H +  \frac{\rho c\hackscore p}{h} T^{(n-1)}
105  \label{eqn:hdgenf}  \label{eqn:hdgenf}
106
107  The PDE is now in a form that satisfies \refEq{eqn:commonform nabla} which is required for \esc to solve our PDE. This can be done by generating a solution for successive increments in the time nodes $t^{(n)}$ where  The PDE is now in a form that satisfies \refEq{eqn:commonform nabla} which is required for \esc to solve our PDE. This can be done by generating a solution for successive increments in the time nodes $t^{(n)}$ where
108  $t^{(0)}=0$ and  $t^{(n)}=t^{(n-1)}+h$ where $h>0$ is the step size and assumed to be constant.  $t^{(0)}=0$ and  $t^{(n)}=t^{(n-1)}+h$ where $h>0$ is the step size and assumed to be constant.
109  In the following the upper index ${(n)}$ refers to a value at time $t^{(n)}$. Finally, by comparing \ref{eqn:hdgenf} with \ref{eqn:commonform} it can be seen that;  In the following the upper index ${(n)}$ refers to a value at time $t^{(n)}$. Finally, by comparing \refEq{eqn:hdgenf} with \refEq{eqn:commonform} it can be seen that;
110
111    u=T^{(n)};
112  A = \kappa; D = \frac{\rho c \hackscore{p}}{h}; f = q \hackscore{H} + \frac{\rho c\hackscore p}{h} T^{(n-1)}  A = \kappa; D = \frac{\rho c \hackscore{p}}{h}; f = q \hackscore{H} + \frac{\rho c\hackscore p}{h} T^{(n-1)}
113
114
115  Now that the general form has been established, it can be submitted to \esc. Note that it is necessary to establish the state of our system at time zero or $T^{(n=0)}$. This is due to the time derivative approximation we have used from \ref{eqn:Tbeuler}. Our model stipulates a starting temperature in the iron bar of 0\textcelsius. Thus the temperature distribution is simply;  Now that the general form has been established, it can be submitted to \esc. Note that it is necessary to establish the state of our system at time zero or $T^{(n=0)}$. This is due to the time derivative approximation we have used from \refEq{eqn:Tbeuler}. Our model stipulates a starting temperature in the iron bar of 0\textcelsius. Thus the temperature distribution is simply;
116
117  T(x,0) = T\hackscore{ref} = 0  T(x,0) = T\hackscore{ref} = 0
118
119  for all $x$ in the domain.  for all $x$ in the domain.
120
121  \subsection{Boundary Conditions}  \subsection{Boundary Conditions}
122  With the PDE sufficiently modified, consideration must now be given to the boundary conditions of our model. Typically there are two main types of boundary conditions known as Neumann and Dirichlet\footnote{More information on Boundary Conditions is available at Wikipedia \url{http://en.wikipedia.org/wiki/Boundary_conditions}}. In this example, we have utilised both conditions. Dirichlet is conceptually simpler and is used to prescribe a known value to the model on its boundary. For this model Dirichlet boundary conditions exist where we have applied our heat source. As the heat source is a constant, we can simulate its presence on that boundary. This is done by continuously resetting the temperature of the boundary, so that it is the same as the heat source.    With the PDE sufficiently modified, consideration must now be given to the boundary conditions of our model. Typically there are two main types of boundary conditions known as Neumann and Dirichlet boundary conditions\footnote{More information on Boundary Conditions is available at Wikipedia \url{http://en.wikipedia.org/wiki/Boundary_conditions}}, respectively.
123    A Dirichlet boundary condition is conceptually simpler and is used to prescribe a known value to the unknown - in our example the temperature - on boundary of the region of interest.
124    \editor{LUTZ: This is not correct!!}
125    For this model Dirichlet boundary conditions exist where we have applied our heat source. As the heat source is a constant, we can simulate its presence on that boundary. This is done by continuously resetting the temperature of the boundary, so that it is the same as the heat source.
126
127  Neumann boundary conditions describe the radiation or flux that is normal to the boundary surface. This aptly describes our insulation conditions as we do not want to exert a constant temperature as with the heat source. However, we do want to prevent any loss of energy from the system.  Neumann boundary conditions describe the radiation or flux that is normal to the boundary surface. This aptly describes our insulation conditions as we do not want to exert a constant temperature as with the heat source. However, we do want to prevent any loss of energy from the system.
128
# Line 111  and simplified to our one dimensional mo Line 139  and simplified to our one dimensional mo
139  where $\eta$ is a given material coefficient depending on the material of the block and the surrounding medium and $\hat{n}\hackscore i$ is the $i$-th component of the outer normal field \index{outer normal field} at the surface of the domain. These two conditions form a boundary value problem that has to be solved for each time step. Due to the perfect insulation in our model we can set $\eta = 0$ which results in zero flux - no energy in or out - we do not need to worry about the Neumann terms of the general form for this example.  where $\eta$ is a given material coefficient depending on the material of the block and the surrounding medium and $\hat{n}\hackscore i$ is the $i$-th component of the outer normal field \index{outer normal field} at the surface of the domain. These two conditions form a boundary value problem that has to be solved for each time step. Due to the perfect insulation in our model we can set $\eta = 0$ which results in zero flux - no energy in or out - we do not need to worry about the Neumann terms of the general form for this example.
140
141  \subsection{A \textit{1D} Clarification}  \subsection{A \textit{1D} Clarification}
142  It is necessary for clarification that we revisit the general PDE from \refeq{eqn:commonform nabla} under the light of a two dimensional domain. \esc is inherently designed to solve problems that are greater than one dimension and so \ref{eqn:commonform nabla} needs to be read as a higher dimensional problem. In the case of two spatial dimensions the \textit{Nabla operator} has in fact two components $\nabla = (\frac{\partial}{\partial x}, \frac{\partial}{\partial y})$. In full, \ref{eqn:commonform nabla} assuming a constant coefficient $A$, takes the form;  It is necessary for clarification that we revisit the general PDE from \refeq{eqn:commonform nabla} under the light of a two dimensional domain. \esc is inherently designed to solve problems that are greater than one dimension and so \refEq{eqn:commonform nabla} needs to be read as a higher dimensional problem. In the case of two spatial dimensions the \textit{Nabla operator} has in fact two components $\nabla = (\frac{\partial}{\partial x}, \frac{\partial}{\partial y})$. In full, \refEq{eqn:commonform nabla} assuming a constant coefficient $A$, takes the form;
143  \label{eqn:commonform2D}  \label{eqn:commonform2D}
144  -A\hackscore{00}\frac{\partial^{2}u}{\partial x^{2}}  -A\hackscore{00}\frac{\partial^{2}u}{\partial x^{2}}
145  -A\hackscore{01}\frac{\partial^{2}u}{\partial x\partial y}  -A\hackscore{01}\frac{\partial^{2}u}{\partial x\partial y}
# Line 122  It is necessary for clarification that w Line 150  It is necessary for clarification that w
150  Notice that for the higher dimensional case $A$ becomes a matrix. It is also  Notice that for the higher dimensional case $A$ becomes a matrix. It is also
151  important to notice that the usage of the Nabla operator creates  important to notice that the usage of the Nabla operator creates
152  a compact formulation which is also independent from the spatial dimension.  a compact formulation which is also independent from the spatial dimension.
153  So to make the general PDE~\ref{eqn:commonform2D} one dimensional as  So to make the general PDE \refEq{eqn:commonform2D} one dimensional as
154  shown in~\ref{eqn:commonform} we need to set  shown in \refEq{eqn:commonform} we need to set
155
156  A\hackscore{00}=A; A\hackscore{01}=A\hackscore{10}=A\hackscore{11}=0  A\hackscore{00}=A; A\hackscore{01}=A\hackscore{10}=A\hackscore{11}=0
157
158
159  \subsection{Developing a PDE Solution Script}  \subsection{Developing a PDE Solution Script}
160  \label{sec:key}  \label{sec:key}
161  To solve the heat diffusion equation (equation \ref{eqn:hd}) we will write a simple \pyt script which uses the \modescript, \modfinley and \modmpl modules. At this point we assume that you have some basic understanding of the \pyt programming language. If not there are some pointers and links available in Section \ref{sec:escpybas} .  To solve the heat diffusion equation (equation \refEq{eqn:hd}) we will write a simple \pyt script which uses the \modescript, \modfinley and \modmpl modules. At this point we assume that you have some basic understanding of the \pyt programming language. If not there are some pointers and links available in Section \ref{sec:escpybas} .
162
163  By developing a script for \esc, the heat diffusion equation can be solved at successive time steps for a predefined period using our general form \ref{eqn:hdgenf}. Firstly it is necessary to import all the libraries\footnote{The libraries contain predefined scripts that are required to solve certain problems, these can be simple like $sine$ and $cosine$ functions or more complicated like those from our \esc library.}  By developing a script for \esc, the heat diffusion equation can be solved at successive time steps for a predefined period using our general form \refEq{eqn:hdgenf}. Firstly it is necessary to import all the libraries\footnote{The libraries contain predefined scripts that are required to solve certain problems, these can be simple like $sine$ and $cosine$ functions or more complicated like those from our \esc library.}
164  that we will require.  that we will require.
165  \begin{python}  \begin{python}
166  from esys.escript import *  from esys.escript import *

Legend:
 Removed from v.2860 changed lines Added in v.2861