ViewVC logotype

Diff of /branches/stage3.0/doc/user/firststep.tex

Parent Directory Parent Directory | Revision Log Revision Log | View Patch Patch

revision 2583 by jfenwick, Fri Jul 31 05:39:36 2009 UTC revision 2584 by jfenwick, Wed Aug 5 02:44:51 2009 UTC
# 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)
309    mypde.setValue(f=1,q=gammaD)    mypde.setValue(f=1,q=gammaD)
310    u = mypde.getSolution()    u = mypde.getSolution()
   # write u to an external file  
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}.
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.  The \MATPLOTLIB module provides a simple and easy to use way to visualize PDE solutions (or other \Data objects).
318  Now you may run the script using the \escript environment  To hand over data from \escript to \MATPLOTLIB the values need to mapped onto a rectangular grid
319  and visualize the solution using \mayavi:  \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.
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.
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()[0].toListOfTuples()
339    y=mydomain.getX()[1].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.
367    \begin{figure}
368    \centerline{\includegraphics[width=\figwidth]{figures/FirstStepResultMATPLOTLIB}}
369    \caption{Visualization of the Poisson Equation Solution for $f=1$ using \MATPLOTLIB.}
370    \label{fig:FirstSteps.3b}
371    \end{figure}
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()
385    gammaD = whereZero(x[0])+whereZero(x[1])
386    # define PDE and get its solution u
387    mypde = Poisson(domain=mydomain)
388    mypde.setValue(f=1,q=gammaD)
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()[0].toListOfTuples()
394    y=mydomain.getX()[1].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    escript 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}.
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.}.
413  \begin{figure}  \begin{figure}
414  \centerline{\includegraphics[width=\figwidth]{figures/FirstStepResult}}  \centerline{\includegraphics[width=\figwidth]{figures/FirstStepResult}}
# Line 329  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}
419    \subsection{Visualization using \VTK}
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.
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".
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()
438      gammaD = whereZero(x[0])+whereZero(x[1])
439      # define PDE and get its solution u
440      mypde = Poisson(domain=mydomain)
441      mypde.setValue(f=1,q=gammaD)
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
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}.

Removed from v.2583  
changed lines
  Added in v.2584

  ViewVC Help
Powered by ViewVC 1.1.26