Welcome to the PySPH documentation!¶

Features¶

• User scripts and equations are written in pure Python.
• Flexibility to define arbitrary SPH equations operating on particles.
• Ability to define your own multi-step integrators in pure Python.
• High-performance: our performance is comparable to hand-written solvers implemented in FORTRAN.
• Seamless multi-core support with OpenMP.
• Seamless parallel integration using Zoltan.
• BSD license.

SPH formulations¶

Currently, PySPH has numerous examples to solve the viscous, incompressible Navier-Stokes equations using the weakly compressible (WCSPH) approach. The following formulations are currently implemented:

• Weakly Compressible SPH (WCSPH) for free-surface flows (Gesteira et al. 2010, Journal of Hydraulic Research, 48, pp. 6–27)

3D dam-break past an obstacle SPHERIC benchmark Test 2

• Transport Velocity Formulation for incompressilbe fluids (Adami et al. 2013, JCP, 241, pp. 292–307).

Streamlines for a driven cavity

• SPH for elastic dynamics (Gray et al. 2001, CMAME, Vol. 190, pp 6641–6662)

Collision of two elastic rings.

• Compressible SPH (Puri et al. 2014, JCP, Vol. 256, pp 308–333)

Credits¶

PySPH is primarily developed at the Department of Aerospace Engineering, IIT Bombay. We are grateful to IIT Bombay for the support. Our primary goal is to build a powerful SPH-based tool for both application and research. We hope that this makes it easy to perform reproducible computational research.

Lead developers:

• Prabhu Ramachandran
• Kunal Puri

Earlier developers:

• Pankaj Pandey (stress solver and improved load balancing, 2011)
• Chandrashekhar Kaushik (original parallel and serial implementation in 2009)

The following have contributed bug-fixes, features, documentation etc.

• Arkopal Dutt, Arpit Agarwal, Sarang Minhas, S Saravanan, Vishnu Sivadasan

Citing PySPH¶

We are in the process of writing up an article on the new PySPH framework as it stands today. In the meanwhile, if you use PySPH and wish to cite it you may use this:

• Prabhu Ramachandran and Kunal Puri, PySPH: A framework for parallel particle simulations, In proceedings of the 3rd International Conference on Particle-Based Methods (Particles 2013), Stuttgart, Germany, 18th September 2013.

History¶

• 2009: PySPH started with a simple Cython based 1D implementation written by Prabhu.
• 2009-2010: Chandrashekhar Kaushik worked on a full 3D SPH implementation with a more general purpose design. The implementation was in a mix of Cython and Python.
• 2010-2012: The previous implementation was a little too complex and was largely overhauled by Kunal and Pankaj. This became the PySPH 0.9beta release. The difficulty with this version was that it was almost entirely written in Cython, making it hard to extend or add new formulations without writing more Cython code. Doing this was difficult and not too pleasant. In addition it was not as fast as we would have liked it. It ended up feeling like we might as well have implemented it all in C++ and exposed a Python interface to that.
• 2011-2012: Kunal also implemented SPH2D and another internal version called ZSPH in Cython which included Zoltan based parallelization using PyZoltan. This was specific to his PhD research and again required writing Cython making it difficult for the average user to extend.
• 2013-present In early 2013, Prabhu reimplemented the core of PySPH to be almost entirely auto-generated from pure Python. The resulting code was faster than previous implementations and very easy to extend entirely from pure Python. Kunal and Prabhu integrated PyZoltan into PySPH and the current version of PySPH was born. Subsequently, OpenMP support was also added in 2015.

Support¶

If you have any questions or are running into any difficulties with PySPH, please email or post your questions on the pysph-users mailing list here: https://groups.google.com/d/forum/pysph-users

Please also take a look at the PySPH issue tracker.

Dependencies¶

$ ./build_zoltan.sh INSTALL_PREFIX

$ export ZOLTAN=$INSTALL_PREFIX

$ export ZOLTAN_INCLUDE=$INSTALL_PREFIX/include
$ export ZOLTAN_LIBRARY=$INSTALL_PREFIX/lib

Installing the dependencies on GNU/Linux¶

GNU/Linux is probably the easiest platform to install PySPH. On Ubuntu one may 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
$ ./build_zoltan.sh ~/zoltan # Replace ~/zoltan with what you want
$ export ZOLTAN=~/zoltan

On Linux it is probably best to install PySPH into its own virtual environment. This will allow you to install PySPH as a user without any superuser priviledges. See the section below on Using a virtualenv for PySPH. In short do the following:

$ virtualenv --system-site-packages pysph_env
$ source pysph_env/bin/activate
$ pip install cython --upgrade # if you have an old version.

You should be set now and should skip to Downloading PySPH and Building and Installing PySPH.

If you are using Enthought Canopy or Anaconda the instructions in the section Installing the dependencies on Mac OS X will be useful as the instructions are similar.

Installing the dependencies on Mac OS X¶

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 if you installed XCode but can’t find clang or gcc.

$ sudo brew install gcc

$ 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 .

$ export CC=gcc-4.9
$ export CXX=g++-4.9

$ pip install cython mako

$ enpkg mayavi

$ conda install cython mayavi
$ pip install mako

$ sudo brew install open-mpi

$ cd /tmp
$ tar xvzf ~/Downloads/mpi4py-1.3.1.tar.gz
$ cd mpi4py-1.3.1
$ export MACOSX_DEPLOYMENT_TARGET=10.7
$ export SDKROOT=/Applications/Xcode.app/Contents/Developer/Platforms/MacOSX.platform/Developer/SDKs/MacOSX10.7.sdk/
$ python setup.py install

$ cd pysph
$ ./build_zoltan.sh ~/zoltan # Replace ~/zoltan with what you want
$ export ZOLTAN=~/zoltan

Installing the dependencies on Windows¶

While it should be possible to use mpi4py and Zoltan on Windows, we do not at this point have much experience with this. Feel free to experiment and let us know if you’d like to share your instructions. The following instructions are all without parallel support.

> pip install cython mako

> enpkg mayavi

$ conda install cython mayavi
$ pip install mako

> cd pysph
> .\windows_env.bat

> SET DISTUTILS_USE_SDK=1
> SET MSSdk=1

Using a virtualenv for PySPH¶

A virtualenv allows you to create an isolated environment for PySPH and its related packages. This is useful in a variety of situations.

• Your OS does not provide a recent enough Cython version (say you are running Debian stable).
• You do not have root access to install any packages PySPH requires.
• You do not want to mess up your system files and wish to localize any installations inside directories you control.
• You wish to use other packages with conflicting requirements.
• You want PySPH and its related packages to be in an “isolated” environment.

You can either install virtualenv (or ask your system administrator to) or just download the virtualenv.py script and use it (run python virtualenv.py after you download the script).

Create a virtualenv like so:

$ virtualenv --system-site-packages pysph_env

This creates a directory called pysph_env which contains all the relevant files for your virtualenv, this includes any new packages you wish to install into it. You can delete this directory if you don’t want it anymore for some reason. This virtualenv will also “inherit” packages from your system. Hence if your system administrator already installed NumPy it may be imported from your virtual environment and you do not need to install it. This is very useful for large packages like Mayavi, Qt etc.

Once you create a virtualenv you can activate it as follows (on a bash shell):

$ source pysph_env/bin/activate

On Windows you run a bat file as follows:

$ pysph_env/bin/activate

This sets up the PATH to point to your virtualenv’s Python. You may now run any normal Python commands and it will use your virtualenv’s Python. For example you can do the following:

$ virtualenv myenv
$ source myenv/bin/activate
(myenv) $ pip install Cython mako nose
(myenv) $ cd pysph
(myenv) $ python setup.py install

Now PySPH will be installed into myenv. You may deactivate your virtualenv using the deactivate command:

(myenv) $ deactivate
$

On Windows, use myenv\Scripts\activate.bat and myenv\Scripts\deactivate.bat.

If for whatever reason you wish to delete myenv just remove the entire directory:

$ rm -rf myenv

$ python `which ipython`

$ venv --help
$ venv --system-site-packages myenv
$ source myenv/bin/activate

Downloading PySPH¶

One way to install PySPH is to use pip

$ pip install PySPH

This will install PySPH, and you should be able to import it and use the modules with your Python scripts that use PySPH. This will not give you access to the examples. The best way to get those is to get PySPH from git or download a tarball or ZIP as described below.

To get PySPH using git type the following

$ git clone https://bitbucket.org/pysph/pysph.git

If you do not have git or do not wish to bother with it, you can get a ZIP or tarball from the pysph site. You can unzip/untar this and use the sources.

In the instructions, we assume that you have the pysph sources in the directory pysph and are inside the root of this directory. For example:

$ unzip pysph-pysph-*.zip
$ cd pysph-pysph-1ce*

or if you cloned the repository:

$ git clone https://bitbucket.org/pysph/pysph.git
$ cd pysph

Once you have downloaded PySPH you should be ready to build and install it, see Building and Installing PySPH.

Building and Installing PySPH¶

Once you have the dependencies installed you can install PySPH with:

$ pip install PySPH

If you downloaded PySPH using git or used a tarball you can do:

$ python setup.py install

You could also do:

$ python setup.py develop

This is useful if you are tracking the latest version of PySPH via git. With git you can update the sources and rebuild using:

$ git pull
$ python setup.py develop

You should be all set now and should next consider Running the tests.

Running the tests¶

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:

$ python -m nose.core -v pysph

This should run all the tests that do not take a long while to complete. If this fails, please contact the pysph-users mailing list or send us email.

Once you run the tests, you should see the section on Running the examples.

Running the examples¶

You can verify the installation by exploring some examples:

$ cd examples
$ python elliptical_drop.py

Try this:

$ python elliptical_drop.py -h

to see the different options supported by each example. You can view the data generated by the simulation (after the simulation is complete or during the simulation) by running the pysph_viewer application. To view the simulated data you may do:

$ pysph_viewer elliptical_drop_output/*.npz

If you have Mayavi installed this should show a UI that looks like:

There are other examples like those in the transport_velocity directory:

$ cd transport_velocity
$ python cavity.py

This runs the driven cavity problem using the transport velocity formulation of Adami et al. You can verify the results for this problem using the helper script examples/transport_velocity/ldcavity_results.py to plot, for example the streamlines:

If you want to use PySPH for elastic dynamics, you can try some of the examples from Gray et al., Comput. Methods Appl. Mech. Engrg. 190 (2001), 6641-6662:

$ cd examples/solid_mech
$ python rings.py

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.

$ python elliptical_drop.py --openmp

$ OMP_NUM_THREADS=8 python elliptical_drop.py --openmp

$ python elliptical_drop.py --disable-output --openmp

$ mpirun -np 4 python dam_break3D.py

Organization of the pysph package¶

PySPH is organized into several sub-packages. These are:

• pysph.base: This subpackage defines the pysph.base.particle_array.ParticleArray, pysph.base.carray.CArray (which are used by the particle arrays), the various SPH Kernels, the nearest neighbor particle search (NNPS) code, and the Cython code generation utilities.
• pysph.sph: Contains the various SPH equations, the Integrator module and associated integration steppers, and the code generation for the SPH looping. pysph.sph.wc contains the equations for the weakly compressible formulation. pysph.sph.solid_mech contains the equations for solid mechanics and pysph.sph.misc has miscellaneous equations.
• pysph.solver: Provides the pysph.solver.solver.Solver, the pysph.solver.application.Application and a convenient way to interact with the solver as it is running.
• pysph.parallel: Provides the parallel functionality.
• pysph.tools: Provides some useful tools including the pysph_viewer which is based on Mayavi.

Imports¶

Taking a look at the example, the first several lines are imports of various modules:

numpy import ones_like, mgrid, sqrt, array, savez
from time import time

# PySPH base and carray imports
from pysph.base.utils import get_particle_array_wcsph
from pysph.base.kernels import CubicSpline
from pyzoltan.core.carray import LongArray

# PySPH solver and integrator
from pysph.solver.application import Application
from pysph.solver.solver import Solver
from pysph.sph.integrator import PECIntegrator

# PySPH sph imports
from pysph.sph.basic_equations import ContinuityEquation, XSPHCorrection
from pysph.sph.wc.basic import TaitEOS, MomentumEquation

Functions for loading/generating the particles¶

Next in the code are two functions called exact_solution and get_circular_patch. The former produces an exact solution for comparison, the latter looks like:

def get_circular_patch(dx=0.025, **kwargs):
    """Create the circular patch of fluid."""
    name = 'fluid'
    x,y = mgrid[-1.05:1.05+1e-4:dx, -1.05:1.05+1e-4:dx]
    x = x.ravel()
    y = y.ravel()

    m = ones_like(x)*dx*dx
    h = ones_like(x)*hdx*dx
    rho = ones_like(x) * ro

    p = ones_like(x) * 1./7.0 * co**2
    cs = ones_like(x) * co

    u = -100*x
    v = 100*y

    # remove particles outside the circle
    indices = []
    for i in range(len(x)):
        if sqrt(x[i]*x[i] + y[i]*y[i]) - 1 > 1e-10:
            indices.append(i)

    pa = get_particle_array_wcsph(x=x, y=y, m=m, rho=rho, h=h, p=p, u=u, v=v,
                                  cs=cs, name=name)

    la = LongArray(len(indices))
    la.set_data(array(indices))

    pa.remove_particles(la)

    print "Elliptical drop :: %d particles"%(pa.get_number_of_particles())

    # add requisite variables needed for this formulation
    for name in ('arho', 'au', 'av', 'aw', 'ax', 'ay', 'az', 'rho0', 'u0',
                 'v0', 'w0', 'x0', 'y0', 'z0'):
        pa.add_property(name)

    return [pa,]

and is used to initialize the particles in Python. In PySPH, we use a ParticleArray object as a container for particles of a given species. You can think of a particle species as any homogenous entity in a simulation. For example, in a two-phase air water flow, a species could be used to represent each phase. A ParticleArray can be conveniently created from the command line using NumPy arrays. For example

>>> from pysph.base.utils import get_particle_array
>>> x, y = numpy.mgrid[0:1:0.01, 0:1:0.01]
>>> x = x.ravel(); y = y.ravel()
>>> pa = sph.get_particle_array(x=x, y=y)

would create a ParticleArray, representing a uniform distribution of particles on a Cartesian lattice in 2D using the helper function get_particle_array() in the base subpackage.

The ParticleArray is highly convenient, supporting methods for insertions, deletions and concatenations. In the get_circular_patch function, we use this convenience to remove a list of particles that fall outside a circular region:

pa.remove_particles(la)

where, a list of indices is provided in the form of a LongArray which, as the name suggests, is an array of 64 bit integers.

Setting up the PySPH framework¶

As we move on, we encounter instantiations of the PySPH framework objects. These are the pysph.solver.application.Application, pysph.sph.integrator.PECIntegrator and pysph.solver.solver.Solver objects:

# Create the application.
app = Application()

kernel = CubicSpline(dim=2)

integrator = PECIntegrator(fluid=WCSPHStep())

# Create and setup a solver.
solver = Solver(kernel=kernel, dim=2, integrator=integrator)

# Setup default parameters.
solver.set_time_step(1e-5)
solver.set_final_time(0.0075)

The Application makes it easy to pass command line arguments to the solver. It is also important for the seamless parallel execution of the same example. To appreciate the role of the Application consider for a moment how might we write a parallel version of the same example. At some point, we would need some MPI imports and the particles should be created in a distributed fashion. All this (and more) is handled through the abstraction of the Application which hides all this detail from the user.

Intuitively, in an SPH simulation, the role of the PECIntegrator should be obvious. In the code, we see that we ask for the “fluid” to be stepped using a WCSPHStep object. Taking a look at the get_circular_patch function once more, we notice that the ParticleArray representing the circular patch was named as fluid. So we’re essentially asking the PySPH framework to step or integrate the properties of the ParticleArray fluid using WCSPHStep. Safe to assume that the framework takes the responsibility to call this integrator at the appropriate time during a time-step.

The Solver is the main driver for the problem. It marshals a simulation and takes the responsibility (through appropriate calls to the integrator) to update the solution to the next time step. It also handles input/output and computing global quantities (such as minimum time step) in parallel.

Specifying the interactions¶

At this stage, we have the particles (represented by the fluid ParticleArray) and the framework to integrate the solution and marshall the simulation. What remains is to define how to actually go about updating properties within a time step. That is, for each particle we must “do something”. This is where the physics for the particular problem comes in.

For SPH, this would be the pairwise interactions between particles. In PySPH, we provide a specific way to define the sequence of interactions which is a list of Equation objects (see SPH equations). For the circular patch test, the sequence of interactions is relatively straightforward:

• Compute pressure from the EOS: \(p = f(\rho)\)
• Compute the rate of change of density: \(\frac{d\rho}{dt}\)
• Compute the rate of change of velocity (accelerations): \(\frac{d\boldsymbol{v}}{dt}\)
• Compute corrections for the velocity (XSPH): \(\frac{d\boldsymbol{x}}{dt}\)

We request this in PySPH like so:

# The equations of motion.
equations = [
    # Equation of state: p = f(rho)
    TaitEOS(dest='fluid', sources=None, rho0=ro, c0=co, gamma=7.0),

    # Density rate: drho/dt
    ContinuityEquation(dest='fluid',  sources=['fluid',]),

    # Acceleration: du,v/dt
    MomentumEquation(dest='fluid', sources=['fluid'], alpha=1.0, beta=1.0),

    # XSPH velocity correction
    XSPHCorrection(dest='fluid', sources=['fluid']),

    ]

Each interaction is specified through an Equation object, which is instantiated with the general syntax:

Equation(dest='array_name', sources, **kwargs)

The dest argument specifies the target or destination ParticleArray on which this interaction is going to operate on. Similarly, the sources argument specifies a list of ParticleArrays from which the contributions are sought. For some equations like the EOS, it doesn’t make sense to define a list of sources and a None suffices. The specification basically tells PySPH that for one time step of the calculation:

• Use the Tait’s EOS to update the properties of the fluid array
• Compute \(\frac{d\rho}{dt}\) for the fluid from the fluid
• Compute accelerations for the fluid from the fluid
• Compute the XSPH corrections for the fluid, using fluid as the source

With the list of equations, our problem is completely defined. PySPH now knows what to do with the particles within a time step. More importantly, this information is enough to generate code to carry out a complete SPH simulation.

Running the example¶

In the last two lines of the example, we use the Application to run the problem:

# Setup the application and solver.  This also generates the particles.
app.setup(solver=solver, equations=equations,
          particle_factory=get_circular_patch)

app.run()

We can see that the Application.setup() method is where we tell PySPH what we want it to do. We pass in the function to create the particles, the list of equations defining the problem and the solver that will be used to marshal the problem.

Many parameters can be configured via the command line, and these will override any parameters setup before the app.setup call. For example one may do the following to find out the various options:

$ python elliptical_drop.py -h

If we run the example without any arguments it will run until a final time of 0.0075 seconds. We can change this for example to 0.005 by the following:

$ python elliptical_drop.py --tf=0.005

When this is run, PySPH will generate Cython code from the equations and integrators that have been provided, compiles that code and runs the simulation. This provides a great deal of convenience for the user without sacrificing performance. The generated code is available in ~/.pysph/source. If the code/equations have not changed, then the code will not be recompiled. This is all handled automatically without user intervention.

If we wish to run the code in parallel (and have compiled PySPH with Zoltan and mpi4py) we can do:

$ mpirun -np 4 /path/to/python elliptical_drop.py

This will automatically parallelize the run. In this example doing this will only slow it down as the number of particles is extremely small.

Visualizing and post-processing¶

You can view the data generated by the simulation (after the simulation is complete or during the simulation) by running the pysph_viewer application. To view the simulated data you may do:

$ pysph_viewer elliptical_drop_output/*.npz

If you have Mayavi installed this should show a UI that looks like:

On the user interface, the right side shows the visualized data. On top of it there are several toolbar icons. The left most is the Mayavi logo and clicking on it will present the full Mayavi user interface that can be used to configure any additional details of the visualization.

On the bottom left of the main visualization UI there is a button which has the text “Launch Python Shell”. If one clicks on this, one obtains a full Python interpreter with a few useful objects available. These are:

>>> dir()
['__builtins__', '__doc__', '__name__', 'interpolator', 'mlab',
 'particle_arrays', 'scene', 'self', 'viewer']
>>> len(particle_arrays)
1
>>> particle_arrays[0].name
'fluid'

The particle_arrays object is a list of ParticleArrays. The interpolator is an instance of pysph.tools.interpolator.Interpolator that is used by the viewer. The other objects can be used to script the user interface if desired.

from pysph.solver.utils import load
data = load('elliptical_drop_100.npz')

particle_arrays = data['arrays']
solver_data = data['solver_data']

fluid = particle_arrays['fluid']
p = fluid.p

from pysph.solver.utils import load
data = load('elliptical_drop_output/elliptical_drop_100.npz')
from pysph.tools.interpolator import Interpolator
parrays = data['arrays']
interp = Interpolator(parrays.values(), num_points=10000)
p = interp.interpolate('p')

from matplotlib import pyplot as plt
plt.contourf(interp.x, interp.y, p)

xmin, xmax, ymin, ymax, zmin, zmax = 0., 1., -1., 1., 0, 1
interp.set_domain((xmin, xmax, ymin, ymax, zmin, zmax), (40, 50, 1))
p = interp.interpolate('p')

interp.set_interpolation_points(x, y, z)

interp = Interpolator(parrays.values(), x=x, y=y, z=z)

The dam-break problem¶

The problem that is used for the illustration is the Weakly Compressible SPH (WCSPH) formulation for free surface flows, applied to a breaking dam problem:

A column of water is initially at rest (presumably held in place by some membrane). The problem simulates a breaking dam in that the membrane is instantly removed and the column is free to fall under it’s own weight and the effect of gravity. This and other variants of the dam break problem can be found in the examples directory of PySPH.

A non-PySPH implementation¶

We first consider the pseudo-code for the non-PySPH implementation. We assume we have been given two ParticleArrays fluid and solid corresponding to the dam-break problem. We also assume that an pysph.base.nnps.NNPS object nps is available and can be used for neighbor queries:

from pysph.base import nnps
fluid = get_particle_array_fluid(...)
solid = get_particle_array_solid(...)
particles = [fluid, solid]
nps = nnps.LinkedListNNPS(dim=2, particles=particles, radius_scale=2.0)

The part of the code responsible for the interactions can be defined as

class SPHCalc:
    def __init__(nnps, particles):
        self.nnps = nnps
        self.particles = particles

    def compute(self):
        self.eos()
        self.accelerations()

    def eos(self):
        for array in self.particles:
            num_particles = array.get_number_of_particles()
            for i in range(num_particles):
                array.p[i] =  # TAIT EOS function for pressure
                array.cs[i] = # TAIT EOS function for sound speed

    def accelerations(self):
        fluid, solid = self.particles[0], self.particles[1]
        nps = self.nps
        nbrs = UIntArray()

        # continuity equation for the fluid
        dst = fluid; dst_index = 0

        # source is fluid
        src = fluid; src_index = 0
        num_particles = dst.get_number_of_particles()
        for i in range(num_particles):

            # get nearest fluid neigbors
            nps.get_nearest_particles(src_index, dst_index, d_idx=i, nbrs)

            for j in nbrs:
                # pairwise quantities
                xij = dst.x[i] - src.x[j]
                yij = dst.y[i] - src.y[j]
                ...

                # kernel interaction terms
                wij = kenrel.function(xi, ...)  # kernel function
                dwij= kernel.gradient(xi, ...)  # kernel gradient

                # compute the interaction and store the contribution
                dst.arho[i] += # interaction term

        # source is solid
        src = solid; src_index = 1
        num_particles = dst.get_number_of_particles()
        for i in range(num_particles):

            # get nearest fluid neigbors
            nps.get_nearest_particles(src_index, dst_index, d_idx=i, nbrs)

            for j in nbrs:
                # pairwise quantities
                xij = dst.x[i] - src.x[j]
                yij = dst.y[i] - src.y[j]
                ...

                # kernel interaction terms
                wij = kenrel.function(xi, ...)  # kernel function
                dwij= kernel.gradient(xi, ...)  # kernel gradient

                # compute the interaction and store the contribution
                dst.arho[i] += # interaction term

        # Destination is solid
        dst = solid; dst_index = 1

        # source is fluid
        src = fluid; src_index = 0

        num_particles = dst.get_number_of_particles()
        for i in range(num_particles):

            # get nearest fluid neigbors
            nps.get_nearest_particles(src_index, dst_index, d_idx=i, nbrs)

            for j in nbrs:
                # pairwise quantities
                xij = dst.x[i] - src.x[j]
                yij = dst.y[i] - src.y[j]
                ...

                # kernel interaction terms
                wij = kenrel.function(xi, ...)  # kernel function
                dwij= kernel.gradient(xi, ...)  # kernel gradient

                # compute the interaction and store the contribution
                dst.arho[i] += # interaction term

We see that the use of multiple particle arrays has forced us to write a fairly long piece of code for the accelerations. In fact, we have only shown the part of the main loop that computes \(a_\rho\) for the continuity equation. Recall that our problem states that the continuity equation should evaluated for all particles, taking influences from all other particles into account. For two particle arrays (fluid, solid), we have four such pairings (fluid-fluid, fluid-solid, solid-fluid, solid-solid). The last one can be eliminated when we consider the that the boundary has zero velocity and hence the contribution will always be trivially zero.

The apparent complexity of the SPHCalc.accelerations method notwithstanding, we notice that similar pieces of the code are being repeated. In general, we can break down the computation for a general source-destination pair like so:

# consider first destination particle array

for all dst particles:
    get_neighbors_from_source()
    for all neighbors:
        compute_pairwise_terms()
        compute_inteactions_for_dst_particle()

# consider next source for this destination particle array
...

# consider the next destination particle array

The predictor-corrector integrator for this problem can be defined as

class Integrator:
    def __init__(self, particles, nps, calc):
        self.particles = particles
        self.nps = nps
        self.calc = calc

    def initialize(self):
        for array in self.particles:
            array.rho0[:] = array.rho[:]
            ...
            array.w0[:] = array.w[:]

   def stage1(self, dt):
       dtb2 = 0.5 * dt
       for array in self.particles:
           array.rho = array.rho0[:] + dtb2*array.arho[:]

           array.u = array.u0[:] + dtb2*array.au[:]
           array.v = array.v0[:] + dtb2*array.av[:]
           ...
           array.z = array.z0[:] + dtb2*array.az[:]

   def stage2(self, dt):
       for array in self.particles:
           array.rho = array.rho0[:] + dt*array.arho[:]

           array.u = array.u0[:] + dt*array.au[:]
           array.v = array.v0[:] + dt*array.av[:]
           ...
           array.z = array.z0[:] + dt*array.az[:]

   def integrate(self, dt):
       self.initialize()
       self.stage1(dt)   # predictor step
       self.nps.update()    # update NNPS structure
       self.calc.compute()  # compute the accelerations
       self.stage2(dt)   # corrector step

The Integrator.integrate method is responsible for updating the solution the next time level. Before the predictor stage, the Integrator.initialize method is called to store the values x0, y0... at the beginning of a time-step. Given the positions of the particles at the half time-step, the NNPS data structure is updated before calling the SPHCalc.compute method. Finally, the corrector step is called once we have the updated accelerations.

This hypothetical implementation can be integrated to the final time by calling the Integrator.integrate method repeatedly. In the next section, we will see how PySPH does this automatically.

PySPH implementation¶

Now that we have a hypothetical implementation outlined, we can proceed to describe the abstractions that PySPH introduces, enabling a highly user friendly and flexible way to define pairwise particle interactions.

We assume that we have the same ParticleArrays (fluid and solid) and NNPS objects as before.

equations = [

    # Equation of state
    Group(equations=[

            TaitEOS(dest='fluid', sources=None, rho0=ro, c0=co, gamma=gamma),
            TaitEOS(dest='boundary', sources=None, rho0=ro, c0=co, gamma=gamma),

            ]),

    Group(equations=[

            # Continuity equation
            ContinuityEquation(dest='fluid', sources=['fluid', 'boundary']),
            ContinuityEquation(dest='boundary', sources=['fluid']),

            # Momentum equation
            MomentumEquation(dest='fluid', sources=['fluid', 'boundary'],
                     alpha=alpha, beta=beta, gy=-9.81, c0=co),

            # Position step with XSPH
            XSPHCorrection(dest='fluid', sources=['fluid'])
            ]),
    ]

class SPHCalc:
    def __init__(nnps, particles):
        ...

    def compute(self):
        self.eos()
        self.accelerations()

for destination in particles:
    for equation in equations:
        equation.initialize(destination)

# This is where bulk of the computation happens.
for destination in particles:
    for source in destination.neighbors:
        for equation in equations:
            equation.loop(source, destination)

for destination in particles:
    for equation in equations:
        equation.post_loop(destination)

# Reduce any properties if needed.
total_mass = reduce_array(particles.m, 'sum')
max_u = reduce_array(particles.u, 'max')

class TaitEOS(Equation):
    def __init__(self, dest, sources=None,
                 rho0=1000.0, c0=1.0, gamma=7.0):
        self.rho0 = rho0
        self.rho01 = 1.0/rho0
        self.c0 = c0
        self.gamma = gamma
        self.gamma1 = 0.5*(gamma - 1.0)
        self.B = rho0*c0*c0/gamma
        super(TaitEOS, self).__init__(dest, sources)

    def loop(self, d_idx, d_rho, d_p, d_cs):
        ratio = d_rho[d_idx] * self.rho01
        tmp = pow(ratio, self.gamma)

        d_p[d_idx] = self.B * (tmp - 1.0)
        d_cs[d_idx] = self.c0 * pow( ratio, self.gamma1 )

class ContinuityEquation(Equation):
    def initialize(self, d_idx, d_arho):
        d_arho[d_idx] = 0.0

    def loop(self, d_idx, d_arho, s_idx, s_m, DWIJ, VIJ):
        vijdotdwij = DWIJ[0]*VIJ[0] + DWIJ[1]*VIJ[1] + DWIJ[2]*VIJ[2]
        d_arho[d_idx] += s_m[s_idx]*vijdotdwij

mat = declare("matrix((3,3))")
point = declare("cPoint")

cdef double[3][3] mat
cdef cPoint point

class FindMaxU(Equation):
    def reduce(self, dst):
        m = serial_reduce_array(dst.array.m, 'sum')
        max_u = serial_reduce_array(dst.array.u, 'max')
        dst.total_mass[0] = parallel_reduce_array(m, 'sum')
        dst.max_u[0] = parallel_reduce_array(u, 'max')

class EulerIntegrator(Integrator):
    def one_timestep(self, t, dt):
        self.initialize()
        self.compute_accelerations()
        self.stage1()
        self.do_post_stage(dt, 1)

class EulerStep(IntegratorStep):
    def initialize(self):
        pass
    def stage1(self, d_idx, d_u, d_v, d_w, d_au, d_av, d_aw, d_x, d_y,
                  d_z, d_rho, d_arho, dt=0.0):
        d_u[d_idx] += dt*d_au[d_idx]
        d_v[d_idx] += dt*d_av[d_idx]
        d_w[d_idx] += dt*d_aw[d_idx]

        d_x[d_idx] += dt*d_u[d_idx]
        d_y[d_idx] += dt*d_v[d_idx]
        d_z[d_idx] += dt*d_w[d_idx]

        d_rho[d_idx] += dt*d_arho[d_idx]

Working With Particles¶

As an object oriented framework for particle methods, PySPH provides convenient data structures to store and manipulate collections of particles. These can be constructed from within Python and are fully compatible with NumPy arrays. We begin with a brief description for the basic data structures for arrays

>>> import numpy
>>> from pyzoltan.core.carray import DoubleArray
>>> array = DoubleArray(10)                      # array of doubles of length 10
>>> array.set_data( numpy.arange(10) )           # set the data from a NumPy array
>>> array.get(3)                                 # get the value at a given index
>>> array.set(5, -1.0)                           # set the value at an index to a value
>>> array[3]                                     # standard indexing
>>> array[5] = -1.0                              # standard indexing

>>> import numpy
>>> from pysph.base.utils import get_particle_array
>>> x, y = numpy.mgrid[0:1:0.1, 0:1:0.1]             # create some data
>>> x = x.ravel(); y = y.ravel()                     # flatten the arrays
>>> pa = get_particle_array(name='array', x=x, y=y)  # create the particle array

>>> pa.remove_tagged_particles(tag=100)

>>> indices = [1,3,5,7]
>>> pa.remove_particles( indices )
>>> extracted = pa.extract_particles(indices, props=['rho', 'x', 'y'])

>>> pa.append_parray(another_array)

>>> props = ['au', 'av', 'aw']
>>> pa.set_to_zero(props)

>>> pa.add_constant('total_mass', 0.0)
>>> pa.add_constant('total_force', [0.0, 0.0, 0.0])
>>> print pa.total_mass, pa.total_force

>>> pa = ParticleArray(
...     name='test', x=x,
...     constants=dict(total_mass=0.0, total_force=[0.0, 0.0, 0.0])
... )

Nearest Neighbour Particle Searching (NNPS)¶

To carry out pairwise interactions for SPH, we need to find the nearest neighbours for a given particle within a specified interaction radius. The NNPS object is responsible for handling these nearest neighbour queries for a list of particle arrays:

>>> from pysph.base import nnps
>>> pa1 = get_particle_array(...)                    # create one particle array
>>> pa2 = get_particle_array(...)                    # create another particle array
>>> particles = [pa1, pa2]
>>> nps = nnps.LinkedListNNPS(dim=3, particles=particles, radius_scale=3)

The above will create an NNPS object that uses the classical linked-list algorithm for nearest neighbour searches. The radius of interaction is determined by the argument radius_scale. The book-keeping cells have a length of \(\text{radius_scale} \times h_{\text{max}}\), where \(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.

With these definitions, we can query for nearest neighbors like so:

>>> nbrs = UIntArray()
>>> nps.get_nearest_particles(src_index, dst_index, d_idx, nbrs)

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). 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 NNPS.update() method:

>>> nps.update()

>>> from pysph.base.nnps import DomainManager
>>> from pysph.base.point import Point
>>> domain = DomainManager(xmin, xmax, ymin, ymax, zmin, zmax,
                           periodic_in_x=True, periodic_in_y=True,
                           periodic_in_z=False)

class ParticleTAGS:
    Local = 0
    Remote = 1
    Ghost = 2

>>> x = numpy.array( [0, 1, 2, 3], dtype=numpy.float64 )
>>> tag = numpy.array( [0, 2, 0, 1], dtype=numpy.int32 )

>>> pa = utils.get_particle_array(x=x, tag=tag)

>>> print pa.get_number_of_particles()                     # total number of particles
>>> 4
>>> print pa.num_real_particles                            # no. of particles with tag 0
>>> 2

>>> x, tag = pa.get('x', 'tag', only_real_particles=True)  # get only real particles (tag == 0)
>>> print x
>>> [0. 2.]
>>> print tag
>>> [0 0]

>>> x, tag = pa.get('x', 'tag', only_real_particles=False) # get all particles
>>> print x
>>> [0. 2. 1. 3.]
>>> print tag
>>> [0 0 2 1]

Parallel NNPS with PyZoltan¶

PySPH uses the Zoltan data management library for dynamic load balancing through a Python wrapper PyZoltan, which provides functionality for parallel neighbor queries in a manner completely analogous to NNPS.

Particle data is managed and exchanged in parallel via a derivative of the abstract base class ParallelManager object. Continuing with our example, we can instantiate a ZoltanParallelManagerGeometric object as:

>>> ... # create particles
>>> from pysph.parallel import ZoltanParallelManagerGeometric
>>> pm = ZoltanParallelManagerGeometric(dim, particles, comm, radius_scale, lb_method)

The constructor for the parallel manager is quite similar to the NNPS constructor, with two additional parameters, comm and lb_method. The first is the MPI communicator object and the latter is the partitioning algorithm requested. The following geometric load balancing algorithms are supported:

• Recursive Coordinate Bisection (RCB)
• Recursive Inertial Bisection (RIB)
• Hilbert Space Filling Curves (HSFC)

The particle distribution can be updated in parallel by a call to the ParallelManager.update() method. Particles across processor boundaries that are needed for neighbor queries are assigned the tag Remote as shown in the figure:

Local and remote particles in the vicinity of a processor boundary (dashed line)

Putting it together: A simple example¶

Now that we know how to work with particles, we will use the data structures to carry out the simplest SPH operation, namely, the estimation of particle density from a given distribution of particles.

We consider particles distributed on a uniform Cartesian lattice ( \(\Delta x = \Delta y = \Delta\)) in a doubly periodic domain \([0,1]\times[0,1]\).

The particle mass is set equal to the “volume” \(\Delta^2\) associated with each particle and the smoothing length is taken as \(1.3\times \Delta\). With this initialization, we have for the estimation for the particle density

We will use the CubicSpline kernel, defined in pysph.base.kernels module. The code to set-up the particle distribution is given below

# PySPH imports
from pyzoltan.core.carray import UIntArray
from pysph.base.utils import utils
from pysph.base.kernels import CubicSpline
from pysph.base.nnps import DomainManager, LinkedListNNPS

# NumPy
import numpy

# Create a particle distribution
dx = 0.01; dxb2 = 0.5 * dx
x, y = numpy.mgrid[dxb2:1:dx, dxb2:1:dx]

x = x.ravel(); y = y.ravel()
h = numpy.ones_like(x) * 1.3*dx
m = numpy.ones_like(x) * dx*dx

# Create the particle array
pa = utils.get_particle_array(x=x,y=y,h=h,m=m)

# 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
k = CubicSpline(dim=2)

The NNPS object will create periodic ghosts for the particles along each periodic axis.

The ghost particles are assigned the tag value 2. For this example, periodic ghosts are created along each coordinate direction as shown in the figure.

class Kernel(object):
    def __init__(self, dim=1):
        self.radius_scale = 2.0
        self.dim = dim

    def kernel(self, xij=[0., 0, 0], rij=1.0, h=1.0):
        ...
        return wij

    def gradient(self, xij=[0., 0, 0], rij=1.0, h=1.0, grad=[0, 0, 0]):
        ...
        grad[0] = dwij_x
        grad[1] = dwij_y
        grad[2] = dwij_z

nbrs = UIntArray()                                                      # array for neighbors
x, y, h, m  = pa.get('x', 'y', 'h', 'm', only_real_particles=False)     # source particles will include ghosts

for i in range( pa.num_real_particles ):                                # iterate over all local particles
    xi = x[i]; yi = y[i]; hi = h[i]

    nps.get_nearest_particles(0, 0, i, nbrs)                            # get neighbors
    neighbors = nbrs.get_npy_array()                                    # numpy array of neighbors

    rho = 0.0
    for j in neighbors:                                                 # iterate over each neighbor

        xij = xi - x[j]                                                 # interaction terms
        yij = yi - y[j]
        rij = numpy.sqrt( xij**2 + yij**2 )
        hij = 0.5 * (h[i] + h[j])

        wij = k.kernel( [xij, yij, 0.0], rij, hij)                      # kernel interaction

        rho += m[j] * wij

    pa.rho[i] = rho                                                    # contribution for this destination

Summary¶

In this document, we introduced the most fundamental data structures in PySPH for working with particles. With these data structures, PySPH can be used as a library for managing particles for your application.

If you are interested in the PySPH framework and want to try out some eaxmples, check out the tutorials: Learning the ropes.

A simple example: Partitioning points¶

The following simple example demonstrates the partitioning of a random collection of points in the unit cube \([0,1]^3\). The objects to be partitioned in this case is therefore the points themselves which can be thought of as particles in an SPH simulation. We begin with some imports:

# Imports
import mpi4py.MPI as mpi
from pyzoltan.core.carray import UIntArray, DoubleArray
from pyzoltan.core import zoltan

from numpy import random
import numpy as np

from mpl_toolkits.mplot3d import Axes3D
import matplotlib.pyplot as plt

comm = mpi.COMM_WORLD
rank = comm.Get_rank()
size = comm.Get_size()

colors = ['r', 'g', 'b', 'y', 'm', 'k', 'c', 'burlywood']

numPoints = 1<<12

x = random.random( numPoints )
y = random.random( numPoints )
z = random.random( numPoints )
gid = np.arange( numPoints*size, dtype=np.uint32 )[rank*numPoints:(rank+1)*numPoints]

X = np.zeros( size * numPoints )
Y = np.zeros( size * numPoints )
Z = np.zeros( size * numPoints )
GID = np.arange( numPoints*size, dtype=np.uint32 )

comm.Gather( sendbuf=x, recvbuf=X, root=0 )
comm.Gather( sendbuf=y, recvbuf=Y, root=0 )
comm.Gather( sendbuf=z, recvbuf=Z, root=0 )

xa = DoubleArray(numPoints); xa.set_data(x)
ya = DoubleArray(numPoints); ya.set_data(y)
za = DoubleArray(numPoints); za.set_data(z)
gida = UIntArray(numPoints); gida.set_data(gid)

pz = zoltan.ZoltanGeometricPartitioner(
    dim=3, comm=comm, x=xa, y=ya, z=za, gid=gida)

pz.set_lb_method('RCB')
pz.Zoltan_Set_Param('DEBUG_LEVEL', '1')

Calling the load balance function¶

Once all the parameters are appropriately set-up, we can ask Zoltan to provide new assignments for the particles:

pz.Zoltan_LB_Balance()

This will call the chosen load balancing function internally and upon return, set a number of lists (arrays) indicating which objects need to be exported and which objects need to be imported. The data attributes for the export lists are:

• numExport : Number of objects to be exported
• exportLocalids : Local indices of the objects to be exported
• exportGlobalids : Global indices of the objects to be exported
• exportProcs : A list of size numExport indicating to which processor each object is exported

And similar arrays for the import lists. The import/export lists returned by Zoltan give an application all the information required to initiate the data transfer.

Note

Zoltan does not perform the data transfer. The data transfer must be done by the application or using the Unstructured communication utilities provided by Zoltan.

Given the new assignments, we once again broadcast this to the root to plot the final partition. The partition generated by this approach is shown below.

Point assignment to 4 processors where color indicates assignment.

We can see that the RCB method has resulted in cuts orthogonal to the domain axes. Each processor has exactly one fourth of the total number of particles.

The code for this example can be found in pyzoltan/core/tests/3d_partition.py

Inverting the Import/Export lists¶

In the example above, Zoltan returned a list of objects that were to be imported and exported. There arise situations in applications however, when only one set of arrays is available. For example, a common scenario is that we might know which objects need to be exported to remote processors but do not know in advance which objects need to be imported. The matter is complicated for dynamic applications because without a knowledge of the number of objects to be imported, we cannot allocate buffers of appropriate size on the receiving end.

For these scenarios when only one set (either import or export) of arrays is available, we use the PyZoltan.Zoltan_Invert_Lists() method to get the other set.

PyZoltan exposes this important utility function from Zoltan by assuming that the export lists are known to the application. Upon return from this method, the relevant import lists are also known. Note that the behaviour of import and export lists can be interchanged from the application.

A simple example demonstrating this is given below:

from pyzoltan.core import carray
from pyzoltan.core import zoltan

import numpy
import mpi4py.MPI as mpi

comm = mpi.COMM_WORLD; rank = comm.Get_rank(); size = comm.Get_size()

if rank == 0:
    proclist = numpy.array( [1, 1, 2, 1], dtype=numpy.int32 )
    locids = numpy.array( [1, 3, 5, 7], dtype=numpy.uint32 )
    glbids = numpy.array( [1, 3, 5, 7], dtype=numpy.uint32 )

if rank == 1:
    proclist = numpy.array( [0, 2, 0], dtype=numpy.int32 )
    locids = numpy.array( [1, 3, 5], dtype=numpy.uint32 )
    glbids = numpy.array( [11, 33, 55], dtype=numpy.uint32 )

if rank == 2:
    proclist = numpy.array( [1, 1], dtype=numpy.int32 )
    locids = numpy.array( [1, 3], dtype=numpy.uint32 )
    glbids = numpy.array( [111, 333], dtype=numpy.uint32 )

# create the Zoltan object
zz = zoltan.PyZoltan(comm)

# set the export lists
numExport = proclist.size; zz.numExport = numExport
zz.exportLocalids.resize(numExport); zz.exportLocalids.set_data(locids)
zz.exportGlobalids.resize(numExport); zz.exportGlobalids.set_data(glbids)
zz.exportProcs.resize(numExport); zz.exportProcs.set_data(proclist)

print 'Proc %d to send %s to %s'%(rank, glbids, proclist)

# Invert the lists
zz.Zoltan_Invert_Lists()

# get the import lists
numImport = zz.numImport
importlocids = zz.importLocalids.get_npy_array()
importglbids = zz.importGlobalids.get_npy_array()
importprocs = zz.importProcs.get_npy_array()

print 'Proc %d to recv %s from %s'%(rank, importglbids, importprocs)

In this example (which is hard coded for up to 3 processors), each processor artificially creates a list of objects it knows it must send to remote processors, which is set-up as the export lists for the PyZoltan object. Thereafter, PyZoltan.Zoltan_Invert_Lists() is called to get the lists that must be imported by each processor. The output from this example is shown below:

$ mpirun -n 3 python invert_lists.py
Proc 2 to send [111 333] to [1 1]
Proc 1 to send [11 33 55] to [0 2 0]
Proc 0 to send [1 3 5 7] to [1 1 2 1]
Proc 2 to recv [ 5 33] from [0 1]
Proc 0 to recv [11 55] from [1 1]
Proc 1 to recv [  1   3   7 111 333] from [0 0 0 2 2]

We can see that after a call to this method, each processor knows of remote data that must be received. In an application, this information can be used to effect the data transfer.

Another option is to use the unstructured communication utilities offered by Zoltan. This is described next.

Unstructured point to point communication¶

In the previous section, we saw how to use the Zoltan library function to invert a set of export indices to get corresponding import indices. Naturally, with a little bit of work, we can structure the requisite communication (MPI send and receives) to exchange the data.

This set of operations is fairly common and Zoltan (PyZoltan) provides a very convenient utility called ZComm to perform such communication. The typical use case for ZComm is when we know the list of local objects to send to remote processors but have no information about the objects to be received. An example (pyzoltan/core/tests/zcomm.py) demonstrating the use of the ZComm object is outlined below.

The example can be run with an arbitrary number of processors. Each processor allocates some data locally and randomly picks \(5\) of these objects to be sent to remote processors. The remote processors are also picked randomly:

import mpi4py.MPI as mpi
import numpy as np
from numpy import random

# import the unstructured communication package
from pyzoltan.core import zoltan_comm
from pyzoltan.core import zoltan

# MPI comm, rank and size
comm = mpi.COMM_WORLD; rank = comm.Get_rank(); size = comm.Get_size()

# each processor creates some random data
numObjectsTotal = 1<<10

x = random.random(numObjectsTotal)
gids = np.array( np.arange(size * numObjectsTotal) )[rank*numObjectsTotal:(rank+1)*numObjectsTotal]
gids = gids.astype(np.uint32)

# arbitrarily assign some objects to be sent to some other processor
nsend = np.int32( random.random_integers(low=1, high=10) )
object_ids = random.random_integers( low=0, high=numObjectsTotal, size=nsend )
proclist = random.random_integers(low=0, high=size-1, size=nsend).astype(np.int32)

my_indices = np.where(proclist == rank)[0]
proclist[my_indices] = (rank+1)%size

This information is sufficient enough to instantiate the ZComm object which will be used as the communication plan to exchange this data. Once the communication plan is setup, each processor knows of the data it must receive with the ZComm.nreturn attribute. This is used to allocate receive buffers:

# create the ZComm object
tag = np.int32(0)
zcomm = zoltan_comm.ZComm(comm, tag=tag, nsend=nsend, proclist=proclist)

# the data to send and receive
senddata = x[ object_ids ]
recvdata = np.ones( zcomm.nreturn )

With the send buffer and the allocated receive buffer, we can perform the communication using the ZComm.Comm_Do() method:

# use zoltan to exchange doubles
print "Proc %d, Sending %s to %s"%(rank, senddata, proclist)
zcomm.Comm_Do(senddata, recvdata)
print "Proc %d, Received %s"%(rank, recvdata)

Note that the user does not need to explicitly write the MPI send and receive calls. All that is required from the user is to correctly allocate the data on the receive side. The output from this example is (it will vary given the use of random numbers):

$ mpirun  -n 3 python zcomm.py
Proc 2, Sending [ 0.83476393  0.07041833  0.20059537  0.7722934   0.4529769 ] to [0 1 0 0 1]
Proc 2, Received [ 0.50391764  0.40028207]
Proc 0, Sending [ 0.50391764] to [2]
Proc 1, Sending [ 0.29755463  0.40028207  0.69433472] to [0 2 0]
Proc 1, Received [ 0.07041833  0.4529769 ]
Proc 0, Received [ 0.29755463  0.69433472  0.83476393  0.20059537  0.7722934 ]

senddata = gids[ object_ids ]
recvdata = np.ones(zcomm.nreturn, dtype=np.uint32)
zcomm.set_nbytes(recvdata.itemsize)

print "Proc %d, Sending %s to %s"%(rank, senddata, proclist)
zcomm.Comm_Do(senddata, recvdata)
print "Proc %d, Received %s"%(rank, recvdata)

$ mpirun  -n 3  python zcomm.py
Proc 1, Sending [1054 1692 2034] to [0 2 0]
Proc 0, Sending [214] to [2]
Proc 2, Sending [2720 3034 2511 2412 2975] to [0 1 0 0 1]
Proc 2, Received [ 214 1692]
Proc 1, Received [3034 2975]
Proc 0, Received [1054 2034 2720 2511 2412]

recvdata[:] = rank

updated_info = np.zeros(zcomm.nsend, dtype=senddata.dtype)
print 'Proc %d, sending updated data %s'%(rank, recvdata)
zcomm.Comm_Do_Reverse(recvdata, updated_info)
print 'Proc %d, received updated data %s'%(rank, updated_info)

$ mpirun  -n 3  python zcomm.py
Proc 1, Sending [1054 1692 2034] to [0 2 0]
Proc 0, Sending [214] to [2]
Proc 2, Sending [2720 3034 2511 2412 2975] to [0 1 0 0 1]
Proc 2, Received [ 214 1692]
Proc 1, Received [3034 2975]
Proc 0, Received [1054 2034 2720 2511 2412]
Proc 0, sending updated data [0 0 0 0 0]
Proc 2, sending updated data [2 2]
Proc 1, sending updated data [1 1]
Proc 0, received updated data [2]
Proc 1, received updated data [0 2 0]
Proc 2, received updated data [0 1 0 0 1]

Distributed data directories¶

The Zoltan Distributed Data Directory utility is a convenient way for a processor to locate remote data. It is implemented as a parallel hash map, keyed on the object identifiers (global indices) and with arbitrary user data associated with each entry.

The use of this feature is highly problem dependent since the user defined data will necessarily change for different applications. We use a simple example demonstrating it’s use. Each processor stores ownership of the object in the distributed directory without any user data associated with each entry.

We begin with the standard set of imports and create some data on each processor and assign each object a unique global identifier:

import numpy
import pyzoltan.api as pz
import mpi4py.MPI as mpi

comm = mpi.COMM_WORLD
rank = comm.Get_rank()
size = comm.Get_size()

# every processor owns some data
numObjectsTotal = 5
my_indices = numpy.array( range( rank*numObjectsTotal,(rank+1)*numObjectsTotal ), dtype=numpy.uint32 )

gid = pz.UIntArray(my_indices.size); gid.set_data( my_indices )

Additionally, each processor has an IntArray which denotes object assignment:

part_assignment = numpy.array( [rank]*numObjectsTotal, dtype=numpy.int32 )
part = pz.IntArray( part_assignment.size ); part.set_data( part_assignment )

This is sufficient data to create the distributed directory:

# create a zoltan dd and store the object assignments
dd = pz.Zoltan_DD(comm)

# update the dd with the data
dd.Zoltan_DD_Update(gid, part)

Note that after instantiation of the Zoltan_DD object, we call the Zoltan_DD.Zoltan_DD_Update() method to update the data associated with this directory. Now, given the shared data available with each processor, we can query for object assignments. In the example below, each processor queries for the objects with global indices numObjectsTotal + rank and numObjectsTotal - rank:

# now we can query the dd
owner_data = pz.IntArray()   # output array for the object data assignment
owner_parts = pz.IntArray()  # output array for the object assignment

# every processor requests for information about some data
query_gids = pz.UIntArray(2); query_gids.set_data( numpy.array([numObjectsTotal+rank,
                                                                numObjectsTotal-rank], dtype=numpy.uint32) )

# use Zoltan_DD_Find to query the data
dd.Zoltan_DD_Find(query_gids, owner_parts, owner_data)

The result from this quey with \(3\) processors is shown below:

$ mpirun  -n 3 python zoltan_dd.py
Processor 0, query_gids = [5 5], owner_parts = [1 1], owner_data = [1 1]
Processor 1, query_gids = [6 4], owner_parts = [1 0], owner_data = [1 0]
Processor 2, query_gids = [7 3], owner_parts = [1 0], owner_data = [1 0]