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

revision 2363 by gross, Fri Apr 3 03:56:19 2009 UTC revision 2881 by jfenwick, Thu Jan 28 02:03:15 2010 UTC
# Line 1  Line 1
1
2  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
3  %  %
4  % Copyright (c) 2003-2008 by University of Queensland  % Copyright (c) 2003-2010 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 153  elements is called a mesh\index{finite e Line 153  elements is called a mesh\index{finite e
153  method!mesh}. \fig{fig:FirstSteps.2} shows an  method!mesh}. \fig{fig:FirstSteps.2} shows an
154  example of a FEM mesh with four elements in the $x_0$ and four elements  example of a FEM mesh with four elements in the $x_0$ and four elements
155  in the $x_1$ direction over the unit square.    in the $x_1$ direction over the unit square.
156  For more details we refer the reader to the literature, for instance  For more details we refer the reader to the literature, for instance \Ref{Zienc,NumHand}.
157  \Ref{Zienc,NumHand}.
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
159    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
161    \program{python} command and the \env{OMP_NUM_THREADS} environment variable to control the number
163    \begin{verbatim}
164      escript
165    \end{verbatim}
166    which will pass you on to the python prompt
167    \begin{verbatim}
168    Python 2.5.2 (r252:60911, Oct  5 2008, 19:29:17)
169    [GCC 4.3.2] on linux2
171    >>>
172    \end{verbatim}
173    Here you can use all available python commands and language features, for instance
174    \begin{python}
175     >>> x=2+3
176    >>> print "2+3=",x
177    2+3= 5
178    \end{python}
179    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 286  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    \subsection{Plotting Using \MATPLOTLIB}
317    The \MATPLOTLIB module provides a simple and easy to use way to visualize PDE solutions (or other \Data objects).
318    To hand over data from \escript to \MATPLOTLIB the values need to mapped onto a rectangular grid
319    \footnote{Users of Debian 5(Lenny) please note: this example makes use of the \function{griddata} method in \module{matplotlib.mlab}.
320    This method is not part of version 0.98.1 which is available with Lenny.
321    If you wish to use contour plots, you may need to install a later version.
322    Users of Ubuntu 8.10 or later should be fine.}. We will make use
323    of the \numpy module.
324
325    First we need to create a rectangular grid. We use the following statements:
326    \begin{python}
327    import numpy
328    x_grid = numpy.linspace(0.,1.,50)
329    y_grid = numpy.linspace(0.,1.,50)
330    \end{python}
331    \var{x_grid} is an array defining the x coordinates of the grids while
332    \var{y_grid} defines the y coordinates of the grid. In this case we use $50$ points over the interval $[0,1]$
333    in both directions.
334
335    Now the values created by \escript need to be interpolated to this grid. We will use the \MATPLOTLIB
336    \function{mlab.griddata} function to do this. We can easily extract spatial coordinates as a \var{list} by
337    \begin{python}
338    x=mydomain.getX().toListOfTuples()
339    y=mydomain.getX().toListOfTuples()
340    \end{python}
341    In principle we can apply the same \member{toListOfTuples} method to extract the values from the
342    PDE solution \var{u}. However, we have to make sure that the \Data object we extract the values from
343    uses the same \FunctionSpace as we have us when extracting \var{x} and \var{y}. We apply the
344    \function{interpolation} to \var{u} before extraction to achieve this:
345    \begin{python}
346    z=interpolate(u,mydomain.getX().getFunctionSpace())
347    \end{python}
348    The values in \var{z} are now the values at the points with the coordinates given by \var{x} and \var{y}. These
349    values are now interpolated to the grid defined by \var{x_grid} and \var{y_grid} by using
350    \begin{python}
351    import matplotlib
352    z_grid = matplotlib.mlab.griddata(x,y,z,xi=x_grid,yi=y_grid )
353    \end{python}
354    \var{z_grid} gives now the values of the PDE solution \var{u} at the grid. The values can be plotted now
355    using the \function{contourf}:
356    \begin{python}
357    matplotlib.pyplot.contourf(x_grid, y_grid, z_grid, 5)
358    matplotlib.pyplot.savefig("u.png")
359    \end{python}
360    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
361    \begin{python}
362    matplotlib.pyplot.contourf(x_grid, y_grid, z_grid, 5)
363    matplotlib.pyplot.show()
364    \end{python}
365    which gives an interactive browser window.
366
367  The last statement writes the solution (tagged with the name "sol") to a file named \file{u.xml} in  \begin{figure}
368  \VTK file format.  \centerline{\includegraphics[width=\figwidth]{figures/FirstStepResultMATPLOTLIB}}
369  Now you may run the script and visualize the solution using \mayavi:  \caption{Visualization of the Poisson Equation Solution for $f=1$ using \MATPLOTLIB.}
370    \label{fig:FirstSteps.3b}
371    \end{figure}
372
373    Now we can write the script to solve our Poisson problem
374    \begin{python}
375    from esys.escript import *
376    from esys.escript.linearPDEs import Poisson
377    from esys.finley import Rectangle
378    import numpy
379    import matplotlib
380    import pylab
381    # generate domain:
382    mydomain = Rectangle(l0=1.,l1=1.,n0=40, n1=20)
383    # define characteristic function of Gamma^D
384    x = mydomain.getX()
386    # define PDE and get its solution u
387    mypde = Poisson(domain=mydomain)
389    u = mypde.getSolution()
390    # interpolate u to a matplotlib grid:
391    x_grid = numpy.linspace(0.,1.,50)
392    y_grid = numpy.linspace(0.,1.,50)
393    x=mydomain.getX().toListOfTuples()
394    y=mydomain.getX().toListOfTuples()
395    z=interpolate(u,mydomain.getX().getFunctionSpace())
396    z_grid = matplotlib.mlab.griddata(x,y,z,xi=x_grid,yi=y_grid )
397    # interpolate u to a rectangular grid:
398    matplotlib.pyplot.contourf(x_grid, y_grid, z_grid, 5)
399    matplotlib.pyplot.savefig("u.png")
400    \end{python}
401    The entire code is available as \file{poisson\hackscore matplotlib.py} in the \ExampleDirectory.
402    You can run the script using the {\it escript} environment
403  \begin{verbatim}  \begin{verbatim}
404    python poisson.py    escript poisson_matplotlib.py
mayavi -d u.xml -m SurfaceMap
405  \end{verbatim}  \end{verbatim}
406  See \fig{fig:FirstSteps.3}.  This will create the \file{u.png}, see Figure~\fig{fig:FirstSteps.3b}.
407    For details on the usage of the \MATPLOTLIB module we refer to the documentation~\cite{matplotlib}.
408
409    As pointed out, \MATPLOTLIB is restricted to the two-dimensional case and
410    should be used for small problems only. It can not be used under \MPI as the \member{toListOfTuples} method is
411    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.}.
412
413  \begin{figure}  \begin{figure}
414  \centerline{\includegraphics[width=\figwidth]{figures/FirstStepResult}}  \centerline{\includegraphics[width=\figwidth]{figures/FirstStepResult}}
# Line 306  See \fig{fig:FirstSteps.3}. Line 416  See \fig{fig:FirstSteps.3}.
416  \label{fig:FirstSteps.3}  \label{fig:FirstSteps.3}
417  \end{figure}  \end{figure}
418
419    \subsection{Visualization using \VTK}
420
421    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.
422
423    To write the solution \var{u} in  Poisson problem to the file \file{u.xml} one need to add the line
424    \begin{python}
425    saveVTK("u.xml",sol=u)
426    \end{python}
427    The solution \var{u} is now available in the \file{u.xml} tagged with the name "sol".
428
429    The Poisson problem script is now
430    \begin{python}
431      from esys.escript import *
432      from esys.escript.linearPDEs import Poisson
433      from esys.finley import Rectangle
434      # generate domain:
435      mydomain = Rectangle(l0=1.,l1=1.,n0=40, n1=20)
436      # define characteristic function of Gamma^D
437      x = mydomain.getX()
439      # define PDE and get its solution u
440      mypde = Poisson(domain=mydomain)
442      u = mypde.getSolution()
443      # write u to an external file
444      saveVTK("u.xml",sol=u)
445    \end{python}
446    The entire code is available as \file{poisson\hackscore VTK.py} in the \ExampleDirectory.
447
448    You can run the script using the {\it escript} environment
449    and visualize the solution using \mayavi:
450    \begin{verbatim}
451      escript poisson\hackscore VTK.py
452      mayavi2 -d u.xml -m SurfaceMap
453    \end{verbatim}
454    The result is shown in Figure~\fig{fig:FirstSteps.3}.

Legend:
 Removed from v.2363 changed lines Added in v.2881