Merged in prabhuramachandran/pysph/openmp (pull request #161)

OpenMP support for PySPH

docs/source/installation.rst

@@ -6,8 +6,9 @@ Installation and getting started
 
 To install PySPH, you need a working Python environment with the required
 dependencies installed. You may use any of the available Python distributions.
-If you are new to Python we recommend `Enthought Canopy`_. PySPH will also
-work fine with Anaconda_ or other environments like WinPython_.  The following
+Note that PySPH is currently tested with Python-2.x and not with 3.x (yet).  If
+you are new to Python we recommend `Enthought Canopy`_. PySPH will work
+fine with Anaconda_ or other environments like WinPython_.  The following
 instructions should help you get started.
 
 Since there is a lot of information here, we suggest that you skim the section
@@ -43,7 +44,7 @@ manager, or using pip_ or ``easy_install``.  For more detailed instructions on
 how to do this for different distributions, see below.
 
 Running PySPH requires a working C/C++ compiler on your machine.  On Linux/OS X
-the gcc toolchain will work well on Windows, you will need to have `Microsoft
+the gcc toolchain will work well.  On Windows, you will need to have `Microsoft
 Visual C++ Compiler for Python 2.7
 <http://www.microsoft.com/en-us/download/details.aspx?id=44266>`_ or an
 equivalent compiler.  More details are available below.
@@ -68,6 +69,9 @@ Optional dependencies
 
 The optional dependencies are:
 
+ - OpenMP_: PySPH can use OpenMP if it is available.  Installation instructions
+   are available below.
+
  - Mayavi_: PySPH provides a convenient viewer to visualize the output of
    simulations. This viewer is called ``pysph_viewer`` and requires Mayavi_ to
    be installed.  Since this is only a viewer it is optional for use, however,
@@ -85,7 +89,7 @@ Zoltan_ is very unlikely to be already packaged and will need to be compiled.
 .. _Mayavi: http://code.enthought.com/projects/mayavi
 .. _mpi4py: http://mpi4py.scipy.org/
 .. _Zoltan: http://www.cs.sandia.gov/zoltan/
-
+.. _OpenMP: http://openmp.org/
 
 Building and linking PyZoltan on OSX/Linux
 -------------------------------------------
@@ -131,6 +135,10 @@ install the dependencies using::
     $ sudo apt-get install build-essential python-dev python-numpy \
         python-mako cython python-nose mayavi2 python-qt4 python-virtualenv
 
+OpenMP_ is typically available but can be installed with::
+
+    $ sudo apt-get install libgomp1
+
 If you need parallel support::
 
     $ sudo apt-get install libopenmpi-dev python-mpi4py
@@ -176,7 +184,30 @@ On OS X, your best bet is to install `Enthought Canopy`_ or Anaconda_ or some
 other Python distribution.  Ensure that you have gcc or clang installed by
 installing XCode.  See `this
 <http://stackoverflow.com/questions/12228382/after-install-xcode-where-is-clang>`_
-if you installed XCode but can't find clang or gcc.
+if you installed XCode but can't find clang or gcc.  
+
+^^^^^^^^^^^^^
+OpenMP on OSX
+^^^^^^^^^^^^^
+
+If you need to use OpenMP_, the default clang compiler on OSX does not support
+it.  There are some experimental versions available.  One easy to install
+option is to use brew to install gcc.  For example you can try::
+
+    $ sudo brew install gcc
+
+The build may not see ``omp.h`` and you can work around this by manually
+linking to it like so (modify this to suit your installation)::
+
+    $ cd /usr/local/include
+    $ sudo ln -s ../Cellar/gcc/4.9.2_1/lib/gcc/4.9/gcc/x86_64-apple-darwin12.6.0/4.9.2/include/omp.h .
+
+Once this is done, you need to use this as your default compiler, you can tell
+the Python to use this by setting::
+
+    $ export CC=gcc-4.9
+    $ export CXX=g++-4.9
+
 
 
 ^^^^^^^^^^^^^
@@ -345,7 +376,13 @@ build PySPH with `Microsoft's Visual C++ for Python 2.7
 <http://www.microsoft.com/en-us/download/details.aspx?id=44266>`_.  We
 recommend that you download and install the ``VCForPython27.msi`` available
 from the `link
-<http://www.microsoft.com/en-us/download/details.aspx?id=44266>`_.
+<http://www.microsoft.com/en-us/download/details.aspx?id=44266>`_.  **Make sure
+you install the system requirements specified on that page**.  For example, on
+Windows 7 you will need to install the Microsoft Visual C++ 2008 SP1
+Redistributable Package for your platform (x86 for 32 bit or x64 for 64 bit)
+and on Windows 8 you will need to install the .NET framework 3.5.  Please look
+at the link given above, it should be fairly straightforward.  Note that doing
+this will also get OpenMP_ working for you.
 
 After you do this, you will find a "Microsoft Visual C++ Compiler Package for
 Python 2.7" in your Start menu.  Choose a suitable command prompt from this
@@ -457,7 +494,7 @@ directory::
     With a virtualenv, one should be careful while running things like
     ``ipython`` or ``nosetests`` as these are sometimes also installed on the
     system in ``/usr/bin``.  If you suspect that you are not running the
-    correct Python, you could simply run (on *nix/OS X)::
+    correct Python, you could simply run (on Linux/OS X)::
 
         $ python `which ipython`
 
@@ -539,10 +576,21 @@ You should be all set now and should next consider :ref:`running-the-tests`.
 Running the tests
 ------------------
 
-To test PySPH you can run::
+If you installed PySPH using ``python setup.py develop`` you can run the tests as::
 
    $ python -m nose.core pysph
 
+If you installed PySPH using ``python setup.py install`` you should run the
+tests as::
+
+   $ python -m nose.core -w docs pysph
+
+This is because nosetests will incorrectly  will pick up the local pysph
+packages instead of the installed version.  The ``-w`` option just changes the
+active directory to ``docs``.  Alternatively, change directory to some other
+directory that does not contain the directory ``pysph`` in it and run the first
+command above.
+
 If you see errors you might want more verbose reporting which you can get
 with::
 
@@ -619,6 +667,43 @@ Which runs the problem of the collision of two elastic rings:
 The auto-generated code for the example resides in the directory
 ``~/.pysph/source``. A note of caution however, it's not for the faint hearted.
 
+^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
+Running the examples with OpenMP
+^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
+
+If you have OpenMP available run any of the examples as follows::
+
+    $ python elliptical_drop.py --openmp
+
+This should run faster if you have multiple cores on your machine.  If
+you wish to change the number of threads to run simultaneously, you can
+try the following::
+
+    $ OMP_NUM_THREADS=8 python elliptical_drop.py --openmp
+
+You may need to set the number of threads to about 4 times the number of
+physical cores on your machine to obtain the most scale-up.  If you wish
+to time the actual scale up of the code with and without OpenMP you may
+want to disable any output (which will be serial), you can do this
+like::
+
+    $ python elliptical_drop.py --disable-output --openmp
+
+
+^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
+Running the examples with MPI
+^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
+
+If you compiled PySPH with Zoltan_ and have mpi4py_ installed you may run any
+of the examples with MPI as follows::
+
+    $ mpirun -np 4 python dam_break3D.py
+
+This may not give you significant speedup if the problem is too small.  You can
+also combine OpenMP and MPI if you wish.  You should take care to setup the MPI
+host information suitably to utilize the processors effectively.
+
+
 --------------------------------------
 Organization of the ``pysph`` package
 --------------------------------------

docs/source/using_pysph.rst

@@ -140,7 +140,7 @@ example if the total mass or total force on a particle array needs to be
 calculated, a fixed size constant can be added.  This can be done by adding a
 ``constant`` to the array as illustrated below:
 
-.. code-block::python
+.. code-block:: python
 
     >>> pa.add_constant('total_mass', 0.0)
     >>> pa.add_constant('total_force', [0.0, 0.0, 0.0])
@@ -156,7 +156,7 @@ The constants can also set in the constructor of the :py:class:`ParticleArray`
 by passing a dictionary of constants as a ``constants`` keyword argument.  For
 example:
 
-.. code-block::python
+.. code-block:: python
 
     >>> pa = ParticleArray(
     ...     name='test', x=x,
@@ -192,6 +192,13 @@ cells have a length of :math:`\text{radius_scale} \times h_{\text{max}}`,
 where :math:`h_{\text{max}}` is the maximum smoothing length of *all*
 particles assigned to the local processor.
 
+Note that the ``NNPS`` classes also support caching the neighbors
+computed.  This is useful if one needs to reuse the same set of
+neighbors.  To enable this, simply pass ``cache=True`` to the
+constructor::
+
+    >>> nps = nnps.LinkedListNNPS(dim=3, particles=particles, cache=True)
+
 Since we allow a list of particle arrays, we need to distinguish
 between *source* and *destination* particle arrays in the neighbor
 queries.
@@ -214,7 +221,19 @@ With these definitions, we can query for nearest neighbors like so:
 where `src_index`, `dst_index` and `d_idx` are integers. This will
 return, for the *d_idx* particle of the *dst_index* particle array
 (species), nearest neighbors from the *src_index* particle array
-(species).
+(species).  Passing the `src_index` and `dst_index` every time is
+repetitive so an alternative API is to call ``set_context`` as done
+below::
+
+    >>> nps.set_context(src_index=0, dst_index=0)
+
+If the ``NNPS`` instance is configured to use caching, then it will also
+pre-compute the neighbors very efficiently.  Once the context is set one
+can get the neighbors as::
+
+    >>> nps.get_nearest_neighbors(d_idx, nbrs)
+
+Where `d_idx` and `nbrs` are as discussed above.
 
 If we want to re-compute the data structure for a new distribution of
 particles, we can call the :py:meth:`NNPS.update` method:
@@ -230,24 +249,25 @@ Periodic domains
 ^^^^^^^^^^^^^^^^^^^^^^
 
 The constructor for the :py:class:`NNPS` accepts an optional argument
-(:py:class:`DomainLimits`) that is used to delimit the maximum
+(:py:class:`DomainManager`) that is used to delimit the maximum
 spatial extent of the simulation domain. Additionally, this argument
 is also used to indicate the extents for a periodic domain. We
-construct a :py:class:`DomainLimits` object like so
+construct a :py:class:`DomainManager` object like so
 
 .. code-block:: python
 
-   >>> from pysph.base.nnps import DomainLimits
+   >>> from pysph.base.nnps import DomainManager
    >>> from pysph.base.point import Point
-   >>> domain = DomainLimits(xmin, xmax, ymin, ymax, zmin, zmax,
-                             periodic_in_x, periodic_in_y, periodic_in_z)
+   >>> domain = DomainManager(xmin, xmax, ymin, ymax, zmin, zmax,
+                              periodic_in_x=True, periodic_in_y=True,
+                              periodic_in_z=False)
 
 where `xmin ... zmax` are floating point arguments delimiting the
 simulation domain and `periodic_in_x,y,z` are bools defining the
 periodic axes.
 
 When the :py:class:`NNPS` object is constructed with this
-:py:class:`DomainLimits`, care is taken to create periodic ghosts for
+:py:class:`DomainManager`, care is taken to create periodic ghosts for
 particles in the vicinity of the periodic boundaries. These *ghost*
 particles are given a special **tag** defined by
 :py:class:`ParticleTAGS`
@@ -407,7 +427,7 @@ distribution is given below
    from pyzoltan.core.carray import UIntArray
    from pysph.base.utils import utils
    from pysph.base.kernels import CubicSpline
-   from pysph.base.nnps import DomainLimits, LinkedListNNPS
+   from pysph.base.nnps import DomainManager, LinkedListNNPS
 
    # NumPy
    import numpy
@@ -423,8 +443,8 @@ distribution is given below
    # Create the particle array
    pa = utils.get_particle_array(x=x,y=y,h=h,m=m)
 
-   # Create the periodic DomainLimits object and NNPS
-   domain = DomainLimits(xmin=0., xmax=1., ymin=0., ymax=1., periodic_in_x=True, periodic_in_y=True)
+   # Create the periodic DomainManager object and NNPS
+   domain = DomainManager(xmin=0., xmax=1., ymin=0., ymax=1., periodic_in_x=True, periodic_in_y=True)
    nps = LinkedListNNPS(dim=2, particles=[pa,], radius_scale=2.0, domain=domain)
 
    # The SPH kernel. The dimension argument is needed for the correct normalization constant
@@ -536,4 +556,4 @@ eaxmples, check out the tutorials: :ref:`tutorials`.
 .. _RIB: http://www.cs.sandia.gov/Zoltan/ug_html/ug_alg_rib.html
 .. _HSFC: http://www.cs.sandia.gov/Zoltan/ug_html/ug_alg_hsfc.html
 
-..  LocalWords:  DomainLimits maximum
+..  LocalWords:  DomainManager maximum

pysph/base/config.py

@@ -0,0 +1,38 @@
+"""Simple configuration options for PySPH.
+
+Do not import any PySPH specific extensions here, if you must, do the import
+inside the function/method.
+"""
+
+import sys
+
+
+class Config(object):
+    def __init__(self):
+        self._use_openmp = None
+
+    @property
+    def use_openmp(self):
+        if self._use_openmp is None:
+            self._use_openmp = self._use_openmp_default()
+        return self._use_openmp
+
+    @use_openmp.setter
+    def use_openmp(self, value):
+        self._use_openmp = value
+
+    def _use_openmp_default(self):
+        return False
+
+
+_config = None
+
+def get_config():
+    global _config
+    if _config is None:
+        _config = Config()
+    return _config
+
+def set_config(config):
+    global _config
+    _config = config

pysph/base/cython_generator.py

@@ -11,7 +11,9 @@
 from textwrap import dedent
 import types
 
-from ast_utils import get_assigned, has_return
+from pysph.base.ast_utils import get_assigned, has_return
+from pysph.base.config import get_config
+
 
 logger = logging.getLogger(__name__)
 
@@ -97,6 +99,7 @@ def __init__(self, known_types=None, python_methods=False):
         # Methods to not wrap.
         self.ignore_methods = ['_cython_code_']
         self.known_types = known_types if known_types is not None else {}
+        self._config = get_config()
 
     ##### Public protocol #####################################################
 
@@ -201,8 +204,13 @@ def _get_c_method_spec(self, name, returns, args):
 
         c_ret = 'double' if returns else 'void'
         c_arg_def = ', '.join(c_args)
-        cdefn = 'cdef inline {ret} {name}({arg_def}):'\
-                     .format(ret=c_ret, name=name, arg_def=c_arg_def)
+        if self._config.use_openmp:
+            gil = " nogil" if name != "reduce" else ""
+        else:
+            gil = ""
+        cdefn = 'cdef inline {ret} {name}({arg_def}){gil}:'.format(
+            ret=c_ret, name=name, arg_def=c_arg_def, gil=gil
+        )
 
         return cdefn