/[escript]/trunk/doc/user/firststep.tex
ViewVC logotype

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

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

revision 2375 by gross, Tue Apr 7 04:41:28 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 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  
   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()[0].toListOfTuples()
341    y=mydomain.getX()[1].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()
387    gammaD = whereZero(x[0])+whereZero(x[1])
388    # define PDE and get its solution u
389    mypde = Poisson(domain=mydomain)
390    mypde.setValue(f=1,q=gammaD)
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()[0].toListOfTuples()
396    y=mydomain.getX()[1].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()
440      gammaD = whereZero(x[0])+whereZero(x[1])
441      # define PDE and get its solution u
442      mypde = Poisson(domain=mydomain)
443      mypde.setValue(f=1,q=gammaD)
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.2375  
changed lines
  Added in v.2575

  ViewVC Help
Powered by ViewVC 1.1.26