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

revision 2548 by jfenwick, Mon Jul 20 06:20:06 2009 UTC revision 2574 by gross, Mon Aug 3 09:09:20 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)
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.  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. We will make use
319  and visualize the solution using \mayavi:  of the \numpy module.
320
321    First we need to create a rectangular grid. We use the following statements:
322    \begin{python}
323    import numpy
324    x_grid = numpy.linspace(0.,1.,50)
325    y_grid = numpy.linspace(0.,1.,50)
326    \end{python}
327    \var{x_grid} is an array defining the x coordinates of the grids while
328    \var{y_grid} defines the y coordinates of the grid. In this case we use $50$ points over the interval $[0,1]$
329    in both directions.
330
331    Now created by \escript need to be interpolated to this grid. We will use the \MATPLOTLIB
332    \function{mlab.griddata} function to do this. We can easily extract spatial coordinates as a \var{list} by
333    \begin{python}
334    x=mydomain.getX().toListOfTuples()
335    y=mydomain.getX().toListOfTuples()
336    \end{python}
337    In principle we can apply the same \member{toListOfTuples} method to extract the values from the
338    PDE solution \var{u}. However, we have to make sure that the \Data object we extract the values from
339    uses the same \FunctionSpace as we have us when extracting \var{x} and \var{y}. By applying the
340    \function{interpolation} to \var{u} before extraction we can achieve this:
341    \begin{python}
342    z=interpolate(u,mydomain.getX().getFunctionSpace())
343    \end{python}
344    The values in \var{z} are now the values at the points with the coordinates given by \var{x} and \var{y}. These
345    values are now interpolated to the grid defined by \var{x_grid} and \var{y_grid} by using
346    \begin{python}
347    import matplotlib
348    z_grid = matplotlib.mlab.griddata(x,y,z,xi=x_grid,yi=y_grid )
349    \end{python}
350    \var{z_grid} gives now the values of the PDE solution \var{u} at the grid. The values can be plotted now
351    using the \function{contourf}:
352    \begin{python}
353    matplotlib.pyplot.contourf(x_grid, y_grid, z_grid, 5)
354    matplotlib.pyplot.savefig("u.png")
355    \end{python}
356    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
357    \begin{python}
358    matplotlib.pyplot.contourf(x_grid, y_grid, z_grid, 5)
359    matplotlib.pyplot.show()
360    \end{python}
361    which gives an interactive browser window.
362
363    \begin{figure}
364    \centerline{\includegraphics[width=\figwidth]{figures/FirstStepResultMATPLOTLIB}}
365    \caption{Visualization of the Poisson Equation Solution for $f=1$ using \MATPLOTLIB.}
366    \label{fig:FirstSteps.3b}
367    \end{figure}
368
369    Now we can write the script to solve our Poisson problem
370    \begin{python}
371    from esys.escript import *
372    from esys.escript.linearPDEs import Poisson
373    from esys.finley import Rectangle
374    import numpy
375    import matplotlib
376    import pylab
377    # generate domain:
378    mydomain = Rectangle(l0=1.,l1=1.,n0=40, n1=20)
379    # define characteristic function of Gamma^D
380    x = mydomain.getX()
381    gammaD = whereZero(x)+whereZero(x)
382    # define PDE and get its solution u
383    mypde = Poisson(domain=mydomain)
385    u = mypde.getSolution()
386    # interpolate u to a matplotlib grid:
387    x_grid = numpy.linspace(0.,1.,50)
388    y_grid = numpy.linspace(0.,1.,50)
389    x=mydomain.getX().toListOfTuples()
390    y=mydomain.getX().toListOfTuples()
391    z=interpolate(u,mydomain.getX().getFunctionSpace())
392    z_grid = matplotlib.mlab.griddata(x,y,z,xi=x_grid,yi=y_grid )
393    # interpolate u to a rectangular grid:
394    matplotlib.pyplot.contourf(x_grid, y_grid, z_grid, 5)
395    matplotlib.pyplot.savefig("u.png")
396    \end{python}
397    The entire code is available as \file{poisson\hackscore matplotlib.py} in the \ExampleDirectory.
398    You can run the script using the {\it escript} environment
399  \begin{verbatim}  \begin{verbatim}
400    escript poisson.py    escript poisson_matplotlib.py
mayavi -d u.xml -m SurfaceMap
401  \end{verbatim}  \end{verbatim}
402  See \fig{fig:FirstSteps.3}.  This will create the \file{u.png}, see Figure~\fig{fig:FirstSteps.3b}.
403    For details on the usage of the \MATPLOTLIB module we refer to the documentation~\cite{matplotlib}.
404
405    As pointed out, \MATPLOTLIB is restricted to the two-dimensional case and
406    should be used for small problems only. It can not be used under \MPI as the \member{toListOfTuples} method is
407    not safe under \MPI\footnote{The term 'safe under \MPI' means that a program will produce correct results when run on more than one processors under \MPI.}.
408
409    \subsection{Visualization using \VTK}
410
411    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.
412
413    To write the solution \var{u} in  Poisson problem to the file \file{u.xml} one need to add the line
414    \begin{python}
415    saveVTK("u.xml",sol=u)
416    \end{python}
417    The solution \var{u} is now available in the \file{u.xml} tagged with the name "sol".
418
419  \begin{figure}  \begin{figure}
420  \centerline{\includegraphics[width=\figwidth]{figures/FirstStepResult}}  \centerline{\includegraphics[width=\figwidth]{figures/FirstStepResult}}
# Line 329  See \fig{fig:FirstSteps.3}. Line 422  See \fig{fig:FirstSteps.3}.
422  \label{fig:FirstSteps.3}  \label{fig:FirstSteps.3}
423  \end{figure}  \end{figure}
424
425    The Poisson problem script is now
426    \begin{python}
427      from esys.escript import *
428      from esys.escript.linearPDEs import Poisson
429      from esys.finley import Rectangle
430      # generate domain:
431      mydomain = Rectangle(l0=1.,l1=1.,n0=40, n1=20)
432      # define characteristic function of Gamma^D
433      x = mydomain.getX()
434      gammaD = whereZero(x)+whereZero(x)
435      # define PDE and get its solution u
436      mypde = Poisson(domain=mydomain)
438      u = mypde.getSolution()
439      # write u to an external file
440      saveVTK("u.xml",sol=u)
441    \end{python}
442    The entire code is available as \file{poisson\hackscore VTK.py} in the \ExampleDirectory
443
444    You can run the script using the {\it escript} environment
445    and visualize the solution using \mayavi:
446    \begin{verbatim}
447      escript poisson\hackscore VTK.py
448      mayavi2 -d u.xml -m SurfaceMap
449    \end{verbatim}
450    The result is shown in Figure~\fig{fig:FirstSteps.3}.

Legend:
 Removed from v.2548 changed lines Added in v.2574