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

revision 2370 by gross, Mon Apr 6 06:41:49 2009 UTC revision 2575 by jfenwick, Mon Aug 3 23:38:33 2009 UTC
# Line 1  Line 1
1
2  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
3  %  %
4  % Copyright (c) 2003-2008 by University of Queensland  % Copyright (c) 2003-2009 by University of Queensland
5  % Earth Systems Science Computational Center (ESSCC)  % Earth Systems Science Computational Center (ESSCC)
6  % http://www.uq.edu.au/esscc  % http://www.uq.edu.au/esscc
7  %  %
# Line 33  is the unit square Line 33  is the unit square
33  \label{eq:FirstSteps.1b}  \label{eq:FirstSteps.1b}
34  \end{equation}  \end{equation}
35  The domain is shown in \fig{fig:FirstSteps.1}.  The domain is shown in \fig{fig:FirstSteps.1}.
36  \begin{figure} [h!]  \begin{figure} [ht]
37  \centerline{\includegraphics[width=\figwidth]{figures/FirstStepDomain}}  \centerline{\includegraphics[width=\figwidth]{figures/FirstStepDomain}}
38  \caption{Domain $\Omega=[0,1]^2$ with outer normal field $n$.}  \caption{Domain $\Omega=[0,1]^2$ with outer normal field $n$.}
39  \label{fig:FirstSteps.1}  \label{fig:FirstSteps.1}
# Line 128  problem\index{boundary value problem} (B Line 128  problem\index{boundary value problem} (B
128  the unknown function~$u$.  the unknown function~$u$.
129
130
131  \begin{figure}[h]  \begin{figure}[ht]
132  \centerline{\includegraphics[width=\figwidth]{figures/FirstStepMesh}}  \centerline{\includegraphics[width=\figwidth]{figures/FirstStepMesh}}
133  \caption{Mesh of $4 \time 4$ elements on a rectangular domain.  Here  \caption{Mesh of $4 \time 4$ elements on a rectangular domain.  Here
134  each element is a quadrilateral and described by four nodes, namely  each element is a quadrilateral and described by four nodes, namely
# Line 156  in the $x_1$ direction over the unit squ Line 156  in the $x_1$ direction over the unit squ
156  For more details we refer the reader to the literature, for instance \Ref{Zienc,NumHand}.  For more details we refer the reader to the literature, for instance \Ref{Zienc,NumHand}.
157
158  The \escript solver we want to use to solve this problem is embedded into the python interpreter language. So you can solve the problem interactively but you will learn quickly  The \escript solver we want to use to solve this problem is embedded into the python interpreter language. So you can solve the problem interactively but you will learn quickly
159  that is more efficient to use scripts which you can edit with your favourite editor.  that is more efficient to use scripts which you can edit with your favorite editor.
160  To enter the escript environment you use \program{escript} command\footnote{\program{escript} is not available under Windows yet. If you run under windows you can just use the  To enter the escript environment you use \program{escript} command\footnote{\program{escript} is not available under Windows yet. If you run under windows you can just use the
161  \program{python} command and the \env{OMP_NUM_THREADS} environment variable to control the number  \program{python} command and the \env{OMP_NUM_THREADS} environment variable to control the number
# Line 176  Here you can use all available python co Line 176  Here you can use all available python co
176  >>> print "2+3=",x  >>> print "2+3=",x
177  2+3= 5  2+3= 5
178  \end{python}  \end{python}
179  We refer to the python users guide if you not familar with python.  We refer to the python users guide if you not familiar with python.
180
181  \escript provides the class \Poisson to define a Poisson equation \index{Poisson equation}.  \escript provides the class \Poisson to define a Poisson equation \index{Poisson equation}.
182  (We will discuss a more general form of a PDE \index{partial differential equation!PDE}  (We will discuss a more general form of a PDE \index{partial differential equation!PDE}
# Line 308  Now we can write the script to solve our Line 308  Now we can write the script to solve our
308    mypde = Poisson(domain=mydomain)    mypde = Poisson(domain=mydomain)
310    u = mypde.getSolution()    u = mypde.getSolution()
# write u to an external file
saveVTK("u.xml",sol=u)
311  \end{python}  \end{python}
312  The entire code is available as \file{poisson.py} in the \ExampleDirectory  The question is what we do with the calculated solution \var{u}. Besides postprocessing, eg. calculating the gradient or the average value, which will be discussed later, plotting the solution is one one things you might want to do. \escript offers two ways to do this, both base on external modules or packages and so data need to converted
313    to hand over the solution. The first option is using the \MATPLOTLIB module which allows plotting 2D results relatively quickly, see~\cite{matplotlib}. However, there are limitations when using this tool, eg. in problem size and when solving 3D problems. Therefore \escript provides a second options based on \VTK files which is especially
314    designed for large scale and 3D problem and which can be read by a variety of software packages such as \mayavi \cite{mayavi}, \VisIt~\cite{VisIt}.
315
316  The last statement writes the solution (tagged with the name "sol") to a file named \file{u.xml} in  \subsection{Plotting Using \MATPLOTLIB}
317  \VTK file format.
318  Now you may run the script using the \escript environment  Users of Debian 5(Lenny) please note: this example makes use of the griddata method in \module{matplotlib.mlab}.
319  and visualize the solution using \mayavi:  This method is not part of version 0.98.1 which is available with Lenny.
320    If you wish to use \MATPLOTLIB, you may need to install a later version.
321    Users of Ubuntu 8.10 or later and people should be fine.
322
323    The \MATPLOTLIB module provides a simple and easy to use way to visualize PDE solutions (or other \Data objects).
324    To hand over data from \escript to \MATPLOTLIB the values need to mapped onto a rectangular grid. We will make use
325    of the \numpy module.
326
327    First we need to create a rectangular grid. We use the following statements:
328    \begin{python}
329    import numpy
330    x_grid = numpy.linspace(0.,1.,50)
331    y_grid = numpy.linspace(0.,1.,50)
332    \end{python}
333    \var{x_grid} is an array defining the x coordinates of the grids while
334    \var{y_grid} defines the y coordinates of the grid. In this case we use $50$ points over the interval $[0,1]$
335    in both directions.
336
337    Now the values created by \escript need to be interpolated to this grid. We will use the \MATPLOTLIB
338    \function{mlab.griddata} function to do this. We can easily extract spatial coordinates as a \var{list} by
339    \begin{python}
340    x=mydomain.getX().toListOfTuples()
341    y=mydomain.getX().toListOfTuples()
342    \end{python}
343    In principle we can apply the same \member{toListOfTuples} method to extract the values from the
344    PDE solution \var{u}. However, we have to make sure that the \Data object we extract the values from
345    uses the same \FunctionSpace as we have us when extracting \var{x} and \var{y}. By applying the
346    \function{interpolation} to \var{u} before extraction we can achieve this:
347    \begin{python}
348    z=interpolate(u,mydomain.getX().getFunctionSpace())
349    \end{python}
350    The values in \var{z} are now the values at the points with the coordinates given by \var{x} and \var{y}. These
351    values are now interpolated to the grid defined by \var{x_grid} and \var{y_grid} by using
352    \begin{python}
353    import matplotlib
354    z_grid = matplotlib.mlab.griddata(x,y,z,xi=x_grid,yi=y_grid )
355    \end{python}
356    \var{z_grid} gives now the values of the PDE solution \var{u} at the grid. The values can be plotted now
357    using the \function{contourf}:
358    \begin{python}
359    matplotlib.pyplot.contourf(x_grid, y_grid, z_grid, 5)
360    matplotlib.pyplot.savefig("u.png")
361    \end{python}
362    Here we use $5$ contours. The last statement writes the plot to the file \file{u.png} in the PNG format. Alternatively, one can use
363    \begin{python}
364    matplotlib.pyplot.contourf(x_grid, y_grid, z_grid, 5)
365    matplotlib.pyplot.show()
366    \end{python}
367    which gives an interactive browser window.
368
369    \begin{figure}
370    \centerline{\includegraphics[width=\figwidth]{figures/FirstStepResultMATPLOTLIB}}
371    \caption{Visualization of the Poisson Equation Solution for $f=1$ using \MATPLOTLIB.}
372    \label{fig:FirstSteps.3b}
373    \end{figure}
374
375    Now we can write the script to solve our Poisson problem
376    \begin{python}
377    from esys.escript import *
378    from esys.escript.linearPDEs import Poisson
379    from esys.finley import Rectangle
380    import numpy
381    import matplotlib
382    import pylab
383    # generate domain:
384    mydomain = Rectangle(l0=1.,l1=1.,n0=40, n1=20)
385    # define characteristic function of Gamma^D
386    x = mydomain.getX()
388    # define PDE and get its solution u
389    mypde = Poisson(domain=mydomain)
391    u = mypde.getSolution()
392    # interpolate u to a matplotlib grid:
393    x_grid = numpy.linspace(0.,1.,50)
394    y_grid = numpy.linspace(0.,1.,50)
395    x=mydomain.getX().toListOfTuples()
396    y=mydomain.getX().toListOfTuples()
397    z=interpolate(u,mydomain.getX().getFunctionSpace())
398    z_grid = matplotlib.mlab.griddata(x,y,z,xi=x_grid,yi=y_grid )
399    # interpolate u to a rectangular grid:
400    matplotlib.pyplot.contourf(x_grid, y_grid, z_grid, 5)
401    matplotlib.pyplot.savefig("u.png")
402    \end{python}
403    The entire code is available as \file{poisson\hackscore matplotlib.py} in the \ExampleDirectory.
404    You can run the script using the {\it escript} environment
405  \begin{verbatim}  \begin{verbatim}
406    escript poisson.py    escript poisson_matplotlib.py
mayavi -d u.xml -m SurfaceMap
407  \end{verbatim}  \end{verbatim}
408  See \fig{fig:FirstSteps.3}.  This will create the \file{u.png}, see Figure~\fig{fig:FirstSteps.3b}.
409    For details on the usage of the \MATPLOTLIB module we refer to the documentation~\cite{matplotlib}.
410
411    As pointed out, \MATPLOTLIB is restricted to the two-dimensional case and
412    should be used for small problems only. It can not be used under \MPI as the \member{toListOfTuples} method is
413    not safe under \MPI\footnote{The phrase 'safe under \MPI' means that a program will produce correct results when run on more than one processor under \MPI.}.
414
415    \subsection{Visualization using \VTK}
416
417    As an alternative {\it escript} supports the usage of visualization tools which base on \VTK, eg. mayavi \cite{mayavi}, \VisIt~\cite{VisIt}. In this case the solution is written to a file in the \VTK format. This file the can read by the tool of choice. Using \VTK file is \MPI safe.
418
419    To write the solution \var{u} in  Poisson problem to the file \file{u.xml} one need to add the line
420    \begin{python}
421    saveVTK("u.xml",sol=u)
422    \end{python}
423    The solution \var{u} is now available in the \file{u.xml} tagged with the name "sol".
424
425  \begin{figure}  \begin{figure}
426  \centerline{\includegraphics[width=\figwidth]{figures/FirstStepResult}}  \centerline{\includegraphics[width=\figwidth]{figures/FirstStepResult}}
# Line 329  See \fig{fig:FirstSteps.3}. Line 428  See \fig{fig:FirstSteps.3}.
428  \label{fig:FirstSteps.3}  \label{fig:FirstSteps.3}
429  \end{figure}  \end{figure}
430
431    The Poisson problem script is now
432    \begin{python}
433      from esys.escript import *
434      from esys.escript.linearPDEs import Poisson
435      from esys.finley import Rectangle
436      # generate domain:
437      mydomain = Rectangle(l0=1.,l1=1.,n0=40, n1=20)
438      # define characteristic function of Gamma^D
439      x = mydomain.getX()
441      # define PDE and get its solution u
442      mypde = Poisson(domain=mydomain)
444      u = mypde.getSolution()
445      # write u to an external file
446      saveVTK("u.xml",sol=u)
447    \end{python}
448    The entire code is available as \file{poisson\hackscore VTK.py} in the \ExampleDirectory
449
450    You can run the script using the {\it escript} environment
451    and visualize the solution using \mayavi:
452    \begin{verbatim}
453      escript poisson\hackscore VTK.py
454      mayavi2 -d u.xml -m SurfaceMap
455    \end{verbatim}
456    The result is shown in Figure~\fig{fig:FirstSteps.3}.

Legend:
 Removed from v.2370 changed lines Added in v.2575