/[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 3277 by caltinay, Thu Oct 14 06:14:30 2010 UTC revision 3278 by caltinay, Fri Oct 15 02:00:45 2010 UTC
# Line 200  the right hand side $f$ of the PDE to co Line 200  the right hand side $f$ of the PDE to co
200    mypde = Poisson(mydomain)    mypde = Poisson(mydomain)
201    mypde.setValue(f=1)    mypde.setValue(f=1)
202  \end{python}  \end{python}
203  We have not specified any boundary condition but the  We have not specified any boundary condition but the \Poisson class implicitly
204  \Poisson class implicitly assumes homogeneous Neuman boundary conditions \index{Neumann  assumes homogeneous Neuman boundary conditions\index{Neumann boundary condition!homogeneous} defined by \eqn{eq:FirstSteps.2}.
205  boundary condition!homogeneous} defined by \eqn{eq:FirstSteps.2}. With this boundary  With this boundary condition the BVP\index{boundary value problem!BVP} we have
206  condition the BVP\index{boundary value problem!BVP} we have defined has no unique solution. In fact, with any solution $u$  defined has no unique solution.
207  and any constant $C$ the function $u+C$ becomes a solution as well. We have to add  In fact, with any solution $u$ and any constant $C$ the function $u+C$ becomes
208  a Dirichlet boundary condition \index{Dirichlet boundary condition}. This is done  a solution as well.
209  by defining a characteristic function \index{characteristic function}  We have to add a Dirichlet boundary condition\index{Dirichlet boundary condition}.
210  which has positive values at locations $x=(x\hackscore{0},x\hackscore{1})$ where Dirichlet boundary condition is set  This is done by defining a characteristic function\index{characteristic function}
211  and $0$ elsewhere. In our case of $\Gamma^D$ defined by \eqn{eq:FirstSteps.2c},  which has positive values at locations $x=(x\hackscore{0},x\hackscore{1})$
212  we need to construct a function \var{gammaD} which is positive for the cases $x\hackscore{0}=0$ or $x\hackscore{1}=0$. To get  where Dirichlet boundary condition is set and $0$ elsewhere.
213  an object \var{x} which contains the coordinates of the nodes in the domain use  In our case of $\Gamma^D$ defined by \eqn{eq:FirstSteps.2c}, we need to
214    construct a function \var{gammaD} which is positive for the cases $x\hackscore{0}=0$ or $x\hackscore{1}=0$.
215    To get an object \var{x} which contains the coordinates of the nodes in the domain use
216  \begin{python}  \begin{python}
217    x=mydomain.getX()    x=mydomain.getX()
218  \end{python}  \end{python}
219  The method \method{getX} of the \Domain \var{mydomain}  The method \method{getX} of the \Domain \var{mydomain} gives access to locations
220  gives access to locations  in the domain defined by \var{mydomain}.
221  in the domain defined by \var{mydomain}. The object \var{x} is actually a \Data object which will be  The object \var{x} is actually a \Data object which will be discussed in
222  discussed in \Chap{ESCRIPT CHAP} in more detail. What we need to know here is that  \Chap{ESCRIPT CHAP} in more detail.
223    What we need to know here is that \var{x} has \Rank (number of dimensions) and
224  \var{x} has \Rank (number of dimensions) and a \Shape (list of dimensions) which can be viewed by  a \Shape (list of dimensions) which can be viewed by calling the \method{getRank} and \method{getShape} methods:
 calling the \method{getRank} and \method{getShape} methods:  
225  \begin{python}  \begin{python}
226    print "rank ",x.getRank(),", shape ",x.getShape()    print "rank ",x.getRank(),", shape ",x.getShape()
227  \end{python}  \end{python}
# Line 237  will print Line 238  will print
238  \begin{python}  \begin{python}
239    Function space type: Finley_Nodes on FinleyMesh    Function space type: Finley_Nodes on FinleyMesh
240  \end{python}  \end{python}
241  which tells us that the coordinates are stored on the nodes of (rather than on points in the interior of) a \finley mesh.  which tells us that the coordinates are stored on the nodes of (rather than on
242  To get the  $x\hackscore{0}$ coordinates of the locations we use the  points in the interior of) a \finley mesh.
243  statement  To get the  $x\hackscore{0}$ coordinates of the locations we use the statement
244  \begin{python}  \begin{python}
245    x0=x[0]    x0=x[0]
246  \end{python}  \end{python}
247  Object \var{x0}  Object \var{x0} is again a \Data object now with \Rank $0$ and \Shape $()$.
248  is again a \Data object now with \Rank $0$ and  It inherits the \FunctionSpace from \var{x}:
 \Shape $()$. It inherits the \FunctionSpace from \var{x}:  
249  \begin{python}  \begin{python}
250    print x0.getRank(),x0.getShape(),x0.getFunctionSpace()    print x0.getRank(), x0.getShape(), x0.getFunctionSpace()
251  \end{python}  \end{python}
252  will print  will print
253  \begin{python}  \begin{python}
254    0 () Function space type: Finley_Nodes on FinleyMesh    0 () Function space type: Finley_Nodes on FinleyMesh
255  \end{python}  \end{python}
256  We can now construct a function \var{gammaD} which is only non-zero on the bottom and left edges  We can now construct a function \var{gammaD} which is only non-zero on the
257  of the domain with  bottom and left edges of the domain with
258  \begin{python}  \begin{python}
259    from esys.escript import whereZero    from esys.escript import whereZero
260    gammaD=whereZero(x[0])+whereZero(x[1])    gammaD=whereZero(x[0])+whereZero(x[1])
261  \end{python}  \end{python}
262    
263  \code{whereZero(x[0])} creates function which equals $1$ where \code{x[0]} is (almost) equal to zero  \code{whereZero(x[0])} creates a function which equals $1$ where \code{x[0]} is (almost) equal to zero and $0$ elsewhere.
264  and $0$ elsewhere.  Similarly, \code{whereZero(x[1])} creates a function which equals $1$ where \code{x[1]} is equal to zero and $0$ elsewhere.
265  Similarly, \code{whereZero(x[1])} creates function which equals $1$ where \code{x[1]} is  The sum of the results of \code{whereZero(x[0])} and \code{whereZero(x[1])}
 equal to zero and $0$ elsewhere.  
 The sum of the results of \code{whereZero(x[0])} and \code{whereZero(x[1])}  
266  gives a function on the domain \var{mydomain} which is strictly 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.
267  Note that \var{gammaD} has the same \Rank, \Shape and \FunctionSpace like \var{x0} used to define it. So from  Note that \var{gammaD} has the same \Rank, \Shape and \FunctionSpace like \var{x0} used to define it.
268    So from
269  \begin{python}  \begin{python}
270    print gammaD.getRank(),gammaD.getShape(),gammaD.getFunctionSpace()    print gammaD.getRank(), gammaD.getShape(), gammaD.getFunctionSpace()
271  \end{python}  \end{python}
272  one gets  one gets
273  \begin{python}  \begin{python}
274    0 () Function space type: Finley_Nodes on FinleyMesh    0 () Function space type: Finley_Nodes on FinleyMesh
275  \end{python}  \end{python}
276  An additional parameter \var{q} of the \code{setValue} method of the \Poisson class defines the  An additional parameter \var{q} of the \code{setValue} method of the \Poisson
277  characteristic function \index{characteristic function} of the locations  class defines the characteristic function\index{characteristic function} of
278  of the domain where homogeneous Dirichlet boundary condition \index{Dirichlet boundary condition!homogeneous}  the locations of the domain where the homogeneous Dirichlet boundary condition\index{Dirichlet boundary condition!homogeneous} is set.
279  are set. The complete definition of our example is now:  The complete definition of our example is now:
280  \begin{python}  \begin{python}
281    from esys.linearPDEs import Poisson    from esys.escript.linearPDEs import Poisson
282    x = mydomain.getX()    x = mydomain.getX()
283    gammaD = whereZero(x[0])+whereZero(x[1])    gammaD = whereZero(x[0])+whereZero(x[1])
284    mypde = Poisson(domain=mydomain)    mypde = Poisson(domain=mydomain)
285    mypde.setValue(f=1,q=gammaD)    mypde.setValue(f=1,q=gammaD)
286  \end{python}  \end{python}
287  The first statement imports the \Poisson class definition from the \linearPDEs module \escript package.  The first statement imports the \Poisson class definition from the \linearPDEs module.
288  To get the solution of the Poisson equation defined by \var{mypde} we just have to call its  To get the solution of the Poisson equation defined by \var{mypde} we just have to call its \method{getSolution} method.
 \method{getSolution}.  
289    
290  Now we can write the script to solve our Poisson problem  Now we can write the script to solve our Poisson problem
291  \begin{python}  \begin{python}
# Line 304  Now we can write the script to solve our Line 302  Now we can write the script to solve our
302    mypde.setValue(f=1,q=gammaD)    mypde.setValue(f=1,q=gammaD)
303    u = mypde.getSolution()    u = mypde.getSolution()
304  \end{python}  \end{python}
305  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  The question is what we do with the calculated solution \var{u}.
306  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  Besides postprocessing, e.g. calculating the gradient or the average value, which will be discussed later, plotting the solution is one of the things you might want to do.
307  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}.  \escript offers two ways to do this, both based on external modules or packages and so data need to converted to hand over the solution.
308    The first option is using the \MATPLOTLIB module which allows plotting 2D results relatively quickly from within the Python script, see~\cite{matplotlib}.
309    However, there are limitations when using this tool, especially for large problems and when solving 3-dimensional problems.
310    Therefore, \escript provides functionality to export data as files which can subsequently be read by third-party software packages such as \mayavi\cite{mayavi} or \VisIt~\cite{VisIt}.
311    
312  \subsection{Plotting Using \MATPLOTLIB}  \subsection{Plotting Using \MATPLOTLIB}
313  The \MATPLOTLIB module provides a simple and easy to use way to visualize PDE solutions (or other \Data objects).  The \MATPLOTLIB module provides a simple and easy-to-use way to visualize PDE solutions (or other \Data objects).
314  To hand over data from \escript to \MATPLOTLIB the values need to mapped onto a rectangular grid  To hand over data from \escript to \MATPLOTLIB the values need to mapped onto
315  \footnote{Users of Debian 5(Lenny) please note: this example makes use of the \function{griddata} method in \module{matplotlib.mlab}.  a rectangular grid\footnote{Users of Debian 5 (Lenny) please note: this example
316    makes use of the \function{griddata} method in \module{matplotlib.mlab}.
317  This method is not part of version 0.98.1 which is available with Lenny.  This method is not part of version 0.98.1 which is available with Lenny.
318  If you wish to use contour plots, you may need to install a later version.  If you wish to use contour plots, you may need to install a later version.
319  Users of Ubuntu 8.10 or later should be fine.}. We will make use  Users of Ubuntu 8.10 or later should be fine.}. We will make use of the \numpy module.
 of the \numpy module.  
320    
321  First we need to create a rectangular grid. We use the following statements:  First we need to create a rectangular grid which is accomplished by the following statements:
322  \begin{python}  \begin{python}
323  import numpy    import numpy
324  x_grid = numpy.linspace(0.,1.,50)    x_grid = numpy.linspace(0., 1., 50)
325  y_grid = numpy.linspace(0.,1.,50)    y_grid = numpy.linspace(0., 1., 50)
326  \end{python}  \end{python}
327  \var{x_grid} is an array defining the x coordinates of the grids while  \var{x_grid} is an array defining the x coordinates of the grid while
328  \var{y_grid} defines the y coordinates of the grid. In this case we use $50$ points over the interval $[0,1]$  \var{y_grid} defines the y coordinates of the grid.
329  in both directions.  In this case we use $50$ points over the interval $[0,1]$ in both directions.
330    
331  Now the values created by \escript need to be interpolated to this grid. We will use the \MATPLOTLIB  Now the values created by \escript need to be interpolated to this grid.
332  \function{mlab.griddata} function to do this. We can easily extract spatial coordinates as a \var{list} by  We will use the \MATPLOTLIB \function{mlab.griddata} function to do this.
333    Spatial coordinates are easily extracted as a \var{list} by
334  \begin{python}  \begin{python}
335  x=mydomain.getX()[0].toListOfTuples()    x=mydomain.getX()[0].toListOfTuples()
336  y=mydomain.getX()[1].toListOfTuples()    y=mydomain.getX()[1].toListOfTuples()
337  \end{python}  \end{python}
338  In principle we can apply the same \member{toListOfTuples} method to extract the values from the  In principle we can apply the same \member{toListOfTuples} method to extract the values from the PDE solution \var{u}.
339  PDE solution \var{u}. However, we have to make sure that the \Data object we extract the values from  However, we have to make sure that the \Data object we extract the values from
340  uses the same \FunctionSpace as we have us when extracting \var{x} and \var{y}. We apply the  uses the same \FunctionSpace as we have used when extracting \var{x} and \var{y}.
341  \function{interpolation} to \var{u} before extraction to achieve this:  We apply the \function{interpolation} to \var{u} before extraction to achieve this:
342  \begin{python}  \begin{python}
343  z=interpolate(u,mydomain.getX().getFunctionSpace())    z=interpolate(u, mydomain.getX().getFunctionSpace())
344  \end{python}  \end{python}
345  The values in \var{z} are now the values at the points with the coordinates given by \var{x} and \var{y}. These  The values in \var{z} are the values at the points with the coordinates given by \var{x} and \var{y}.
346  values are now interpolated to the grid defined by \var{x_grid} and \var{y_grid} by using  These values are interpolated to the grid defined by \var{x_grid} and \var{y_grid} by using
347  \begin{python}  \begin{python}
348  import matplotlib    import matplotlib
349  z_grid = matplotlib.mlab.griddata(x,y,z,xi=x_grid,yi=y_grid )    z_grid = matplotlib.mlab.griddata(x, y, z, xi=x_grid, yi=y_grid)
350  \end{python}  \end{python}
351  \var{z_grid} gives now the values of the PDE solution \var{u} at the grid. The values can be plotted now  Now \var{z_grid} gives the values of the PDE solution \var{u} at the grid which can be plotted using \function{contourf}:
 using the \function{contourf}:  
352  \begin{python}  \begin{python}
353  matplotlib.pyplot.contourf(x_grid, y_grid, z_grid, 5)    matplotlib.pyplot.contourf(x_grid, y_grid, z_grid, 5)
354  matplotlib.pyplot.savefig("u.png")    matplotlib.pyplot.savefig("u.png")
355  \end{python}  \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  Here we use 5 contours. The last statement writes the plot to the file \file{u.png} in the PNG format.
357    Alternatively, one can use
358  \begin{python}  \begin{python}
359  matplotlib.pyplot.contourf(x_grid, y_grid, z_grid, 5)    matplotlib.pyplot.contourf(x_grid, y_grid, z_grid, 5)
360  matplotlib.pyplot.show()    matplotlib.pyplot.show()
361  \end{python}  \end{python}
362  which gives an interactive browser window.  which gives an interactive browser window.
363    
364  \begin{figure}  \begin{figure}
365  \centerline{\includegraphics[width=\figwidth]{figures/FirstStepResultMATPLOTLIB}}  \centerline{\includegraphics[width=\figwidth]{figures/FirstStepResultMATPLOTLIB}}
366  \caption{Visualization of the Poisson Equation Solution for $f=1$ using \MATPLOTLIB.}  \caption{Visualization of the Poisson Equation Solution for $f=1$ using \MATPLOTLIB}
367  \label{fig:FirstSteps.3b}  \label{fig:FirstSteps.3b}
368  \end{figure}  \end{figure}
369    
370  Now we can write the script to solve our Poisson problem  Now we can write the script to solve our Poisson problem
371  \begin{python}  \begin{python}
372  from esys.escript import *    from esys.escript import *
373  from esys.escript.linearPDEs import Poisson    from esys.escript.linearPDEs import Poisson
374  from esys.finley import Rectangle    from esys.finley import Rectangle
375  import numpy    import numpy
376  import matplotlib    import matplotlib
377  import pylab    import pylab
378  # generate domain:    # generate domain:
379  mydomain = Rectangle(l0=1.,l1=1.,n0=40, n1=20)    mydomain = Rectangle(l0=1.,l1=1.,n0=40, n1=20)
380  # define characteristic function of Gamma^D    # define characteristic function of Gamma^D
381  x = mydomain.getX()    x = mydomain.getX()
382  gammaD = whereZero(x[0])+whereZero(x[1])    gammaD = whereZero(x[0])+whereZero(x[1])
383  # define PDE and get its solution u    # define PDE and get its solution u
384  mypde = Poisson(domain=mydomain)    mypde = Poisson(domain=mydomain)
385  mypde.setValue(f=1,q=gammaD)    mypde.setValue(f=1,q=gammaD)
386  u = mypde.getSolution()    u = mypde.getSolution()
387  # interpolate u to a matplotlib grid:    # interpolate u to a matplotlib grid:
388  x_grid = numpy.linspace(0.,1.,50)    x_grid = numpy.linspace(0.,1.,50)
389  y_grid = numpy.linspace(0.,1.,50)    y_grid = numpy.linspace(0.,1.,50)
390  x=mydomain.getX()[0].toListOfTuples()    x=mydomain.getX()[0].toListOfTuples()
391  y=mydomain.getX()[1].toListOfTuples()    y=mydomain.getX()[1].toListOfTuples()
392  z=interpolate(u,mydomain.getX().getFunctionSpace())    z=interpolate(u,mydomain.getX().getFunctionSpace())
393  z_grid = matplotlib.mlab.griddata(x,y,z,xi=x_grid,yi=y_grid )    z_grid = matplotlib.mlab.griddata(x,y,z,xi=x_grid,yi=y_grid )
394  # interpolate u to a rectangular grid:    # interpolate u to a rectangular grid:
395  matplotlib.pyplot.contourf(x_grid, y_grid, z_grid, 5)    matplotlib.pyplot.contourf(x_grid, y_grid, z_grid, 5)
396  matplotlib.pyplot.savefig("u.png")    matplotlib.pyplot.savefig("u.png")
397  \end{python}  \end{python}
398  The entire code is available as \file{poisson\hackscore matplotlib.py} in the \ExampleDirectory.  The entire code is available as \file{poisson\hackscore matplotlib.py} in the \ExampleDirectory.
399  You can run the script using the {\it escript} environment  You can run the script using the {\it escript} environment
400  \begin{verbatim}  \begin{verbatim}
401    run-escript poisson_matplotlib.py  run-escript poisson_matplotlib.py
402  \end{verbatim}  \end{verbatim}
403  This will create the \file{u.png}, see Figure~\fig{fig:FirstSteps.3b}.  This will create the \file{u.png}, see \fig{fig:FirstSteps.3b}.
404  For details on the usage of the \MATPLOTLIB module we refer to the documentation~\cite{matplotlib}.  For details on the usage of the \MATPLOTLIB module we refer to the documentation~\cite{matplotlib}.
405    
406  As pointed out, \MATPLOTLIB is restricted to the two-dimensional case and  As pointed out, \MATPLOTLIB is restricted to the two-dimensional case and
407  should be used for small problems only. It can not be used under \MPI as the \member{toListOfTuples} method is  should be used for small problems only.
408  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.}.  It can not be used under \MPI as the \member{toListOfTuples} method is not
409    safe under \MPI\footnote{The phrase 'safe under \MPI' means that a program
410    will produce correct results when run on more than one processor under \MPI.}.
411    
412  \begin{figure}  \begin{figure}
413  \centerline{\includegraphics[width=\figwidth]{figures/FirstStepResult}}  \centerline{\includegraphics[width=\figwidth]{figures/FirstStepResult}}
# Line 411  not safe under \MPI\footnote{The phrase Line 415  not safe under \MPI\footnote{The phrase
415  \label{fig:FirstSteps.3}  \label{fig:FirstSteps.3}
416  \end{figure}  \end{figure}
417    
418  \subsection{Visualization using \VTK}  \subsection{Visualization using export files}
419    
420  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.  As an alternative to \MATPLOTLIB, {\it escript} supports exporting data to
421    \VTK and \SILO files which can be read by visualization tools such as
422    mayavi\cite{mayavi} and \VisIt~\cite{VisIt}. This method is \MPI safe and
423    works with large 2D and 3D problems.
424    
425  To write the solution \var{u} in  Poisson problem to the file \file{u.xml} one need to add the line  To write the solution \var{u} of the Poisson problem in the \VTK file format
426    to the file \file{u.vtu} one needs to add:
427  \begin{python}  \begin{python}
428  saveVTK("u.xml",sol=u)    saveVTK("u.vtu", sol=u)
429  \end{python}  \end{python}
430  The solution \var{u} is now available in the \file{u.xml} tagged with the name "sol".  This file can then be opened in a \VTK compatible visualization tool where the
431    solution is accessible by the name {\it sol}.
432    
433  The Poisson problem script is now  The Poisson problem script is now
434  \begin{python}  \begin{python}
# Line 436  The Poisson problem script is now Line 445  The Poisson problem script is now
445    mypde.setValue(f=1,q=gammaD)    mypde.setValue(f=1,q=gammaD)
446    u = mypde.getSolution()    u = mypde.getSolution()
447    # write u to an external file    # write u to an external file
448    saveVTK("u.xml",sol=u)    saveVTK("u.vtu",sol=u)
449  \end{python}  \end{python}
450  The entire code is available as \file{poisson\hackscore VTK.py} in the \ExampleDirectory.  The entire code is available as \file{poisson\hackscore vtk.py} in the \ExampleDirectory.
451    
452  You can run the script using the {\it escript} environment  You can run the script using the {\it escript} environment and visualize the
453  and visualize the solution using \mayavi:  solution using \mayavi:
454  \begin{verbatim}  \begin{verbatim}
455    run-escript poisson_VTK.py  run-escript poisson_VTK.py
456    mayavi2 -d u.xml -m SurfaceMap  mayavi2 -d u.vtu -m SurfaceMap
457  \end{verbatim}  \end{verbatim}
458  The result is shown in Figure~\fig{fig:FirstSteps.3}.  The result is shown in \fig{fig:FirstSteps.3}.
459    

Legend:
Removed from v.3277  
changed lines
  Added in v.3278

  ViewVC Help
Powered by ViewVC 1.1.26