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

revision 106 by jgs, Fri Dec 17 07:43:12 2004 UTC revision 107 by jgs, Thu Jan 27 06:21:48 2005 UTC
# Line 1  Line 1
1  % $Id$  % $Id$
2  \chapter{How to Solve The Diffusion Equation}  \chapter{How to Solve A Diffusion Problem}
3  \label{DIFFUSION CHAP}  \label{DIFFUSION CHAP}
4
5  \begin{figure}  \begin{figure}
# Line 8  Line 8
8  \label{DIFFUSION FIG 1}  \label{DIFFUSION FIG 1}
9  \end{figure}  \end{figure}
10
\begin{figure}
\centerline{\includegraphics[width=\figwidth]{DiffusionRes1}}
\centerline{\includegraphics[width=\figwidth]{DiffusionRes16}}
\centerline{\includegraphics[width=\figwidth]{DiffusionRes32}}
\centerline{\includegraphics[width=\figwidth]{DiffusionRes48}}
\caption{Results of the Temperture Diffusion Problem for Time Steps $1$ $16$, $32$ and $48$.}
\label{DIFFUSION FIG 2}
\end{figure}

11  \section{\label{DIFFUSION OUT SEC}Outline}  \section{\label{DIFFUSION OUT SEC}Outline}
12  In this chapter we will discuss how to solve the time depeneded-temperature diffusion\index{diffusion equation} within  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.  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.  On the surface, energy can radiate into the surrounding environment.
15  \fig{DIFFUSION FIG 1} shows the configuration.  \fig{DIFFUSION FIG 1} shows the configuration.
16
17  In the next \Sec{DIFFUSION TEMP SEC} we will present the relevant model. A  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)}$.  time integration scheme is introduced to calculate the temperature at given time nodes $t^{(n)}$.
19  We will see that at time step a so-called Helmholtz equation \index{Helmholtz equation} has to be solved.  We will see that at each time step a Helmholtz equation \index{Helmholtz equation}
20  The implementation iof a Helmholtz equation solver will be discussed in \Sec{DIFFUSION HELM SEC}.  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  In Section~\ref{DIFFUSION TRANS SEC} the solver of the Helmholtz equation is used to build a
23  solver for the temperature diffusion problem.  solver for the temperature diffusion problem.
24
25  \section{\label{DIFFUSION TEMP SEC}Temperature Diffusion}  \section{\label{DIFFUSION TEMP SEC}Temperature Diffusion}
26
27  The temperature $T$ is a function of its location in the domain and time $t>0$. The governing equation  The unknown temperature $T$ is a function of its location in the domain and time $t>0$. The governing equation
28  in the interior of the domain is given by  in the interior of the domain is given by
29
30  \rho c\hackscore p T\hackscore{,t} - (\kappa T\hackscore{,i})\hackscore{,i} = q  \rho c\hackscore p T\hackscore{,t} - (\kappa T\hackscore{,i})\hackscore{,i} = q
31  \label{DIFFUSION TEMP EQ 1}  \label{DIFFUSION TEMP EQ 1}
32
33  where $\rho c\hackscore p$ and $\kappa$ are given material constants. In case of a composite  where $\rho c\hackscore p$ and $\kappa$ are given material constants. In case of a composite
34  material the parameters are depending on the their location in the domain. $q$ is  material the parameters depend on their location in the domain. $q$ is
35  a heat source (or sink) within the domain. We are using Einstein summation convention \index{summation convention}  a heat source (or sink) within the domain. We are using Einstein summation convention \index{summation convention}
36  as introduced in \Chap{FirstSteps}. In our case we assume $q$ to be equal to a constant $q^{c}$  as introduced in \Chap{FirstSteps}. In our case we assume $q$ to be equal to a constant heat production rate
37  on a circle or sphere with center $x^c$ and radius $r$ and $0$ elsewhere:  $q^{c}$ on a circle or sphere with center $x^c$ and radius $r$ and $0$ elsewhere:
38
39  q(x,t)=  q(x,t)=
40  \left\{  \left\{
# Line 58  q^c  & & \|x-x^c\| \le r \\ Line 49  q^c  & & \|x-x^c\| \le r \\
49  for all $x$ in the domain and all time  $t>0$.  for all $x$ in the domain and all time  $t>0$.
50
51  On the surface of the domain we are  On the surface of the domain we are
53  which precribes the normal component of the flux $J\hackscore i= \kappa T\hackscore{,i}$ to be proportional  which precribes the normal component of the flux $\kappa T\hackscore{,i}$ to be proportional
54  to the difference of the current temperature to the surrounding temperature $T\hackscore{ref}$:      to the difference of the current temperature to the surrounding temperature $T\hackscore{ref}$:
55
56   \kappa T\hackscore{,i} n\hackscore i = \eta (T\hackscore{ref}-T)   \kappa T\hackscore{,i} n\hackscore i = \eta (T\hackscore{ref}-T)
57  \label{DIFFUSION TEMP EQ 2}  \label{DIFFUSION TEMP EQ 2}
58
59  $\eta$ is a given material coefficient depending on the material and the surrounding medium.  $\eta$ is a given material coefficient depending on the material of the block and the surrounding medium.
60  As usual $n_i$ is the $i$-th component of the outer normal field \index{outer normal field}  As usual $n\hackscore i$ is the $i$-th component of the outer normal field \index{outer normal field}
61  at the surface of the domain.  at the surface of the domain.
62
63  To solve the the time depended \eqn{DIFFUSION TEMP EQ 1} the initial temperature at time  To solve the time dependent \eqn{DIFFUSION TEMP EQ 1} the initial temperature at time
64  $t=0$ has to be given. Here we assume that the initial temperature is the surrounding temperature:  $t=0$ has to be given. Here we assume that the initial temperature is the surrounding temperature:
65
66  T(x,0)=T\hackscore{ref}  T(x,0)=T\hackscore{ref}
67  \label{DIFFUSION TEMP EQ 4}  \label{DIFFUSION TEMP EQ 4}
68
69  for all $x$ in the domain. It is pointed out that  for all $x$ in the domain. It is pointed out that
70  the initial conditions is fullfilling the  the initial conditions satisfy the
71  boundary condition defined by \eqn{DIFFUSION TEMP EQ 2}.  boundary condition defined by \eqn{DIFFUSION TEMP EQ 2}.
72
73  The temperature is calculated discrete time nodes $t^{(n)}$ where  The temperature is calculated at discrete time nodes $t^{(n)}$ where
74  $t^{(0)}=0$ and  $t^{(n)}=t^{(n-1)}+h$ where $h>0$ is the step size which is assumed to be constant.  $t^{(0)}=0$ and  $t^{(n)}=t^{(n-1)}+h$ where $h>0$ is the step size which is assumed to be constant.
75  In the following the upper index ${(n)}$ is refering to a value at time $t^{(n)}$. The simplest  In the following the upper index ${(n)}$ refers to a value at time $t^{(n)}$. The simplest
76  and most robust scheme to approximate the time derivative of the the temperature is the  and most robust scheme to approximate the time derivative of the the temperature is
77  \index{backward Euler} scheme, see~\cite{XXX} for alternatives. The backward Euler scheme bases  the backward Euler
78  on the Taylor expansion  \index{backward Euler} scheme, see~\cite{XXX} for alternatives. The backward Euler
79    scheme is based
80    on the Taylor expansion of $T$ at time $t^{(n)}$:
81
82  T^{(n-1)}\approx T^{(n)}+T\hackscore{,t}^{(n)}(t^{(n-1)}-t^{(n)})  T^{(n-1)}\approx T^{(n)}+T\hackscore{,t}^{(n)}(t^{(n-1)}-t^{(n)})
83  =T^{(n-1)} + h \cdot T\hackscore{,t}^{(n)}  =T^{(n-1)} - h \cdot T\hackscore{,t}^{(n)}
84  \label{DIFFUSION TEMP EQ 6}  \label{DIFFUSION TEMP EQ 6}
85
86  which is inserted into \eqn{DIFFUSION TEMP EQ 1}. By separating the terms at  This is inserted into \eqn{DIFFUSION TEMP EQ 1}. By separating the terms at
87  $t^{(n)}$ and  $t^{(n-1)}$ one gets for $n=1,2,3\ldots$  $t^{(n)}$ and  $t^{(n-1)}$ one gets for $n=1,2,3\ldots$
88
89  \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)}  \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)}
90  \label{DIFFUSION TEMP EQ 7}  \label{DIFFUSION TEMP EQ 7}
91
92  where $T^{(0)}=T\hackscore{ref}$ is taken form the initial condition given by \eqn{DIFFUSION TEMP EQ 4}.  where $T^{(0)}=T\hackscore{ref}$ is taken form the initial condition given by \eqn{DIFFUSION TEMP EQ 4}.
93  Together with the natural boundary condition from \eqn{DIFFUSION TEMP EQ 2}.  Together with the natural boundary condition
94
95     \kappa T\hackscore{,i}^{(n)} n\hackscore i = \eta (T\hackscore{ref}-T^{(n)})
96    \label{DIFFUSION TEMP EQ 2222}
97
98    taken from \eqn{DIFFUSION TEMP EQ 2}
99  this forms a boundary value problem that has to be solved for each time step.  this forms a boundary value problem that has to be solved for each time step.
100  As a first step to implement a solver for the temperature diffusion problem we will  As a first step to implement a solver for the temperature diffusion problem we will
101  first implement a solver for the  boundary value problem that has to be solved at each time step.  first implement a solver for the  boundary value problem that has to be solved at each time step.
# Line 113  where $u$ plays the role of $T^{(n)}$ an Line 111  where $u$ plays the role of $T^{(n)}$ an
111  \omega=\frac{\rho c\hackscore p}{h} \mbox{ and } f=q+\frac{\rho c\hackscore p}{h}T^{(n-1)} \;.  \omega=\frac{\rho c\hackscore p}{h} \mbox{ and } f=q+\frac{\rho c\hackscore p}{h}T^{(n-1)} \;.
112  \label{DIFFUSION HELM EQ 1b}  \label{DIFFUSION HELM EQ 1b}
113
114  With $g=\eta T\hackscore{ref}$ the radiation condition defined by \eqn{DIFFUSION TEMP EQ 2}  With $g=\eta T\hackscore{ref}$ the radiation condition defined by \eqn{DIFFUSION TEMP EQ 2222}
115  takes the form  takes the form
116
117  \kappa u\hackscore{,i} n\hackscore{i} =  g - \eta u\mbox{ on } \Gamma  \kappa u\hackscore{,i} n\hackscore{i} =  g - \eta u\mbox{ on } \Gamma
# Line 121  takes the form Line 119  takes the form
119
120  The partial differential  The partial differential
121  \eqn{DIFFUSION HELM EQ 1} together with boundary conditions of \eqn{DIFFUSION HELM EQ 2}  \eqn{DIFFUSION HELM EQ 1} together with boundary conditions of \eqn{DIFFUSION HELM EQ 2}
122  is called a Helmholtz equation \index{Helmholtz equation}.  is called the Helmholtz equation \index{Helmholtz equation}.
123
124  We want to use the \LinearPDE class provided by \escript to define and solve a general linear PDE such as the  We want to use the \LinearPDE class provided by \escript to define and solve a general linear PDE such as the
125  Helmholtz equation. We have used a special case of the \LinearPDE class, namely the  Helmholtz equation. We have used a special case of the \LinearPDE class, namely the
127  Here we will write our own specialized class of the \LinearPDE to solve the Helmholtz equation.  Here we will write our own specialized sub-class of the \LinearPDE to define the Helmholtz equation
128    and use the \method{getSolution} method of parent class to actually solve the problem.
129
130  The general form of a single PDE that can be handeled by the \LinearPDE class is  The form of a single PDE that can be handled by the \LinearPDE class is
131  \label{EQU.FEM.1}  \label{EQU.FEM.1}
132  -(A\hackscore{jl} u\hackscore{,l})\hackscore{,j}+D u =Y \; .  -(A\hackscore{jl} u\hackscore{,l})\hackscore{,j}+D u =Y \; .
133
134    We show here the terms which are relevant for the Helmholtz problem.
135  The general form and systems is discussed in \Sec{SEC LinearPDE}.    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}.  $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.  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  They may be constant or may depend on their
139  location in the domain but must not depend on the unknown solution $u$.  location in the domain but must not depend on the unknown solution $u$.
# Line 145  n\hackscore{j}A\hackscore{jl} u\hackscor Line 145  n\hackscore{j}A\hackscore{jl} u\hackscor
145  where, as usual, $n$ denotes the outer normal field on the surface of the domain. Notice that  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.  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}  By inspecting the Helmholtz equation \index{Helmholtz equation}
149  we can easily assign values to the coefficients in the  we can easily assign values to the coefficients in the
150  general PDE of the \LinearPDE class:  general PDE of the \LinearPDE class:
151  \label{DIFFUSION HELM EQ 3}  \label{DIFFUSION HELM EQ 3}
# Line 154  A\hackscore{ij}=\kappa \delta\hackscore{ Line 154  A\hackscore{ij}=\kappa \delta\hackscore{
154  d=\eta & y= g &  \\  d=\eta & y= g &  \\
155  \end{array}  \end{array}
156
157  $\delta\hackscore{ij}$ is the Kronecker symbol \index{Kronecker symbol} defined by $\delta=\hackscore{ij}=1$ for  $\delta\hackscore{ij}$ is the Kronecker symbol \index{Kronecker symbol} defined by $\delta\hackscore{ij}=1$ for
158  $i=j$ and $0$ otherwise.  $i=j$ and $0$ otherwise.
159
160  We want to implement a  We want to implement a
161  new class which we will call \class{Helmholtz} that provides the same methods like the \LinearPDE class but  new class which we will call \class{Helmholtz} that provides the same methods as the \LinearPDE class but
162  is defined through the coefficeints $\kappa$, $\omega$, $f$, $\eta$,  is described using the coefficients $\kappa$, $\omega$, $f$, $\eta$,
$g$ rather than the general form given by \eqn{EQU.FEM.1}.
Python's
mechanism of inhertence allows doing this in a very easy way.
The advantage is that our new \class{Helmholtz} can be used in any context
that works with a \LinearPDE class but with an easier interface to define the PDE.
This improves reuasablity as well as maintainability of program codes.

We want to implement a
new class which we will call \class{Helmholtz} that provides the same methods like the \LinearPDE class but
is defined through the coefficeints $\kappa$, $\omega$, $f$, $\alpha$,
163  $g$ rather than the general form given by \eqn{EQU.FEM.1}.  $g$ rather than the general form given by \eqn{EQU.FEM.1}.
164  Python's mechanism of subclasses allows doing this in a very easy way.  Python's mechanism of subclasses allows us to do this in a very simple way.
165  The \Poisson class of the \linearPDEsPack module,  The \Poisson class of the \linearPDEsPack module,
166  which we have already used in \Chap{FirstSteps}, is in fact a subclass of the  which we have already used in \Chap{FirstSteps}, is in fact a subclass of the more general
167  \LinearPDE class. That means that it all methods (such as the \method{getSolution})  \LinearPDE class. That means that all methods (such as the \method{getSolution})
168  from the parent class \LinearPDE are defined for any \Poisson object. However, new  from the parent class \LinearPDE are available for any \Poisson object. However, new
169  methods can be added and methods of the parent class can be redefined. In fact,  methods can be added and methods of the parent class can be redefined. In fact,
170  the \Poisson class redefines the \method{setValue} of the \LinearPDE class which is called to assign  the \Poisson class redefines the \method{setValue} of the \LinearPDE class which is called to assign
171  values to the coefficients of the PDE. This is exactly what we will do when we define  values to the coefficients of the PDE. This is exactly what we will do when we define
# Line 183  our new \class{Helmholtz} class: Line 173  our new \class{Helmholtz} class:
173  \begin{python}  \begin{python}
174  from esys.linearPDEs import LinearPDE  from esys.linearPDEs import LinearPDE
175  import numarray  import numarray
176  class Helmholtz(LinearPDE)  class Helmholtz(LinearPDE):
177     def setValue(self,kappa=0,omega=1,f=0,eta=0,g=0)     def setValue(self,kappa=0,omega=1,f=0,eta=0,g=0)
178          self._setValue(A=kappa*numarray.identity(self.getDim()),D=omega,Y=f,d=eta,y=g)          ndim=self.getDim()
179            kronecker=numarray.identity(ndim)
180            self._setValue(A=kappa*kronecker,D=omega,Y=f,d=eta,y=g)
181  \end{python}  \end{python}
182  \code{class Helmholtz(linearPDE)} declares the new \class{Helmholtz} class as a subclass  \code{class Helmholtz(linearPDE)} declares the new \class{Helmholtz} class as a subclass
183  of the \LinearPDE which have imported in the first line of the script.  of \LinearPDE which we have imported in the first line of the script.
184  We add the method \method{setValue} to the class which overwrites the  We add the method \method{setValue} to the class which overwrites the
185  \method{setValue} method of the \LinearPDE class. The new methods which has the  \method{setValue} method of the \LinearPDE class. The new method which has the
186  parameters of the Helmholtz \eqn{DIFFUSION HELM EQ 1} as arguments  parameters of the Helmholtz \eqn{DIFFUSION HELM EQ 1} as arguments
187  maps the parameters of the coefficients of the general PDE defined  maps the parameters of the coefficients of the general PDE defined
188  in \eqn{EQU.FEM.1}. The coefficient \var{A} is defined by Kroneckers symbol. we use the  in \eqn{EQU.FEM.1}. We are actually using the \method{_setValue} of
189  \numarray function \function{identity} which return a square matrix which has ones on the  the \LinearPDE class (notice the leeding underscoure).
190  main diagonal and zeros out side the main diagonal. The argument of \function{identity} gives the order of the matrix.  The coefficient \var{A} involves the Kronecker symbol. We use the
191  Here we use  \numarray function \function{identity} which returns a square matrix with ones on the
192  the \method{getDim} of the \LinearPDE class object \var{self} to get the spatial dimension of the domain of the  main diagonal and zeros off the main diagonal. The argument of \function{identity} gives the order of the matrix.
193  PDE. As we will make use of the \class{Helmholtz} class several times, it is convient to  The \method{getDim} of the \LinearPDE class object \var{self} to get the spatial dimensions of the domain of the
194  put is definition into a file which we name \file{mytools.py} available in the \ExampleDirectory.  PDE. As we will make use of the \class{Helmholtz} class several times, it is convenient to
195    put its definition into a file which we name \file{mytools.py} available in the \ExampleDirectory.
196  You can use your favourite editor to create and edit the file.    You can use your favourite editor to create and edit the file.
197
198  An object of the \class{Helmholtz} class is created through the statments:  An object of the \class{Helmholtz} class is created through the statments:
199  \begin{python}  \begin{python}
200  from mytools import *  from mytools import Helmholtz
201  mypde=Helmholtz(mydomain)  mypde = Helmholtz(mydomain)
202  mypde.setValue(kappa=10.,omega=0.1,f=12)  mypde.setValue(kappa=10.,omega=0.1,f=12)
203  u=mypde.getSolution()  u = mypde.getSolution()
204  \end{python}  \end{python}
205  In the first statement we import all definition from the \file{mytools.py}  \index{scripts!\file{mytools.py}}. Make sure  In the first statement we import all definition from the \file{mytools.py} \index{scripts!\file{mytools.py}}. Make sure
206  that \file{mytools.py} is in directory from where you have started Python.  that \file{mytools.py} is in the directory from where you started Python.
207  \var{mydomain} is the \Domain of the PDE. In the third statment values are  \var{mydomain} is the \Domain of the PDE which we have undefined here. In the third statment values are
208  assigned to the PDE parameters. As no values for arguments \var{eta} and \var{g} are  assigned to the PDE parameters. As no values for arguments \var{eta} and \var{g} are
209  specified the default values $0$ are used \footnote{It would be better to use the default value  specified, the default values $0$ are used. \footnote{It would be better to use the default value
210  \var{escript.Data()} rather then $0$ as then the coefficient would be defined as being not present and  \var{escript.Data()} rather then $0$ as then the coefficient would be defined as being not present and
211  would not be processed when the PDE is evaluated.}. In the forth statement the solution of the  would not be processed when the PDE is evaluated}. In the fourth statement the solution of the
212  PDE is returned.  PDE is returned.
213
214  We want to test our \class{Helmholtz} class on a rectangular domain  To test our \class{Helmholtz} class on a rectangular domain
215  of length $l$ and height $h$. We do this by choosing a simple test solution,  of length $l\hackscore{0}=5$ and height $l\hackscore{1}=1$, we choose a simple test solution. Actually, we
216  here we take $u=x\hackscore{0}$ and then calculate the right hand side terms $f$ and $g$ such that  we take $u=x\hackscore{0}$ and then calculate the right hand side terms $f$ and $g$ such that
217  the test solution becomes the solution of the problem. If we take $\kappa=1$  the test solution becomes the solution of the problem. If we assume $\kappa$ as being constant,
218  an easy calculation shows that we have to choose $f=\omega \cdot x\hackscore{0}$. On the boundary we get  an easy calculation shows that we have to choose $f=\omega \cdot x\hackscore{0}$. On the boundary we get
219  $\kappa n\hackscore{i} u\hackscore{,i}=n\hackscore{0}$.    $\kappa n\hackscore{i} u\hackscore{,i}=\kappa n\hackscore{0}$.
220  So we have to set $g=n\hackscore{0}+\eta x\hackscore{0}$. The following script \file{helmholtztest.py}  So we have to set $g=\kappa n\hackscore{0}+\eta x\hackscore{0}$. The following script \file{helmholtztest.py}
221  \index{scripts!\file{helmholtztest.py}} which is available in the \ExampleDirectory  \index{scripts!\file{helmholtztest.py}} which is available in the \ExampleDirectory
222  implements this test problem using the \finley PDE solver:  implements this test problem using the \finley PDE solver:
223  \begin{python}  \begin{python}
224  from mytools import *  from mytools import Helmholtz
225  from esys.escript import *  from esys.escript import Lsup
226  import esys.finley  from esys.finley import Rectangle
227  #... set some parameters ...  #... set some parameters ...
228    kappa=1.
229  omega=0.1  omega=0.1
230  eta=10.  eta=10.
231  #... generate domain ...  #... generate domain ...
232  mydomain = esys.finley.Rectangle(l0=5.,l1=1.,n0=50, n1=10)  mydomain = Rectangle(l0=5.,l1=1.,n0=50, n1=10)
233  #... open PDE and set coefficients ...  #... open PDE and set coefficients ...
234  mypde=Helmholtz(mydomain)  mypde=Helmholtz(mydomain)
235  n=mydomain.getNormal()  n=mydomain.getNormal()
236  x=mydomain.getX()  x=mydomain.getX()
237  mypde.setValue(1,omega,omega*x[0],eta,n[0]+eta*x[0])  mypde.setValue(kappa,omega,omega*x[0],eta,kappa*n[0]+eta*x[0])
238  #... calculate error of the PDE solution ...  #... calculate error of the PDE solution ...
239  u=mypde.getSolution()  u=mypde.getSolution()
240  print "error is ",Lsup(u-x[0])  print "error is ",Lsup(u-x[0])
241  \end{python}  \end{python}
242  The script is similar to the script \file{mypoisson.py} dicussed in \Chap{FirstSteps}.  The script is similar to the script \file{mypoisson.py} dicussed in \Chap{FirstSteps}.
243  \code{mydomain.getNormal()} returns the outer normal field on the surface of the domain. The function \function{Lsup}  \code{mydomain.getNormal()} returns the outer normal field on the surface of the domain. The function \function{Lsup}
244  is imported by the \code{from escript import *} statement and returns the maximum absulute value of it argument. To run  is imported by the \code{from esys.escript import Lsup} statement and returns the maximum absulute value of its argument.
245  The error shown by the print statement should be in the order of $10^{-7}$. As piecewise linear interpolation is  The error shown by the print statement should be in the order of $10^{-7}$. As piecewise bi-linear interpolation is
246  used to approximate the solution and our solution is a linear function of the spatial coordinates one may  used to approximate the solution and our solution is a linear function of the spatial coordinates one might
247  expect that the error is zero. However, as most PDE packages uses an iterative solver which is terminated  expect that the error would be zero or in the order of machine precision (typically $\approx 10^{-15}$).
248  when a given tolerance has been reached. The default tolerance is $10^{-8}$. Thsi value can be altered by using the  However most PDE packages use an iterative solver which is terminated
249    when a given tolerance has been reached. The default tolerance is $10^{-8}$. This value can be altered by using the
250  \method{setTolerance} of the \LinearPDE class.  \method{setTolerance} of the \LinearPDE class.
251
252  \section{The Transition Problem}  \section{The Transition Problem}
253  \label{DIFFUSION TRANS SEC}  \label{DIFFUSION TRANS SEC}
254  Now we are ready to solve the original time dependent problem. The main  Now we are ready to solve the original time dependent problem. The main
255  part of the script is the loop over the time $t$ which takes the following form:  part of the script is the loop over time $t$ which takes the following form:
256  \begin{python}  \begin{python}
257    t=0
258    T=Tref
259  mypde=Helmholtz(mydomain)  mypde=Helmholtz(mydomain)
260  while t<t_end:  while t<t_end:
261        mypde.setValue(kappa,rhocp/h,q+rhocp/h*T,eta,eta*Tref)        mypde.setValue(kappa,rhocp/h,q+rhocp/h*T,eta,eta*Tref)
# Line 266  while t<t_end: Line 263  while t<t_end:
263        t+=h        t+=h
264  \end{python}  \end{python}
265  \var{kappa}, \var{rhocp}, \var{eta} and \var{Tref} are input parameters of the model. \var{q} is the heat source  \var{kappa}, \var{rhocp}, \var{eta} and \var{Tref} are input parameters of the model. \var{q} is the heat source
266  in the domain and \var{h} is the time step size which has to be chosen. Notice that the \class{Hemholtz}  in the domain and \var{h} is the time step size. Notice that the \class{Hemholtz} class object \var{mypde}
267  is created before the loop over time is entered while in each time step only the coefficients  is created before the loop over time is entered while in each time step only the coefficients
268  are reset in each time step. This way some information about the reperesentation of the PDE can be reused  are reset in each time step. This way some information about the representation of the PDE can be reused
269  \footnote{The efficience can be improved further by setting the coefficients in the operator  \footnote{The efficience can be improved further by setting the coefficients in the operator
270  \var{kappa}, \var{omega} and \var{eta} before entering the \code{while}-loop and only update the coefficients  \var{kappa}, \var{omega} and \var{eta} before entering the \code{while}-loop and only updating the coefficients
271  in the right hand side \var{f} and \var{g}. This needs a more careful implementation of the \method{setValue}  in the right hand side \var{f} and \var{g}. This needs a more careful implementation of the \method{setValue}
272  method but gives the advantage that the \LinearPDE class can save rebuilding the PDE operator}. The variable \var{T}  method but gives the advantage that the \LinearPDE class can save rebuilding the PDE operator.}. The variable \var{T}
273  holds the current temperature. It is used to calculate the right hand side coefficient \var{f} in the  holds the current temperature. It is used to calculate the right hand side coefficient \var{f} in the
274  Helmholtz \eqn{DIFFUSION HELM EQ 1}. Statement \code{T=mypde.getSolution()} overwrites \var{T} with the  Helmholtz equation in \eqn{DIFFUSION HELM EQ 1}. The statement \code{T=mypde.getSolution()} overwrites \var{T} with the
275  temperature of the new time step $\var{t}+\var{h}$. To get this iterative process going we need to sepcify the  temperature of the new time step $\var{t}+\var{h}$. To get this iterative process going we need to specify the
276  initial temperature distribution, which equal to $T\hackscore{ref}$.  initial temperature distribution, which equal to $T\hackscore{ref}$.
277
278  The heat source \var{q} which is defined in \eqn{DIFFUSION TEMP EQ 1b} shall be \var{q0}  The heat source \var{q} which is defined in \eqn{DIFFUSION TEMP EQ 1b} is \var{qc}
279  at an area defined as a circle of radius \var{r} and center \var{xc} and zero outside this circle.  in an area defined as a circle of radius \var{r} and center \var{xc} and zero outside this circle.
280  \var{q0} is a fixed constant. The following script defines \var{q} as desired:    \var{q0} is a fixed constant. The following script defines \var{q} as desired:
281  \begin{python}  \begin{python}
282    from esys.escript import length
283  xc=[0.02,0.002]  xc=[0.02,0.002]
284  r=0.001  r=0.001
285  x=mydomain.getX()  x=mydomain.getX()
286  q=q0*(length(x-xc)-r).whereNegative()  q=q0*(length(x-xc)-r).whereNegative()
287  \end{python}  \end{python}
288  \var{x} is a \Data class object of  \var{x} is a \Data class object of
289  the \escript module defining the locations of points in the \Domain \var{mydomain}.  the \escript module defining locations in the \Domain \var{mydomain}.
290  \code{length(x-xc)} calculates the distances in the Euclidean norm  The \function{length()} imported from the \escript module returns the
291  of the locations \var{x} to the center of the circle \var{xc} where the heat source is acting.  Euclidean norm:
292  Notice that the coordinates of \var{xc} are defined as a list of floating point numbers. It is independently  \label{DIFFUSION HELM EQ 3aba}
293  converted into a \Data class object before subtracted from \var{x}. The method \method{whereNegative} of  \|y\|=\sqrt{
294    y\hackscore{i}
295    y\hackscore{i}
296    } = \function{esys.escript.length}(y)
297
298    So \code{length(x-xc)} calculates the distances
299    of the location \var{x} to the center of the circle \var{xc} where the heat source is acting.
300    Note that the coordinates of \var{xc} are defined as a list of floating point numbers. It is independently
301    converted into a \Data class object before being subtracted from \var{x}. The method \method{whereNegative} of
302  a \Data class object, in this case the result of the expression  a \Data class object, in this case the result of the expression
303  \code{length(x-xc)-r}, returns a \Data class which is equal one where the object negative and  \code{length(x-xc)-r}, returns a \Data class which is equal to one where the object is negative and
304  zero elsewhere. After multiplication with \var{q0} we get a function with the deired property.  zero elsewhere. After multiplication with \var{qc} we get a function with the desired property.
305
306  Now we can put the components together to the script \file{diffusion.py} which is available in the \ExampleDirectory:  Now we can put the components together to create the script \file{diffusion.py} which is available in the \ExampleDirectory:
307  \index{scripts!\file{diffusion.py}}:  \index{scripts!\file{diffusion.py}}:
308  \begin{python}  \begin{python}
309  from mytools import *  from mytools import Helmholtz
310  from esys.escript import *  from esys.escript import Lsup
311  import esys.finley  from esys.finley import Rectangle
312  #... set some parameters ...  #... set some parameters ...
313  x_c=[0.02,0.002]  xc=[0.02,0.002]
314  r=0.001  r=0.001
315  q0=50.e6  qc=50.e6
316  Tref=0.  Tref=0.
317  rhocp=2.6e6  rhocp=2.6e6
318  eta=75.  eta=75.
319  kappa=240.  kappa=240.
320  t_end=5.  tend=5.
321  # ...time step size and counter ...  # ... time, time step size and counter ...
322    t=0
323  h=0.1  h=0.1
324  i=0  i=0
t=0
325  #... generate domain ...  #... generate domain ...
326  mydomain = esys.finley.Rectangle(l0=0.05,l1=0.01,n0=250, n1=50)  mydomain = Rectangle(l0=0.05,l1=0.01,n0=250, n1=50)
327  #... open PDE ...  #... open PDE ...
328  mypde=Helmholtz(mydomain)  mypde=Helmholtz(mydomain)
329  # ... set heat source: ....  # ... set heat source: ....
330  x=mydomain.getX()  x=mydomain.getX()
331  q=q0*(length(x-x_c)-r).whereNegative()  q=qc*(length(x-xc)-r).whereNegative()
332  # ... set initial temperature ....  # ... set initial temperature ....
333  T=Tref  T=Tref
334  # ... start iteration:  # ... start iteration:
335  while t<t_end:  while t<tend:
336        i+=1        i+=1
337        t+=h        t+=h
338        print "time step :",t        print "time step :",t
# Line 337  while t<t_end: Line 343  while t<t_end:
343  The script will create the files \file{T.1.dx},  The script will create the files \file{T.1.dx},
344   \file{T.2.dx}, $\ldots$, \file{T.50.dx} in the directory where the script has been started. The files give the   \file{T.2.dx}, $\ldots$, \file{T.50.dx} in the directory where the script has been started. The files give the
345  temperature distributions at time steps $1$, $2$, $\ldots$, $50$ in the \OpenDX file format.  temperature distributions at time steps $1$, $2$, $\ldots$, $50$ in the \OpenDX file format.
346
347    \begin{figure}
348    \centerline{\includegraphics[width=\figwidth]{DiffusionRes1}}
349    \centerline{\includegraphics[width=\figwidth]{DiffusionRes16}}
350    \centerline{\includegraphics[width=\figwidth]{DiffusionRes32}}
351    \centerline{\includegraphics[width=\figwidth]{DiffusionRes48}}
352    \caption{Results of the Temperture Diffusion Problem for Time Steps $1$ $16$, $32$ and $48$.}
353    \label{DIFFUSION FIG 2}
354    \end{figure}
355
356  An easy way to visualize the results is the command  An easy way to visualize the results is the command
357  \begin{verbatim}  \begin{verbatim}
358  dx -edit diffusion.net  dx -edit diffusion.net &
359  \end{verbatim}  \end{verbatim}
360  where \file{diffusion.net} is an \OpenDX script available in the \ExampleDirectory.  where \file{diffusion.net} is an \OpenDX script available in the \ExampleDirectory.
361  \fig{DIFFUSION FIG 2} shows the result for some selected time steps.  Use the \texttt{Sequencer} to move forward and and backwards in time.
362    \fig{DIFFUSION FIG 2} shows the result for some selected time steps.

Legend:
 Removed from v.106 changed lines Added in v.107