Overview

Modern computer architectures work most efficiently with memory references that are localized in time and predictable in space. Particle simulations typically do neither. Because the particles can move in unpredictable ways, the memory reference pattern they generate is often random, resulting in inefficient use of the hardware. Many optimizations are known for this problem, but they often do not get incorporated in real codes because they are complex and may obscure the physics being modeled by the codes.

The goal of the PST is to simplify the construction of complex particle codes by providing a high-level interface and to allow the time and space localization optimizations to be separated from the numerical algorithms and physical models that are represented by the code.

This page has a quick summary of the project, and these links have more detailed information about the PST. This document is focused on describing functionality rather than syntax. A detailed description of the syntax is available with the library itself.

• Interoperabilty
• Computations Using Particles
• Particles object
• Attributes object
• Calcuating with particles
• Layouts and Load Balancing

Approach

The key components of the PST are the Particles base class and the Attribute class. The Particles class contains a list of Attribute objects describing the particles to be modeled. Each Attribute is a distributed vector of elements of arbitrary type, with one element to hold the value of that variable for each particle. The Particles class provides methods to access each Attribute and to add or remove individual particles. The user can create his or her own unique type of particles representation by creating a new class derived from Particles that has specific types of Attribute objects as data members. This new user-defined class can register its Attribute objects and then all of the methods of Particles will operate on each of the Attributes of this new class, making this representation completely general. For some common cases, we will provide sample classes derived from Particles that provide arrays of built-in scalar types and can be constructed directly in C++ or via a C or Fortran interface.

The Attribute data distribution across multiple processors will be managed by one of several Layout classes that we will provide. The Layout classes will implement different parallel decomposition strategies. The most common of these strategies is a "spatial layout", in which the particles are assumed to move within a block-decomposed spatial domain and we keep particle data locality with any field data that exists on this domain. This is done so that when particles interact with field quantities (i.e., via a gather or scatter operation), such operations are local to a processor and therefore fast. Once Layout objects are set up, data distribution can handled automatically with a simple function call, and we will provide several load balancing classes that can repartition the domain based on different optimization strategies.

There are a few different ways that the user can interact with Attribute data in the PST.

• As far as the Particles class is concerned, it just contains a set of Attributes, and it manages them to make sure that they all represent the same number of particles with the same distribution.
• The Particles class then provides access to individual Attribute elements so the user can manipulate those arrays however they wish. In many situations, including all C and FORTRAN uses and some C++ uses, the user will simply loop over these local arrays to manipulate the particle data.
• In C++, the Attribute class can be "expression enabled" using, for example, the Portable Expression Template Engine, or PETE, so that one can write simple data-parallel mathematical expressions involving Attributes and scalar types that are evaluated very efficiently.
• The user may impose "boundary conditions" on an Attribute that take a particular action when a particle's value for that Attribute goes outside a specified range.
• For users writing particle-in-cell (PIC) simulations, we will provide a simple interface for interpolating data between particle positions and mesh locations with a choice of interpolation schemes for both writing to a mesh and reading from a mesh.
• For simulations involving the computation of interactions between nearby particles, an Attribute is capable of maintaining a separate list of data for particles stored with sub-domains in non-local memory. That way, interaction force calculations involving spatially-neighboring particles can be performed efficiently on parallel computers.

Tool Availability

Because the PST and PETE make extensive use of the template mechanism in C++ and templates are a relatively advanced feature, not all compilers will be able to compile the library. The KCC compiler from KAI and the EGCS compiler from Cygnus work well on all Unix platforms for advanced C++ source code, along with the SGI CC compiler on Irix platforms, and CodeWarrior from Metrowerks and the Intel C++ compiler on Macintosh and Windows platforms. For users that only link against the library (such as C and Fortran users) there are no restrictions on the compiler to use once the library itself is compiled.

Participants

Why is interoperability so hard and so rare?

• Data formats. The LINPACK and LAPACK libraries have been very successful, and a big reason for that is that everybody knows and understands the data format - FORTRAN order two dimensional arrays - and that format is well suited to the computation at hand. In many cases the fact that different packages have their data in conflicting formats is a fatal barrier to having them interoperate.
• Tools. If dense arrays of numbers in a single memory space with a single thread of control what all you need, the primary tool of the scientific developer - the FORTRAN language - does a good job. Other languages (such as Matlab, C++ and Mathematica for different types of tasks) have only recently gotten to the point where they are really useful for mainstream scientific computing. The state of tools for scientific computing is vastly more primitive than it is for GUI or even Web programming.
• Production versus prototype software. Producing code that someone else can use out of the box on a problem different from one that the original author has done is much harder than producing code only for the developer's own use. Industrial estimates of the cost of moving prototype code to production range from a factor of three to a factor of ten above producing the prototype.
• Reward systems. Individual scientists are usually not trained in writing reusable software, are rarely rewarded for that skill, and projects are rarely rewarded for taking a software project past the prototype stage. A prototype will get you a paper, and taking the code from there to production quality will rarely generate a new paper.

External Polymorphism and Traits

The crux is that the library generally doesn't just specify the generic operations, it specifies very specific syntax for how to do those operations. There are many subtle choices to be made in the syntax and the precise meaning of each syntax construction, and it is on this sort of mundane issue that much interoperability founders.

External polymorphism and traits are specific techniques for encapsulating an external data structure and allowing the user to supply their own data structure along with instructions on how to use it. They greatly simplify the task of coming up with a consistent interface for the PST to use and which a user can supply. By using these techniques, the actual data format can be a user choice, not a library choice, and we think that is a key step toward routine library interoperability.

Leveraging

In order to leverage off of our previous work in the POOMA, SMARTS, and other libraries in the ACL, we are in the process of factoring those libraries so that the PST (and other packages) can use the capabilities it needs without importing capabilities it does not need. For example, the PST does not need the complexities of expression evaluation in C++, but it does need the domain calculus and communication layer. By factoring POOMA so it can use only what it needs, we can leverage previous work without burdening customers with more than they need.

Object Oriented Programming and Generics

Modern C++ provides object oriented mechanisms (inheritance, virtual functions, etc.) and generic mechanism (templates) as well as the low level elements necessary for programming close to the machine (pointers, bitwise operations) necessary for representing complex data structures and algorithms efficiently. The sort of separation that is achieved with Arrays and Engines cannot be done without appropriate tools.

That is not to say that only sophisticated C++ programmers can take advantage of a library such as the PST. The library itself can be written using these techniques, and then export a C or FORTRAN interface. The calling program then need not care what language the library is written in. The C++ user would have a richer or more convenient interface with which to manipulate the PST, but the functionality is available more widely. 

Particle-Particle

Some typical examples of this sort of computation are Smooth Particle Hydrodynamics and many kinds of semi-classical molecular simulations.

Particle-Mesh

Plasma simulations are a very common form of PIC because the interaction length scale is typically much longer than mean distance between particles, so the force on any given particle is the sum of many small interactions and is well represented by a smooth field represented on a mesh.

Combinations and Complexities

At high enough plasma densities, the hard interactions in a plasma simulation are important because they represent fast momentum and energy transfers that otherwise would not be properly handled. Direct Simulation Monte Carlo (DSMC) is a technique for smoothly matching a fluid simulation taking place purely on a grid to a surface with a length scale comparable to the mean free path, and through that transition both grid and particle quantities must be handled in consistent ways.

Because there are a huge number of different ways in which these computations need to be constructed, a library to facilitate them must be careful not to box in the user. It must provide generic mechanisms for fundamental operations (gather/scatter, select subsets of particles, loop over particles, define attributes of particles, etc) rather than seemingly high level but inflexible "deposit charge" functions. These low level operations can then be combined flexibly and recursively to build common higher level operations while allowing users to open the hood write arbitrary algorithms while still taking advantage of the bookkeeping functions for operations like moving particles between processors using Layouts. 

There are a couple of possible ways to handle that, and the PST team has explored the implications of them for performance, flexibility and generality. Our decision has been to use what is often called the "multiple arrays" technique, which means that the Particles object holds an array for each attribute rather than a single array of particle objects, where each particle object has all of the attributes for each particle.

This technique has the following advantages:

• Performance. Many particle calculations don't involve all of the particle attributes in any one loop. On a cache based machine if the attributes for each particle are stored together, data that will not get used in a particular calculation will have to be loaded into cache, which is a performance bottleneck. On vector machines, keeping the attributes for each particle together imply a non-unit stride for any given attribute, which often results in performance degradation.
• Flexibility. If the attributes for each particle were stored together you would not be able to vary the specific attributes in a Particles class at run time. With separate arrays for each attribute, you can add and delete attributes at will.
• Generality. If a given code really would benefit from storing the particle attributes together, this behavior can be recovered by having a single attribute which is a composite object that contains all of the attributes for each particle.

Synchronization

This synchronization is maintained by making sure that all of the operations available in the Particles class (add, delete, reorder, etc.) are synchronization preserving.

Data Registration and Access

Several of these concrete subclasses are provided, both to be used directly and as examples of how to design your own if the provided ones do not have the necessary functionality.

• Staticly defined attributes. If a given application knows that it will never have to add and delete attribute arrays, simplifications can be made in the implementation and the interface.
• Dynamically defined attributes. When the user requires the ability to dynamically add and delete attributes, the more complex subclass would be used.
• Different attribute representations. There are FORTRAN arrays that cannot be resized, C/C++ arrays that can be resized, there are standard C++ classes like std::vector, there are any number of more complex array classes that internally handle parallelism, such as the POOMA Array classes. Different subclasses of Particles can work with different data representations.

The Attribute object provides two kinds of services:

• Synchronization control. This is the interface the Particles object sees, and which it uses to keep the various Attributes it is controlling in synch. Array resizing operations are fundamental here, and the numerical user is not intended to use these controls directly.
• Element access. This is the interface the user sees to manipulate the data in the Attribute to perform the computation. This usually comes down to random access on a one dimensional array. This interface is also used by the Particles class to insert and delete elements.

Generic Interface

The PST resolves this tension by giving the user the ability to describe the format for the data, and prebuilding descriptions for some common cases. It places few restrictions on the details of the interface for the Attribute object the user supplies, and instead has the user provide a traits class to describe the interface.

For example, some arrays are not resizable, and those that are can have different interfaces to actually resize it. The traits class the user supplies records whether it is resizable and if it is, what is the function to call to resize it.

C and FORTRAN Interface

Conventional FORTRAN 77 arrays, for example, are not resizable, and allocatable arrays have a well defined interface. In C most arrays are manipulated with malloc/free/realloc, and you could allow the user to supply their own functions for each of those. In C++ there are the containers in the standard library which have their own interface, and parallel POOMA arrays with still another. By packaging the interface specification in a user supplied traits class and using that to supply a number of prebuilt ones, ease of use and extensibility are provided. 

• Get out of the way. A primary goal of any advanced library must be to be sure to get out of the way of doing simple things. For the PST that means that it has to be able to only add capabilities to the users' data structures, so that the they can do everything they could do without the PST, plus new things.
• Encapsulate the bookkeeping. Often the most complex part of a code that does particle computations is the bookkeeping of packing messages, unpacking them, sorting the particles for cache performance, looping over subsets of the particles, and so on. There is no way of hiding from the user that these things have to happen, but you can encapsulate the mechanics of doing it.
• Prefer several small orthogonal tools to one monolithic integrated one. Whenever you can break a large tool up into several smaller ones that work together but don't overlap, you increase the flexibility available to the user.

Sort, Load Balance, Create, Destroy

Gather/Scatter with Stencils

DIY: They're your arrays. Use them.

Particle Layout

A couple of prepackaged particle layout objects are the trivial one which never moves any particles, and a spacial layout which defines mappings between regions of physical space and CPUs. These are fairly abstract representations of how to map particles to CPUs.

Attribute Layout