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

revision 567 by gross, Mon Feb 27 00:04:17 2006 UTC revision 568 by gross, Tue Feb 28 03:59:06 2006 UTC
# Line 23  In Section~\ref{DIFFUSION TRANS SEC} the Line 23  In Section~\ref{DIFFUSION TRANS SEC} the
23  solver for the temperature diffusion problem.  solver for the temperature diffusion problem.
24
25  \subsection{\label{DIFFUSION TEMP SEC}Temperature Diffusion}  \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  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  in the interior of the domain is given by
28  \begin{equation}  \begin{equation}
# Line 117  takes the form Line 116  takes the form
116  \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
117  \label{DIFFUSION HELM EQ 2}  \label{DIFFUSION HELM EQ 2}
118  \end{equation}  \end{equation}
119  The partial differential  The partial differential \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}
120  is called the Helmholtz equation \index{Helmholtz equation}.  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 PDE such as the  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. We have used a special case of the \LinearPDE class, namely the  Helmholtz equation. For a single PDE the \LinearPDE class supports the following form:
124  \Poisson class already in \Chap{FirstSteps}.  \begin{equation}\label{LINEARPDE.SINGLE.1}
125  Here we will write our own specialized sub-class of the \LinearPDE to define the Helmholtz equation  -(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  and use the \method{getSolution} method of parent class to actually solve the problem.  \end{equation}
127    The coefficients $A$, $B$, $C$, $D$, $X$ and $Y$ have to be specified through \Data objects in the
128  The form of a single PDE that can be handled by the \LinearPDE class is  \Function on the PDE or objects that can be converted into such \Data objects.
129  \begin{equation}\label{EQU.FEM.1}  $A$ is a \RankTwo, $B$, $C$ and $X$ are \RankOne and $D$ and $Y$ are scalar.
130  -(A\hackscore{jl} u\hackscore{,l})\hackscore{,j}+D u =Y \; .  The following natural
131  \end{equation}  boundary conditions are considered \index{boundary condition!natural} on $\Gamma$:
132  We show here the terms which are relevant for the Helmholtz problem.  \begin{equation}\label{LINEARPDE.SINGLE.2}
133  The general form and systems is discussed in \Sec{SEC LinearPDE}.    n\hackscore{j}(A\hackscore{jl} u\hackscore{,l}+B\hackscore{j} u)+d u=n\hackscore{j}X\hackscore{j} + y  \;.
134  $A$, $D$ and $Y$ are the known coeffecients of the PDE. \index{partial differential equation!coefficients}  \end{equation}
135  Notice that $A$ is a matrix or tensor of order 2 and $D$ and $Y$ are scalar.  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  They may be constant or may depend on their  each a \Scalar in the \FunctionOnBoundary.  Constraints \index{constraint} for the solution prescribing the value of the
137  location in the domain but must not depend on the unknown solution $u$.  solution at certain locations in the domain. They have the form
138  The following natural boundary conditions \index{boundary condition!natural} that  \begin{equation}\label{LINEARPDE.SINGLE.3}
139  are used in the \LinearPDE class have the form  u=r \mbox{ where } q>0
140  \begin{equation}\label{EQU.FEM.2}  \end{equation}
141  n\hackscore{j}A\hackscore{jl} u\hackscore{,l}+du=y  \;.  $r$ and $q$ are each \Scalar where $q$ is the characteristic function
142  \end{equation}  \index{characteristic function} defining where the constraint is applied.
143  where, as usual, $n$ denotes the outer normal field on the surface of the domain. Notice that  The constraints defined by \eqn{LINEARPDE.SINGLE.3} override any other condition set by \eqn{LINEARPDE.SINGLE.1}
144  the coefficient $A$ is already used in the PDE in \eqn{EQU.FEM.1}. $d$ and $y$ are given scalar coefficients.  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}  By inspecting the Helmholtz equation \index{Helmholtz equation}
151  we can easily assign values to the coefficients in the  we can easily assign values to the coefficients in the
# Line 155  d=\eta & y= g &  \\ Line 157  d=\eta & y= g &  \\
157  \end{array}  \end{array}
158  \end{equation}  \end{equation}
159  $\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
160  $i=j$ and $0$ otherwise.  $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. \code{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 \code{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))
202    print mypde.checkSymmetry() # returns False
203    mypde.setValue(C=kronecker(mydomain))
204    print mypde.checkSymmetry() # returns True
205    \end{python}
206    Unfortunately, a \code{checkSymmetry} is very expensive and is recommendable to use for
207    testing and debugging purposes only. The \code{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
215
216
217
218
219
220
221    =====
222
223    Here we will write our own specialized sub-class of the \LinearPDE to define the Helmholtz equation
224    and use the \method{getSolution} method of parent class to actually solve the problem.
225
226  We want to implement a  We want to implement a
227  new class which we will call \class{Helmholtz} that provides the same methods as the \LinearPDE class but  new class which we will call \class{Helmholtz} that provides the same methods as the \LinearPDE class but
228  is described using the coefficients $\kappa$, $\omega$, $f$, $\eta$,  is described using the coefficients $\kappa$, $\omega$, $f$, $\eta$,
229  $g$ rather than the general form given by \eqn{EQU.FEM.1}.  $g$ rather than the general form given by \eqn{EQU.FEM.1}.
230  Python's mechanism of subclasses allows us to do this in a very simple way.  Python's mechanism of subclasses allows us to do this in a very simple way.
231  The \Poisson class of the \linearPDEs module,  That means that all methods (such as the \method{getSolution})
which we have already used in \Chap{FirstSteps}, is in fact a subclass of the more general
\LinearPDE class. That means that all methods (such as the \method{getSolution})
232  from the parent class \LinearPDE are available for any \Poisson object. However, new  from the parent class \LinearPDE are available for any \Poisson object. However, new
233  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,
234  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
235  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
236  our new \class{Helmholtz} class:  our new \class{Helmholtz} class:
\begin{python}
from esys.linearPDEs import LinearPDE
import numarray
class Helmholtz(LinearPDE):
def setValue(self,kappa=0,omega=1,f=0,eta=0,g=0)
ndim=self.getDim()
kronecker=numarray.identity(ndim)
self._LinearPDE_setValue(A=kappa*kronecker,D=omega,Y=f,d=eta,y=g)
\end{python}
\code{class Helmholtz(linearPDE)} declares the new \class{Helmholtz} class as a subclass
of \LinearPDE which we have imported in the first line of the script.
We add the method \method{setValue} to the class which overwrites the
\method{setValue} method of the \LinearPDE class. The new method which has the
parameters of the Helmholtz \eqn{DIFFUSION HELM EQ 1} as arguments
maps the parameters of the coefficients of the general PDE defined
in \eqn{EQU.FEM.1}. We are actually using the \method{_LinearPDE__setValue} of
the \LinearPDE class.
The coefficient \var{A} involves the Kronecker symbol. We use the
\numarray function \function{identity} which returns a square matrix with ones on the
main diagonal and zeros off the main diagonal. The argument of \function{identity} gives the order of the matrix.
The \method{getDim} of the \LinearPDE class object \var{self} to get the spatial dimensions of the domain of the
PDE. As we will make use of the \class{Helmholtz} class several times, it is convenient to
put its definition into a file which we name \file{mytools.py} available in the \ExampleDirectory.
You can use your favourite editor to create and edit the file.

An object of the \class{Helmholtz} class is created through the statments:
\begin{python}
from mytools import Helmholtz
mypde = Helmholtz(mydomain)
mypde.setValue(kappa=10.,omega=0.1,f=12)
u = mypde.getSolution()
\end{python}
In the first statement we import all definition from the \file{mytools.py} \index{scripts!\file{mytools.py}}. Make sure
that \file{mytools.py} is in the directory from where you started Python.
\var{mydomain} is the \Domain of the PDE which we have undefined here. In the third statment values are
assigned to the PDE parameters. As no values for arguments \var{eta} and \var{g} are
specified, the default values $0$ are used. \footnote{It would be better to use the default value
\var{escript.Data()} rather then $0$ as then the coefficient would be defined as being not present and
would not be processed when the PDE is evaluated}. In the fourth statement the solution of the
PDE is returned.
237
238  To test our \class{Helmholtz} class on a rectangular domain  To test our \class{Helmholtz} class on a rectangular domain
239  of length $l\hackscore{0}=5$ and height $l\hackscore{1}=1$, we choose a simple test solution. Actually, we  of length $l\hackscore{0}=5$ and height $l\hackscore{1}=1$, we choose a simple test solution. Actually, we
# Line 218  the test solution becomes the solution o Line 242  the test solution becomes the solution o
242  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
243  $\kappa n\hackscore{i} u\hackscore{,i}=\kappa n\hackscore{0}$.    $\kappa n\hackscore{i} u\hackscore{,i}=\kappa n\hackscore{0}$.
244  So we have to set $g=\kappa 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}
245  \index{scripts!\file{helmholtztest.py}} which is available in the \ExampleDirectory  \index{scripts!\file{helmholtz.py}} which is available in the \ExampleDirectory
246  implements this test problem using the \finley PDE solver:  implements this test problem using the \finley PDE solver:
247  \begin{python}  \begin{python}
248  from mytools import Helmholtz  from esys.escript import *
249  from esys.escript import Lsup  from linearPDEs import LinearPDE
250  from esys.finley import Rectangle  from esys.finley import Rectangle
251  #... set some parameters ...  #... set some parameters ...
252  kappa=1.  kappa=1.
# Line 231  eta=10. Line 255  eta=10.
255  #... generate domain ...  #... generate domain ...
256  mydomain = Rectangle(l0=5.,l1=1.,n0=50, n1=10)  mydomain = Rectangle(l0=5.,l1=1.,n0=50, n1=10)
257  #... open PDE and set coefficients ...  #... open PDE and set coefficients ...
258  mypde=Helmholtz(mydomain)  mypde=LinearPDE(mydomain)
259    mypde.setSymmetryOn()
260  n=mydomain.getNormal()  n=mydomain.getNormal()
261  x=mydomain.getX()  x=mydomain.getX()
262  mypde.setValue(kappa,omega,omega*x,eta,kappa*n+eta*x)  mypde.setValue(A=kappa*kronecker(mydomain),D=omega,Y=omega*x,d=eta,y=kappa*n+eta*x)
263  #... calculate error of the PDE solution ...  #... calculate error of the PDE solution ...
264  u=mypde.getSolution()  u=mypde.getSolution()
265  print "error is ",Lsup(u-x)  print "error is ",Lsup(u-x)

Legend:
 Removed from v.567 changed lines Added in v.568