--- temp_trunk_copy/doc/user/firststep.tex 2008/01/11 02:29:38 1384
+++ trunk/doc/user/firststep.tex 2009/08/03 23:38:33 2575
@@ -1,27 +1,21 @@
+
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
-% $Id$
-%
-%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
-%
-% Copyright 2003-2007 by ACceSS MNRF
-% Copyright 2007 by University of Queensland
-%
-% http://esscc.uq.edu.au
-% Primary Business: Queensland, Australia
-% Licensed under the Open Software License version 3.0
-% http://www.opensource.org/licenses/osl-3.0.php
+% Copyright (c) 2003-2009 by University of Queensland
+% Earth Systems Science Computational Center (ESSCC)
+% http://www.uq.edu.au/esscc
%
-%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+% Primary Business: Queensland, Australia
+% Licensed under the Open Software License version 3.0
+% http://www.opensource.org/licenses/osl-3.0.php
%
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+
\section{The First Steps}
\label{FirstSteps}
-\begin{figure}
-\centerline{\includegraphics[width=\figwidth]{figures/FirstStepDomain}}
-\caption{Domain $\Omega=[0,1]^2$ with outer normal field $n$.}
-\label{fig:FirstSteps.1}
-\end{figure}
+
In this chapter we give an introduction how to use \escript to solve
a partial differential equation \index{partial differential equation} (PDE \index{partial differential equation!PDE}). We assume you are at least a little familiar with Python. The knowledge presented at the Python tutorial at \url{http://docs.python.org/tut/tut.html}
@@ -39,6 +33,11 @@
\label{eq:FirstSteps.1b}
\end{equation}
The domain is shown in \fig{fig:FirstSteps.1}.
+\begin{figure} [ht]
+\centerline{\includegraphics[width=\figwidth]{figures/FirstStepDomain}}
+\caption{Domain $\Omega=[0,1]^2$ with outer normal field $n$.}
+\label{fig:FirstSteps.1}
+\end{figure}
$\Delta$ denotes the Laplace operator\index{Laplace operator}, which is defined by
\begin{equation}
@@ -60,7 +59,7 @@
\end{equation*}
}
Basically, in the subindex of a function, any index to the left of the comma denotes a spatial derivative with respect
-to the index. To get a more compact form we will write $w\hackscore{,ij}=(w\hackscore {,i})\hackscore{,j}$
+to the index. To get a more compact form we will write $u\hackscore{,ij}=(u\hackscore {,i})\hackscore{,j}$
which leads to
\begin{equation}
\Delta u = u\hackscore{,00}+u\hackscore{,11}=\sum\hackscore{i=0}^2 u\hackscore{,ii}
@@ -99,9 +98,9 @@
$n\hackscore{i} u\hackscore{,i}= n\hackscore{0} u\hackscore{,0} +
n\hackscore{1} u\hackscore{,1}$.
\footnote{Some readers may familiar with the notation
-\begin{equation*}
+$
\frac{\partial u}{\partial n} = n\hackscore{i} u\hackscore{,i}
-\end{equation*}
+$
for the normal derivative.}
The Neumann boundary condition of \eqn{eq:FirstSteps.2} should be fulfilled on the
set $\Gamma^N$ which is the top and right edge of the domain:
@@ -126,12 +125,11 @@
Dirichlet boundary condition in \eqn{eq:FirstSteps.2d} form a so
called boundary value
problem\index{boundary value problem} (BVP\index{boundary value problem!BVP}) for
-the unknown
-function $u$.
+the unknown function~$u$.
-\begin{figure}
-\centerline{\includegraphics[width=\figwidth]{figures/FirstStepMesh.eps}}
+\begin{figure}[ht]
+\centerline{\includegraphics[width=\figwidth]{figures/FirstStepMesh}}
\caption{Mesh of $4 \time 4$ elements on a rectangular domain. Here
each element is a quadrilateral and described by four nodes, namely
the corner points. The solution is interpolated by a bi-linear
@@ -155,8 +153,30 @@
method!mesh}. \fig{fig:FirstSteps.2} shows an
example of a FEM mesh with four elements in the $x_0$ and four elements
in the $x_1$ direction over the unit square.
-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}.
+
+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
+that is more efficient to use scripts which you can edit with your favorite editor.
+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
+\program{python} command and the \env{OMP_NUM_THREADS} environment variable to control the number
+of threads.}:
+\begin{verbatim}
+ escript
+\end{verbatim}
+which will pass you on to the python prompt
+\begin{verbatim}
+Python 2.5.2 (r252:60911, Oct 5 2008, 19:29:17)
+[GCC 4.3.2] on linux2
+Type "help", "copyright", "credits" or "license" for more information.
+>>>
+\end{verbatim}
+Here you can use all available python commands and language features, for instance
+\begin{python}
+ >>> x=2+3
+>>> print "2+3=",x
+2+3= 5
+\end{python}
+We refer to the python users guide if you not familiar with python.
\escript provides the class \Poisson to define a Poisson equation \index{Poisson equation}.
(We will discuss a more general form of a PDE \index{partial differential equation!PDE}
@@ -250,7 +270,7 @@
Similarly, \code{whereZero(x[1])} creates function which equals $1$ where \code{x[1]} is
equal to zero and $0$ elsewhere.
The sum of the results of \code{whereZero(x[0])} and \code{whereZero(x[1])}
-gives a function on the domain \var{mydomain} which is exactly positive where $x\hackscore{0}$ or $x\hackscore{1}$ is equal to zero.
+gives a function on the domain \var{mydomain} which is strictly positive where $x\hackscore{0}$ or $x\hackscore{1}$ is equal to zero.
Note that \var{gammaD} has the same \Rank, \Shape and \FunctionSpace like \var{x0} used to define it. So from
\begin{python}
print gammaD.getRank(),gammaD.getShape(),gammaD.getFunctionSpace()
@@ -288,23 +308,149 @@
mypde = Poisson(domain=mydomain)
mypde.setValue(f=1,q=gammaD)
u = mypde.getSolution()
- # write u to an external file
- saveVTK("u.xml",sol=u)
\end{python}
-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
+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
+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}.
+
+\subsection{Plotting Using \MATPLOTLIB}
+
+Users of Debian 5(Lenny) please note: this example makes use of the griddata method in \module{matplotlib.mlab}.
+This method is not part of version 0.98.1 which is available with Lenny.
+If you wish to use \MATPLOTLIB, you may need to install a later version.
+Users of Ubuntu 8.10 or later and people should be fine.
-The last statement writes the solution (tagged with the name "sol") to a file named \file{u.xml} in
-\VTK file format.
-Now you may run the script and visualize the solution using \mayavi:
+The \MATPLOTLIB module provides a simple and easy to use way to visualize PDE solutions (or other \Data objects).
+To hand over data from \escript to \MATPLOTLIB the values need to mapped onto a rectangular grid. We will make use
+of the \numpy module.
+
+First we need to create a rectangular grid. We use the following statements:
+\begin{python}
+import numpy
+x_grid = numpy.linspace(0.,1.,50)
+y_grid = numpy.linspace(0.,1.,50)
+\end{python}
+\var{x_grid} is an array defining the x coordinates of the grids while
+\var{y_grid} defines the y coordinates of the grid. In this case we use $50$ points over the interval $[0,1]$
+in both directions.
+
+Now the values created by \escript need to be interpolated to this grid. We will use the \MATPLOTLIB
+\function{mlab.griddata} function to do this. We can easily extract spatial coordinates as a \var{list} by
+\begin{python}
+x=mydomain.getX()[0].toListOfTuples()
+y=mydomain.getX()[1].toListOfTuples()
+\end{python}
+In principle we can apply the same \member{toListOfTuples} method to extract the values from the
+PDE solution \var{u}. However, we have to make sure that the \Data object we extract the values from
+uses the same \FunctionSpace as we have us when extracting \var{x} and \var{y}. By applying the
+\function{interpolation} to \var{u} before extraction we can achieve this:
+\begin{python}
+z=interpolate(u,mydomain.getX().getFunctionSpace())
+\end{python}
+The values in \var{z} are now the values at the points with the coordinates given by \var{x} and \var{y}. These
+values are now interpolated to the grid defined by \var{x_grid} and \var{y_grid} by using
+\begin{python}
+import matplotlib
+z_grid = matplotlib.mlab.griddata(x,y,z,xi=x_grid,yi=y_grid )
+\end{python}
+\var{z_grid} gives now the values of the PDE solution \var{u} at the grid. The values can be plotted now
+using the \function{contourf}:
+\begin{python}
+matplotlib.pyplot.contourf(x_grid, y_grid, z_grid, 5)
+matplotlib.pyplot.savefig("u.png")
+\end{python}
+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
+\begin{python}
+matplotlib.pyplot.contourf(x_grid, y_grid, z_grid, 5)
+matplotlib.pyplot.show()
+\end{python}
+which gives an interactive browser window.
+
+\begin{figure}
+\centerline{\includegraphics[width=\figwidth]{figures/FirstStepResultMATPLOTLIB}}
+\caption{Visualization of the Poisson Equation Solution for $f=1$ using \MATPLOTLIB.}
+\label{fig:FirstSteps.3b}
+\end{figure}
+
+Now we can write the script to solve our Poisson problem
+\begin{python}
+from esys.escript import *
+from esys.escript.linearPDEs import Poisson
+from esys.finley import Rectangle
+import numpy
+import matplotlib
+import pylab
+# generate domain:
+mydomain = Rectangle(l0=1.,l1=1.,n0=40, n1=20)
+# define characteristic function of Gamma^D
+x = mydomain.getX()
+gammaD = whereZero(x[0])+whereZero(x[1])
+# define PDE and get its solution u
+mypde = Poisson(domain=mydomain)
+mypde.setValue(f=1,q=gammaD)
+u = mypde.getSolution()
+# interpolate u to a matplotlib grid:
+x_grid = numpy.linspace(0.,1.,50)
+y_grid = numpy.linspace(0.,1.,50)
+x=mydomain.getX()[0].toListOfTuples()
+y=mydomain.getX()[1].toListOfTuples()
+z=interpolate(u,mydomain.getX().getFunctionSpace())
+z_grid = matplotlib.mlab.griddata(x,y,z,xi=x_grid,yi=y_grid )
+# interpolate u to a rectangular grid:
+matplotlib.pyplot.contourf(x_grid, y_grid, z_grid, 5)
+matplotlib.pyplot.savefig("u.png")
+\end{python}
+The entire code is available as \file{poisson\hackscore matplotlib.py} in the \ExampleDirectory.
+You can run the script using the {\it escript} environment
\begin{verbatim}
- python poisson.py
- mayavi -d u.xml -m SurfaceMap
+ escript poisson_matplotlib.py
\end{verbatim}
-See \fig{fig:FirstSteps.3}.
+This will create the \file{u.png}, see Figure~\fig{fig:FirstSteps.3b}.
+For details on the usage of the \MATPLOTLIB module we refer to the documentation~\cite{matplotlib}.
+
+As pointed out, \MATPLOTLIB is restricted to the two-dimensional case and
+should be used for small problems only. It can not be used under \MPI as the \member{toListOfTuples} method is
+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.}.
+
+\subsection{Visualization using \VTK}
+
+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.
+
+To write the solution \var{u} in Poisson problem to the file \file{u.xml} one need to add the line
+\begin{python}
+saveVTK("u.xml",sol=u)
+\end{python}
+The solution \var{u} is now available in the \file{u.xml} tagged with the name "sol".
\begin{figure}
-\centerline{\includegraphics[width=\figwidth]{figures/FirstStepResult.eps}}
+\centerline{\includegraphics[width=\figwidth]{figures/FirstStepResult}}
\caption{Visualization of the Poisson Equation Solution for $f=1$}
\label{fig:FirstSteps.3}
\end{figure}
+The Poisson problem script is now
+\begin{python}
+ from esys.escript import *
+ from esys.escript.linearPDEs import Poisson
+ from esys.finley import Rectangle
+ # generate domain:
+ mydomain = Rectangle(l0=1.,l1=1.,n0=40, n1=20)
+ # define characteristic function of Gamma^D
+ x = mydomain.getX()
+ gammaD = whereZero(x[0])+whereZero(x[1])
+ # define PDE and get its solution u
+ mypde = Poisson(domain=mydomain)
+ mypde.setValue(f=1,q=gammaD)
+ u = mypde.getSolution()
+ # write u to an external file
+ saveVTK("u.xml",sol=u)
+\end{python}
+The entire code is available as \file{poisson\hackscore VTK.py} in the \ExampleDirectory
+
+You can run the script using the {\it escript} environment
+and visualize the solution using \mayavi:
+\begin{verbatim}
+ escript poisson\hackscore VTK.py
+ mayavi2 -d u.xml -m SurfaceMap
+\end{verbatim}
+The result is shown in Figure~\fig{fig:FirstSteps.3}.