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

revision 568 by gross, Tue Feb 28 03:59:06 2006 UTC revision 569 by gross, Tue Feb 28 05:34:37 2006 UTC
# Line 173  where we assume that \code{mydomain} is Line 173  where we assume that \code{mydomain} is
173  \code{kappa} \code{omega} \code{eta} and \code{g} are given scalar values  \code{kappa} \code{omega} \code{eta} and \code{g} are given scalar values
174  typically \code{float} or \Data objects. The \method{setValue} method  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  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  the PDE and returns the solution \code{u} of the PDE. \function{kronecker} is \escript function
177  returning the Kronecker symbol.  returning the Kronecker symbol.
178
179  The coefficients can set by several calls of \method{setValue} where the order can be chosen arbitrarily.  The coefficients can set by several calls of \method{setValue} where the order can be chosen arbitrarily.
# Line 193  A\hackscore{jl}=A\hackscore{lj} \mbox{ a Line 193  A\hackscore{jl}=A\hackscore{lj} \mbox{ a
193
194  Note that $D$ and $d$ may have any value and the right hand side $X$, $Y$, $y$ as well as the constraints  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.  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.  The \LinearPDE class provides the method \method{checkSymmetry} method to check if the given PDE is symmetric.
197  \begin{python}  \begin{python}
198  mypde=LinearPDE(mydomain)  mypde=LinearPDE(mydomain)
199  mypde.setValue(A=kappa*kronecker(mydomain),d=eta)  mypde.setValue(A=kappa*kronecker(mydomain),d=eta)
# Line 203  print mypde.checkSymmetry() # returns Fa Line 203  print mypde.checkSymmetry() # returns Fa
203  mypde.setValue(C=kronecker(mydomain)[0])  mypde.setValue(C=kronecker(mydomain)[0])
204  print mypde.checkSymmetry() # returns True  print mypde.checkSymmetry() # returns True
205  \end{python}  \end{python}
206  Unfortunately, a \code{checkSymmetry} is very expensive and is recommendable to use for  Unfortunately, a \method{checkSymmetry} is very expensive and is recommendable to use for
207  testing and debugging purposes only. The \code{setSymmetryOn} method is used to  testing and debugging purposes only. The \method{setSymmetryOn} method is used to
208  declare a PDE symmetric:  declare a PDE symmetric:
209  \begin{python}  \begin{python}
210  mypde = LinearPDE(mydomain)  mypde = LinearPDE(mydomain)
211  mypde.setValue(A=kappa*kronecker(mydomain))  mypde.setValue(A=kappa*kronecker(mydomain))
212  mypde.setSymmetryOn()  mypde.setSymmetryOn()
213  \end{python}  \end{python}
214    Now we want to see how we actually solve the Helmholtz equation.
215    on a rectangular domain
216    of length $l\hackscore{0}=5$ and height $l\hackscore{1}=1$. We choose a simple test solution such that we
217    can verify the returned solution against the exact answer. Actually, we

=====

Here we will write our own specialized sub-class of the \LinearPDE to define the Helmholtz equation
and use the \method{getSolution} method of parent class to actually solve the problem.

We want to implement a
new class which we will call \class{Helmholtz} that provides the same methods as the \LinearPDE class but
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 subclasses allows us to do this in a very simple way.
That means that all methods (such as the \method{getSolution})
from the parent class \LinearPDE are available for any \Poisson object. However, new
methods can be added and methods of the parent class can be redefined. In fact,
the \Poisson class redefines the \method{setValue} of the \LinearPDE class which is called to assign
values to the coefficients of the PDE. This is exactly what we will do when we define
our new \class{Helmholtz} class:

To test our \class{Helmholtz} class on a rectangular domain
of length $l\hackscore{0}=5$ and height $l\hackscore{1}=1$, we choose a simple test solution. Actually, we
218  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
219  the test solution becomes the solution of the problem. If we assume $\kappa$ as being constant,  the test solution becomes the solution of the problem. If we assume $\kappa$ as being constant,
220  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
221  $\kappa n\hackscore{i} u\hackscore{,i}=\kappa n\hackscore{0}$.    $\kappa n\hackscore{i} u\hackscore{,i}=\kappa n\hackscore{0}$.
222  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{helmholtz.py}
223  \index{scripts!\file{helmholtz.py}} which is available in the \ExampleDirectory  \index{scripts!\file{helmholtz.py}} which is available in the \ExampleDirectory
224  implements this test problem using the \finley PDE solver:  implements this test problem using the \finley PDE solver:
225  \begin{python}  \begin{python}
226  from esys.escript import *  from esys.escript import *
227  from linearPDEs import LinearPDE  from esys.escript.linearPDEs import LinearPDE
228  from esys.finley import Rectangle  from esys.finley import Rectangle
229  #... set some parameters ...  #... set some parameters ...
230  kappa=1.  kappa=1.
# Line 259  mypde=LinearPDE(mydomain) Line 237  mypde=LinearPDE(mydomain)
237  mypde.setSymmetryOn()  mypde.setSymmetryOn()
238  n=mydomain.getNormal()  n=mydomain.getNormal()
239  x=mydomain.getX()  x=mydomain.getX()
240  mypde.setValue(A=kappa*kronecker(mydomain),D=omega,Y=omega*x[0],d=eta,y=kappa*n[0]+eta*x[0])  mypde.setValue(A=kappa*kronecker(mydomain),D=omega,Y=omega*x[0], \
241                   d=eta,y=kappa*n[0]+eta*x[0])
242  #... calculate error of the PDE solution ...  #... calculate error of the PDE solution ...
243  u=mypde.getSolution()  u=mypde.getSolution()
244  print "error is ",Lsup(u-x[0])  print "error is ",Lsup(u-x[0])
245  \end{python}  \end{python}
246  The script is similar to the script \file{mypoisson.py} dicussed in \Chap{FirstSteps}.  The script is similar to the script \file{poisson.py} dicussed in \Chap{FirstSteps}.
247  \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}
248  is imported by the \code{from esys.escript import Lsup} statement and returns the maximum absulute value of its argument.  imported by the \code{from esys.escript import *} statement and returns the maximum absulute value of its argument.
249  The error shown by the print statement should be in the order of $10^{-7}$. As piecewise bi-linear interpolation is  The error shown by the print statement should be in the order of $10^{-7}$. As piecewise bi-linear interpolation is
250  used to approximate the solution and our solution is a linear function of the spatial coordinates one might  used by \finley approximate the solution and our solution is a linear function of the spatial coordinates one might
251  expect that the error would be zero or in the order of machine precision (typically $\approx 10^{-15}$).  expect that the error would be zero or in the order of machine precision (typically $\approx 10^{-15}$).
252  However most PDE packages use an iterative solver which is terminated  However most PDE packages use an iterative solver which is terminated
253  when a given tolerance has been reached. The default tolerance is $10^{-8}$. This value can be altered by using the  when a given tolerance has been reached. The default tolerance is $10^{-8}$. This value can be altered by using the
# Line 281  part of the script is the loop over time Line 260  part of the script is the loop over time
260  \begin{python}  \begin{python}
261  t=0  t=0
262  T=Tref  T=Tref
263  mypde=Helmholtz(mydomain)  mypde=LinearPDE(mydomain)
264    mypde.setValue(A=kappa*kronecker(mydomain),D=rhocp/h,d=eta,y=eta*Tref)
265  while t<t_end:  while t<t_end:
266        mypde.setValue(kappa,rhocp/h,q+rhocp/h*T,eta,eta*Tref)        mypde.setValue(Y=q+rhocp/h*T)
267        T=mypde.getSolution()        T=mypde.getSolution()
268        t+=h        t+=h
269  \end{python}  \end{python}
270  \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
271  in the domain and \var{h} is the time step size. Notice that the \class{Hemholtz} class object \var{mypde}  in the domain and \var{h} is the time step size.
272  is created before the loop over time is entered while in each time step only the coefficients  The variable \var{T}
are reset in each time step. This way some information about the representation of the PDE can be reused
\footnote{The efficience can be improved further by setting the coefficients in the operator
\var{kappa}, \var{omega} and \var{eta} before entering the \code{while}-loop and only updating the coefficients
in the right hand side \var{f} and \var{g}. This needs a more careful implementation of the \method{setValue}
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 equation in \eqn{DIFFUSION HELM EQ 1}. The 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 specify 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    The \LinearPDE class object \var{mypde}
278    and coefficients not changing over time are set up before the loop over time is entered. In each time step only the coefficient
279    $Y$ is reset as it depends on the temperature of the previous time step. This allows the PDE solver to reuse information
280    from previous time steps as much as possible.
281
282  The heat source \var{q} which is defined in \eqn{DIFFUSION TEMP EQ 1b} is \var{qc}  The heat source \var{q} which is defined in \eqn{DIFFUSION TEMP EQ 1b} is \var{qc}
283  in 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.
# Line 308  from esys.escript import length Line 287  from esys.escript import length
287  xc=[0.02,0.002]  xc=[0.02,0.002]
288  r=0.001  r=0.001
289  x=mydomain.getX()  x=mydomain.getX()
290  q=q0*(length(x-xc)-r).whereNegative()  q=q0*whereNegative(length(x-xc)-r)
291  \end{python}  \end{python}
292  \var{x} is a \Data class object of  \var{x} is a \Data class object of
293  the \escript module defining locations in the \Domain \var{mydomain}.  the \escript module defining locations in the \Domain \var{mydomain}.
# Line 323  y\hackscore{i} Line 302  y\hackscore{i}
302  So \code{length(x-xc)} calculates the distances    So \code{length(x-xc)} calculates the distances
303  of the location \var{x} to the center of the circle \var{xc} where the heat source is acting.  of the location \var{x} to the center of the circle \var{xc} where the heat source is acting.
304  Note that the coordinates of \var{xc} are defined as a list of floating point numbers. It is independently  Note that the coordinates of \var{xc} are defined as a list of floating point numbers. It is independently
305  converted into a \Data class object before being subtracted from \var{x}. The method \method{whereNegative} of  converted into a \Data class object before being subtracted from \var{x}. The function \function{whereNegative}
306  a \Data class object, in this case the result of the expression  applied to
307  \code{length(x-xc)-r}, returns a \Data class which is equal to one where the object is negative and  \code{length(x-xc)-r}, returns a \Data class which is equal to one where the object is negative and
308  zero elsewhere. After multiplication with \var{qc} we get a function with the desired property.  zero elsewhere. After multiplication with \var{qc} we get a function with the desired property.
309
310  Now we can put the components together to create 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:
311  \index{scripts!\file{diffusion.py}}:  \index{scripts!\file{diffusion.py}}:
312  \begin{python}  \begin{python}
313  from mytools import Helmholtz  from esys.escript import *
314  from esys.escript import Lsup  from esys.escript.linearPDEs import LinearPDE
315  from esys.finley import Rectangle  from esys.finley import Rectangle
316  #... set some parameters ...  #... set some parameters ...
317  xc=[0.02,0.002]  xc=[0.02,0.002]
# Line 350  i=0 Line 329  i=0
329  #... generate domain ...  #... generate domain ...
330  mydomain = Rectangle(l0=0.05,l1=0.01,n0=250, n1=50)  mydomain = Rectangle(l0=0.05,l1=0.01,n0=250, n1=50)
331  #... open PDE ...  #... open PDE ...
332  mypde=Helmholtz(mydomain)  mypde=LinearPDE(mydomain)
333    mypde.setSymmetryOn()
334    mypde.setValue(A=kappa*kronecker(mydomain),D=rhocp/h,d=eta,y=eta*Tref)
335  # ... set heat source: ....  # ... set heat source: ....
336  x=mydomain.getX()  x=mydomain.getX()
337  q=qc*(length(x-xc)-r).whereNegative()  q=qc*whereNegative(length(x-xc)-r)
338  # ... set initial temperature ....  # ... set initial temperature ....
339  T=Tref  T=Tref
340  # ... start iteration:  # ... start iteration:
# Line 361  while t<tend: Line 342  while t<tend:
342        i+=1        i+=1
343        t+=h        t+=h
344        print "time step :",t        print "time step :",t
345        mypde.setValue(kappa=kappa,omega=rhocp/h,f=q+rhocp/h*T,eta=eta,g=eta*Tref)        mypde.setValue(Y=q+rhocp/h*T)
346        T=mypde.getSolution()        T=mypde.getSolution()
347        T.saveDX("T%d.dx"%i)        saveVTK("T.%d.xml"%i,temp=T)
348  \end{python}  \end{python}
349  The script will create the files \file{T.1.dx},  The script will create the files \file{T.1.xml},
350   \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.xml}, $\ldots$, \file{T.50.xml} in the directory where the script has been started. The files give the
351  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 \VTK file format.
352
353  \begin{figure}  \begin{figure}
354  \centerline{\includegraphics[width=\figwidth]{DiffusionRes1}}  \centerline{\includegraphics[width=\figwidth]{DiffusionRes1}}
# Line 377  temperature distributions at time steps Line 358  temperature distributions at time steps
358  \caption{Results of the Temperture Diffusion Problem for Time Steps $1$ $16$, $32$ and $48$.}  \caption{Results of the Temperture Diffusion Problem for Time Steps $1$ $16$, $32$ and $48$.}
359  \label{DIFFUSION FIG 2}  \label{DIFFUSION FIG 2}
360  \end{figure}  \end{figure}

361  An easy way to visualize the results is the command  An easy way to visualize the results is the command
362  \begin{verbatim}  \begin{verbatim}
363  dx -edit diffusion.net &  mayavi -d T.1.xml -m SurfaceMap &
364  \end{verbatim}  \end{verbatim}
365  where \file{diffusion.net} is an \OpenDX script available in the \ExampleDirectory.  Use the \texttt{Configure Data}
366  Use the \texttt{Sequencer} to move forward and and backwards in time.  to move forward and and backwards in time.
367  \fig{DIFFUSION FIG 2} shows the result for some selected time steps.  \fig{DIFFUSION FIG 2} shows the result for some selected time steps.

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