Contents:
Introduction
Overview
Attributes
Layout
Derivation
Synchronization and Related Issues
Example: Simple Harmonic Oscillator
Boundary Conditions
Example: Elastic Collision
Summary
Particles are primarily used in one of two ways in large scientific applications. The first is to track sample particles using Monte Carlo techniques, for example, to gather statistics that describe the conditions of a complex physical system. Particles of this kind are often referred to as tracers. The second is to perform direct numerical simulation of systems that contain discrete point-like entities such as ions or molecules.
In both scenarios, the application contains one or more sets of particles. Each set has some data associated with it that describes its members' characteristics, such as mass and charge. Particles typically exist in a spatial domain, and they may interact directly with one another or with field quantities defined on that domain.
This tutorial gives an overview of POOMA's support for particles, then discusses some implementation details. The classes introduced in this tutorial are illustrated by two short programs: one that tracks particles under the influence of a simple one-dimensional harmonic oscillator potential, and another that models particles bouncing off the walls of a closed three-dimensional box. The next tutorial then shows how particles and fields can be combined to create complete simulation applications.
POOMA's Particles class is a container for a heterogeneous collection of particle attributes. The class uses dynamic storage for particle data (in the form of DynamicArrays), so that particles can be added or deleted as necessary. It contains a layout object that manages the distribution of particle data across multiple patches, and it applies boundary conditions to particles when attribute data values exceed a prescribed range. In addition, global functions are provided for interpolating data between particle and field element positions.
Each Particles object keeps a list of pointers to its elements' attributes. When an application wants to add or delete a particle, it invokes a method on the Particles object, which delegates the call to the layout object for the contained attributes. Particles also provides a member function called sync(), which the application invokes in order to update the particle count and data distribution across contexts and to apply the boundary conditions.
Applications can define a new type of particles collection by deriving from the Particles class. The derived class declares data members for the attributes needed to characterize this type of particle; the types of these data members are discussed below. The constructor for this class calls Particles::addAttribute to register each attribute and add it to the list. In this way, the Particles class can be extended by the application to accommodate any sort of particle description.
The distribution of particle data stored in DynamicArrays is directed by a particle layout class. (The details of the mechanism used to specify layout and other information for Particles classes are discussed below.) Each particle layout class employs a particular strategy to determine the patch in which a particle's data should be stored. For instance, SpatialLayout keeps each particle in the patch that contains field data for elements that are nearest to the particle's current spatial position. This strategy is useful for cases where the particles need to interact with field data or with particles nearby to them.
Each particle attribute is implemented as a DynamicArray, a class derived from the one-dimensional specialization of POOMA's Array class. DynamicArray extends the notion of a one-dimensional array to allow applications to add or delete elements at will. When particles are destroyed, the empty slots left behind can be filled by moving elements from the end of the list (backfill) or by sliding all the remaining elements over and preserving the existing order (shift up). At the same time, DynamicArrays can be used in data-parallel expressions in the same way as ordinary Arrays, so that the application can update particle attributes such as position and velocity using either a single statement or a loop over individual particles.
At first glance, it might seem more sensible to have applications define a type T that stores all the attribute data for one particle in a single data structure, and then use this as a template argument to the Particles class, which would store a DynamicArray of values of this type. POOMA's designers considered this option, but discarded it. The reason is that most compute-intensive operations in scientific applications are implemented as loops in which one or more separate attributes are read or written. In order to make the evaluation of expressions involving attributes as efficient as possible, it is therefore important to ensure that data are arranged as separate one-dimensional arrays for each attribute, rather than as a single array of structures with one structure per particle. This arrangement makes common cases such as:
for (int i=0; i<n; ++i)
{
x[i] += dt * vx[i];
y[i] += dt * vy[i];
}
run more quickly, as it makes much better use of the cache.
As mentioned above, each Particles object uses a layout object to determine in which patch a particle's data should be stored. The layout manages the program's requests to re-arrange particle data. With SpatialLayout, for example, the application provides a particle position attribute which is used to determine how particle data should be distributed. The particle layout then directs the Particles object to move particle data from one patch to another as dictated by its strategy. The Particles object in turn delegates this task to the layout object for the particle attributes, which tells each of the attributes using this layout to move their data as needed. All of this is handled by a single call to Particles::sync(), which in turn calls Particles::swap() to actually move particle data around.
In general, creating a new Particles class is a three-step process. The first step is to declare a traits class whose typedefs specify the type of engine the particle attributes are to use and the way the data for those attributes is to be distributed. An example of such a traits class is the following:
struct MyParticleTraits
{
typedef MultiPatch<GridTag,Brick> AttributeEngineTag_t;
typedef UniformLayout ParticleLayout_t;
};
This traits class will be used to specialize the Particles class template when an application class representing a concrete set of particles is derived from it. Particles uses public typedefs to give sensible names to these traits parameters, so that the derived application-level class can access them (as shown below). For the application developer's convenience, a set of pre-defined particle traits classes with specific choices of attribute engine and particle layout type are provided in the header file src/Particles/CommonParticleTraits.h. These define combinations of shared brick and multi-patch brick engines with both uniform and spatial layouts, and include the following:
Name AttributeEngineTag_t ParticleLayout_t SharedBrickUniform SharedBrick UniformLayout SharedBrickSpatial <Cent,Mesh,FieldLayout> SharedBrick SpatialLayout< DiscreteGeometry<Cent,Mesh>, FieldLayout> MPBrickUniform MultiPatch<GridTag,Brick> UniformLayout MPBrickSpatial <Cent,Mesh,FieldLayout> MultiPatch<GridTag,Brick> SpatialLayout< DiscreteGeometry<Cent,Mesh>, FieldLayout>
The SharedBrick engine type is just like Brick, except that the engine's layout can be shared by other engines constructed with the same layout argument. The effect of this is that the layout of all of the attributes remains synchronized. SharedBrick should only be used when running serially; otherwise, applications should use MultiPatch.
The second step is to derive a class from Particles. The new class can be templated on whatever the developer desires, as long as a traits class type is provided for the template parameter of the Particles base class. In the example below, the new class being derived from Particles is templated on the same traits class as Particles. For the sake of convenience, typedefs may be provided for the instantiated parent class and for its layout type. The constructor for the application class then usually takes a concrete layout object of the type specified in the typedef above as a constructor argument:
template <class PT>
class MyParticles : public Particles<PT>
{
public:
// instantiated type of parent class
typedef Particles<PT> Base_t;
// type of layout (from traits class via parent class)
typedef typename Base_t::ParticleLayout_t ParticleLayout_t;
// type of attribute engine tag (from traits class via parent class)
typedef typename Base_t::AttributeEngineTag_t EngineTag_t;
// some particle attributes as public data members
DynamicArray<double, EngineTag_t> charge;
DynamicArray<double, EngineTag_t> mass;
DynamicArray<int, EngineTag_t> count;
// constructor invokes Particles(layout) to cache layout
MyParticles(const ParticleLayout_t &layout)
: Particles<PT>(layout)
{
// register attributes
addAttribute(charge);
addAttribute(mass);
addAttribute(count);
}
};
Note that the attribute elements in this example have different element types, i.e., charge and mass are double, while count is int. Attribute elements may in general have any type, including any user-defined type.
Finally, the application class MyParticles is instantiated with the traits class MyParticleTraits to create an actual set of particles. An actual layout is declared first, and it is passed as a constructor argument to the instance of the application-level class to control the distribution of particle data between patches. This layout object typically has one or more constructor arguments that specify such things as the number of patches the particles are to be distributed over:
int main()
{
const int numPatches = 10;
MyParticleTraits::ParticleLayout_t layout(numPatches);
MyParticles<MyParticleTraits> particles(layout);
}
While this may seem complex at first, each level of indirection or generalization is needed in order to provide flexibility. The type of engine and layout to be used, for example, could be passed directly as template parameters to Particles, rather than being combined together in a traits class. However, this would make user-level code fragile in the face of future changes to the library: if other traits are needed later, they can be added to the traits class in one place, rather than needing to be specified every time something is derived from Particles. This bundling also makes it easier to specify the same basic properties (engine and layout) for two or more interacting Particles-derived classes.
For efficiency reasons, Particles does not automatically move particle data between patches after every operation, but instead waits for the application to call the method sync(). Particles can also be configured to cache requests to delete particles, rather than deleting them immediately.
Particles::sync() is a member template, i.e., it is templated on its single argument. This argument must be one of the particle set's attributes. SpatialLayout assumes that the attribute given to sync() is the particles' positions, and uses it to update the distribution of particle data so that particles are located on the same patch as nearby field data. Applications must therefore be careful not to mistakenly pass a non-spatial attribute, such as temperature or pressure, to SpatialLayout.
UniformLayout, which divides particles as evenly as possible between patches, without regard for spatial position, only uses the attribute passed to sync() as a template for the current distribution of particle data. Any attribute with the same distribution as the actual particle data can therefore be used.
The use of a parameter in Particles::sync() is one important difference between the implementation of particles in this version of POOMA and its predecessor. In the old design, all Particles classes came with a pre-defined attribute R that was the particles' position, which synchronization always referred to. The new scheme allows applications to switch the attribute that is used to represent the position, e.g., to switch back and forth between a "current" position attribute currpos and a "new" position attribute newpos. It also allows particles to be weighted according to some attribute, so that the distribution scheme load-balances by weight.
Of course, before particle data can be (re-)distributed, the particles themselves must be created. Particles provides two methods for doing this. The first, globalCreate(num,renum), creates a specified number of particles, spread as evenly as possible across all patches. The particles are normally renumbered after the creation operation, although this can be overridden by passing false as a second parameter to the method.
Particles::create(num, patch, renum), on the other hand, creates a specified number of particles within the local context, and adds them to either the last local patch (if the patch argument is negative) or to a specific patch (if patch is non-negative). The particles are renumbered after this operation unless false is passed as a third parameter to this method.
After particles have been created (or destroyed), they must be renumbered to ensure that each has a unique ID. In general, the renumber() method surveys all the patches to find out what the current local domain of each patch is. It then reconstructs a global domain across all the patches, effectively renumbering the particles from 0 to N-1, where N is the total number of particles. The more complex sync() method applies the particle boundary conditions, performs any deferred particle destroy requests, swaps particles between patches according to the particle layout strategy, and then renumbers the particles by calling renumber(). Programs should therefore call renumber() if they have only created or destroyed particles, but have not done deferred destroy requests, modified particle attributes in a way that would require applying boundary conditions (or have no boundary conditions), and do not need to swap particles.
If a program does not (implicitly or explicitly) call renumber() after creating or destroying particles, the global domain for the particles will be incorrect. If the program then tries to read or write a view of a particle attribute by indexing with some domain object, it will not get the right section of the data. This failure could be silent if the view that the program requests exists. Alternatively, the requested view could be outside of the global domain (because renumber() was not called to update the global domain), in which case the layout object for the particle attribute will suffer a run-time assertion failure.
There are also two ways to destroy particles. The first way, which always destroys the particles immediately, is implemented by the method Particles::destroy(domain, patchId, renum). If the patchId parameter is negative (which is the default), the domain is assumed to specify a global numbering of particles. If patchId is non-negative, then domain is assumed to be a local numbering for that patch, i.e., one in which the first particle in the patch has index 0.
Since this method modifies the Particles object right away, the default behavior of this method is to renumber particles after it has finished destroying the specified particles. This can be overridden by passing false as the last parameter to the call.
The second particle destruction method is Particles::deferredDestroy(domain, patch). This is new in this release, and only does deferred destruction, i.e., only caches the requested indices for use later when performDestroy() is called. (Since this method doesn't actually destroy particles right away, there is no need for it to call renumber(). The performDestroy() method, which causes the cached destruction requests to be executed, always performs renumbering.)
As noted above, Particles::globalCreate() normally calls renumber() to update the global domain of the particle attributes after the particles have been created, but before the program tries to do computations involving their attributes. The reason for this is that while globalCreate() allocates space for the new particle data and updates the local domain of the patch or patches on which creation was done, the global domain across all the patches of data is not updated until the call to renumber(). If the global domain is not up to date, the program cannot correctly access the ith particle's data or evaluate a data-parallel expression.
The example for this tutorial simulates the motion of particles under the influence of a simple one-dimensional harmonic oscillator potential. The code, which is included in the release in the examples/Particles/Oscillation directory, is as follows:
001 #include <iostream>
002 #include <stdlib.h>
003
004 #include "Pooma/Particles.h"
005 #include "Pooma/DynamicArrays.h"
006
007 // Dimensionality of this problem
008 static const int PDim = 1;
009
010 // A traits class specifying the engine and layout of a Particles class.
011 template <class EngineTag>
012 struct PTraits
013 {
014 // The type of engine to use in the particle attributes.
015 typedef EngineTag AttributeEngineTag_t;
016
017 // The type of particle layout to use. Here we use a UniformLayout,
018 // which divides the particle data up so as to have an equal number
019 // on each patch.
020 typedef UniformLayout ParticleLayout_t;
021 };
022
023 // A Particles subclass that defines position and velocity as
024 // attributes.
025 template <class PT>
026 class Quanta : public Particles<PT>
027 {
028 public:
029 // Useful things to extract from the base class
030 typedef Particles<PT> Base_t;
031 typedef double AxisType_t;
032 typedef typename Base_t::ParticleLayout_t ParticleLayout_t;
033 typedef typename Base_t::AttributeEngineTag_t AttributeEngineTag_t;
034 enum { dimensions = PDim };
035
036 // Constructor sets up layouts and registers attributes
037 Quanta(const ParticleLayout_t &pl)
038 : Particles<PT>(pl)
039 {
040 addAttribute(x);
041 addAttribute(v);
042 }
043
044 // X position and velocity are attributes (would normally be
045 // private, with accessor methods)
046 DynamicArray< Vector<dimensions, AxisType_t>, AttributeEngineTag_t > x;
047 DynamicArray< Vector<dimensions, AxisType_t>, AttributeEngineTag_t > v;
048 };
049
050 // Engine tag type for attributes. Here we use a MultiPatch engine
051 // with the patches being Bricks of data, and a GridTag, which allows
052 // the patches to possibly be of differing sizes. This is important
053 // since we may not have the same number of particles in each patch.
054 typedef MultiPatch<GridTag, Brick> AttrEngineTag_t;
055
056 // The particle traits class and layout type for this application
057 typedef PTraits<AttrEngineTag_t> PTraits_t;
058 typedef PTraits_t::ParticleLayout_t PLayout_t;
059
060 // Simulation control constants
061 const double omega = 2.0;
062 const double dt = 1.0 / (50.0 * omega);
063 const int nParticle = 100;
064 const int nPatch = 4;
065 const int nIter = 500;
066
067 // Main simulation routine.
068 int main(int argc, char *argv[])
069 {
070 // Initialize POOMA and Inform object for output to terminal
071 Pooma::initialize(argc, argv);
072 Inform out(argv[0]);
073 out << "Begin Oscillation example code" << std::endl;
074
075 // Create a uniform layout object to control particle positions.
076 PLayout_t layout(nPatch);
077
078 // Create Particles, using our special subclass and the layout
079 typedef Quanta<PTraits_t> Particles_t;
080 Particles_t p(layout);
081
082 // Create particles on one patch, then re-distribute (just to show off)
083 p.create(nParticle, 0);
084 for (int ip=0; ip<nPatch; ++ip)
085 {
086 out << "Current size of patch " << ip << " = "
087 << p.attributeLayout().patchDomain(ip).size()
088 << std::endl;
089 }
090
091 out << "Resyncing particles object ... " << std::endl;
092 p.sync(p.x);
093
094 // Show re-balanced distribution.
095 for (int ip=0; ip<nPatch; ++ip)
096 {
097 out << "Current size of patch " << ip << " = "
098 << p.attributeLayout().patchDomain(ip).size()
099 << std::endl;
100 }
101
102 // Randomize positions in domain [-1,+1], and set velocities to zero.
103 // This is done with a loop because POOMA does not yet have RNGs.
104 typedef Particles_t::AxisType_t Coordinate_t;
105 Vector<PDim, Coordinate_t> initPos;
106 srand(12345U);
107 Coordinate_t ranmax = static_cast<Coordinate_t>(RAND_MAX);
108 for (int ip=0; ip<nParticle; ++ip)
109 {
110 for (int idim=0; idim<PDim; ++idim)
111 {
112 initPos(idim) = 2 * (rand() / ranmax) - 1;
113 }
114 p.x(ip) = initPos;
115 p.v(ip) = Vector<PDim, Coordinate_t>(0.0);
116 }
117
118 // print initial state
119 out << "Time = 0.0:" << std::endl;
120 out << "Quanta positions:" << std::endl << p.x << std::endl;
121 out << "Quanta velocities:" << std::endl << p.v << std::endl;
122
123 // Advance particles in each time step according to:
124 // dx/dt = v
125 // dv/dt = -omega^2 * x
126 for (int it=0; it<numit; ++it)
127 {
128 p.x = p.x + dt * p.v;
129 p.v = p.v - dt * omega * omega * p.x;
130 out << "Time = " << (it+1)*dt << ":" << std::endl;
131 out << "Quanta positions:" << std::endl << p.x << std::endl;
132 out << "Quanta velocities:" << std::endl << p.v << std::endl;
133 }
134
135 // Finalize POOMA
136 Pooma::finalize();
137 return 0;
138 }
As discussed earlier, the program begins by creating a traits class that typedefs the names AttributeEngineTag_t and ParticleLayout_t (lines 11-21). An application-specific class called Quanta is then derived from Particles, without specifying the traits to be used (lines 25-48). This class declares two attributes, to store the particles' x coordinate and velocity. The body of its constructor (lines 40-41) adds these attributes to its attribute list, while passing the actual layout object specified by the application up to Particles.
Lines 54, 57 and 58 create some convenience typedefs for the engine and layout that the application will use. Lines 61-65 then define constants describing both the physical parameters to the problem (such as the oscillation frequency) and the computational parameters (the number of particles, the number of patches, etc.). In a real application, many of these values would be variables, rather than hard-wired constants.
After the POOMA library is initialized (line 71), an Inform object is created to manage output. (See the appendix on I/O for a description of this class.) An actual layout is then created (line 76), and used to create an actual set of particles (line 80). The particles themselves are created by the call to Particles::create() on line 83. The output on lines 84-89 shows that all particles are initially created in the zeroth patch.
The sync() call on line 92 redistributes particles across the available patches according to their x coordinates. As the output from lines 95-100 shows, this load-balances the particles as evenly as possible.
The particle positions are randomized on lines 108-116. (A loop is used here because random number generation has not yet been integrated into the expression evaluation machinery in this release of POOMA.) After some more output to show the particles' initial positions, the application finally enters the main timestep loop (lines 126-133). In each time step, particle positions and velocities are updated under the influence of a simple harmonic oscillator force, and then printed out. Once the specified number of timesteps has been executed, the library is shut down (line 136) and the application exits.
In addition to an AttributeList, each Particles object also stores a ParticleBCList of boundary conditions to be applied to the attributes. These are generalized boundary conditions in the sense that they can be applied not only to a particle position attribute, but to any sort of attribute or expression involving attributes. POOMA provides typical particle boundary conditions including periodicity, reflection, absorption, reversal (reflection of one attribute and negation of another), and kill (destroying a particle). Boundary conditions can be updated explicitly by calling Particles::applyBoundaryConditions(), or implicitly by calling Particles::sync() (which performs the same operations, along with several others).
Each boundary condition is assembled by first constructing an instance of the type of boundary condition desired, then invoking the addBoundaryCondition() member function of Particles with three parameters: the subject of the boundary condition (i.e., the attribute or expression to be checked against the range), its object (the attribute to be modified when the subject is outside the range), and the actual boundary condition object. The boundary condition is then applied each time the sync() function is invoked.
The subject and object of a boundary condition are usually the same, but this is not required. In one common case, the subject is an expression involving particle attributes, while the object is the Particles object itself. For example, an application's boundary condition might specify that particles are to be deleted if their kinetic energy goes above some limit. The subject would be the expression 0.5*m*v*v, and the object could be either one of the particle attributes (because deleting a particle from one attribute automatically deletes it from all the others) or the Particles object itself. The object cannot be the expression 0.5*m*v*v because that is a ConstArray and cannot be modified.
Another case involves the reversal boundary condition, which is used to make particles bounce off walls. Bouncing not only reflects the particle position back inside the wall, but also reverses the particle's velocity component in that direction. The reversal boundary condition therefore needs an additional object besides the original subject.
POOMA provides the pre-defined boundary condition classes listed in the table below.
| Class | Behavior |
| AbsorbBC<T>(T min, T max) | Keeps attributes within given limits min or max. If they cross the given boundaries, their values are changed to the given limiting value. |
| KillBC<T>(T min, T max) | If particles cross outside the given boundary, they are destroyed by putting their index in the deferred destroy list. |
| PeriodicBC<T>(T min, T max) | Keeps attributes within a given periodic domain. |
| ReflectBC<T>(T min, T max) | Reflects an attribute back if it crosses outside of the given boundary. |
| ReverseBC<T>(T min, T max) | Reverses (negates) the value of the object attribute if it crosses outside the given domain, and reflects the value of the subject attribute. |
As an example of how particle boundary conditions are used, consider a set of particles bouncing around in a box in three dimensions. examples/Particles/Bounce/Bounce.cpp shows how this can be implemented using POOMA for the case of perfectly elastic collisions. The code is:
001 #include "Pooma/Particles.h"
002 #include "Pooma/DynamicArrays.h"
003 #include "Tiny/Vector.h"
004 #include "Utilities/Inform.h"
005 #include <iostream>
006 #include <stdlib.h>
007
008
009 // Dimensionality of this problem
010 static const int PDim = 3;
011
012 // Particles subclass with position and velocity
013 template <class PT>
014 class Balls : public Particles<PT>
015 {
016 public:
017 // Typedefs
018 typedef Particles<PT> Base_t;
019 typedef typename Base_t::AttributeEngineTag_t AttributeEngineTag_t;
020 typedef typename Base_t::ParticleLayout_t ParticleLayout_t;
021 typedef double AxisType_t;
022 typedef Vector<PDim,AxisType_t> PointType_t;
023
024 // Constructor: set up layouts, register attributes
025 Balls(const ParticleLayout_t &pl)
026 : Particles<PT>(pl)
027 {
028 addAttribute(pos);
029 addAttribute(vel);
030 }
031
032 // Position and velocity attributes (as public members)
033 DynamicArray< PointType_t, AttributeEngineTag_t > pos;
034 DynamicArray< PointType_t, AttributeEngineTag_t > vel;
035 };
036
037 // Use canned traits class from CommonParticleTraits.h
038 // MPBrickUniform provides MultiPatch Brick Engine for
039 // particle attributes and UniformLayout for particle data.
040 typedef MPBrickUniform PTraits_t;
041
042 // Type of particle layout
043 typedef PTraits_t::ParticleLayout_t ParticleLayout_t;
044
045 // Type of actual particles
046 typedef Balls<PTraits_t> Particle_t;
047
048 // Number of particles in simulation
049 const int NumPart = 100;
050
051 // Number of timesteps in simulation
052 const int NumSteps = 100;
053
054 // Number of patches to distribute particles across.
055 // Typically one would use one patch per processor.
056 const int numPatches = 16;
057
058 // Main simulation routine
059 int main(int argc, char *argv[])
060 {
061 // Initialize POOMA and output stream
062 Pooma::initialize(argc, argv);
063 Inform out(argv[0]);
064
065 out << "Begin Bounce example code" << std::endl;
066 out << "-------------------------" << std::endl;
067
068 // Create a particle layout object for our use
069 ParticleLayout_t particleLayout(numPatches);
070
071 // Create the actual Particles object (but not the particles as yet)
072 Particle_t balls(particleLayout);
073
074 // Create some particles, recompute the global domain, and initialize
075 // the attributes randomly. The globalCreate call will create an equal
076 // number of particles on each patch. The particle positions are initialized
077 // within a 12 X 20 X 28 domain, and the velocity components are all
078 // in the range -4 to +4.
079 balls.globalCreate(NumPart);
080 srand(12345U);
081 Particle_t::PointType_t initPos, initVel;
082 typedef Particle_t::AxisType_t Coordinate_t;
083 Coordinate_t ranmax = static_cast<Coordinate_t>(RAND_MAX);
084 for (int i = 0; i < NumPart; ++i)
085 {
086 for (int d = 0; d < PDim; ++d)
087 {
088 initPos(d) = 4 * (2 * (d+1) + 1) * (rand() / ranmax);
089 initVel(d) = 4 * (2 * (rand() / ranmax) - 1);
090 }
091 balls.pos(i) = initPos;
092 balls.vel(i) = initVel;
093 }
094
095 // Display the particle positions and velocities.
096 out << "Timestep 0: " << std::endl;
097 out << "Ball positions: " << balls.pos << std::endl;
098 out << "Ball velocities: " << balls.vel << std::endl;
099
100 // Set up a reversal boundary condition, so that particles will
101 // bounce off the domain boundaries.
103 Particle_t::PointType_t lower, upper;
104 for (int d = 0; d < PDim; ++d)
105 {
106 lower(d) = 0.0;
107 upper(d) = (d+1) * 8.0 + 4.0;
108 }
109 ReverseBC<Particle_t::PointType_t> bounce(lower, upper);
110 balls.addBoundaryCondition(balls.pos, balls.vel, bounce);
111
112 // Advance simulation stepwise
113 for (int it=1; it <= NumSteps; ++it)
114 {
115 // Advance ball positions (timestep dt = 1)
116 balls.pos += balls.vel;
117
118 // Invoke boundary conditions
119 balls.applyBoundaryConditions();
120
121 // Print out the current particle data
122 out << "Timestep " << it << ": " << std::endl;
123 out << "Ball positions: " << balls.pos << std::endl;
124 out << "Ball velocities: " << balls.vel << std::endl;
125 }
126
127 // Shut down POOMA and exit
128 Pooma::finalize();
129 return 0;
130 }
After defining the dimension of the problem (line 10), this program defines a class Balls, which represents the set of particles (lines 13-35). Its two attributes represent the particles' positions and velocities (lines 33-34). Note how the type of engine used for evaluating these attributes is defined in terms of the types exported by the traits class with which Balls is instantiated (AttributeEngineTag_t, line 19), while the type used to represent the points is defined in terms of the dimension of the problem (line 22), rather than being made 1-, 2-, or 3-dimensional explicitly. This style of coding makes it much easier to change the simulation as the program evolves.
Rather than defining a particle traits class explicitly, as the oscillation example above did, this program uses one of the pre-defined traits class given in src/Particles/CommonParticleTraits.h. For the purposes of this example, a multipatch brick engine is used for particle attributes, and particle data is laid out uniformly. Once again, a typedef is used to create a symbolic name for this combination, so that the program can be updated by making a single change in a single location.
Lines 43-56 then define the types used in the simulation, and the constants that control the simulation's evolution. It would be possible to shorten this part of the program by combining some of these type definitions (as on line 43), but readability would suffer.
The main body of the program follows; as usual, it begins by initializing the POOMA library, and creating an output handler of type Inform (lines 62-63). Line 69 then creates a layout object describing the domain of the problem.
The particles object itself comes into being on line 72, although the actual particles aren't created until line 79. Recall that by default, globalCreate() (re-)numbers the particles by calling Particles' renumber() method. As discussed earlier, this could be prevented by passing false as a second parameter to globalCreate(), i.e., by calling globalCreate(N, false). Lines 80-93 then randomize the balls' initial positions and velocities.
Lines 103-110 are the most novel part of this simulation, as they create reflecting boundary conditions for the simulation, and add them to the balls object. Lines 103-108 define where particles bounce; again, this is done in a dimension-independent fashion in order to make code evolution as easy as possible. Line 109 turns upper and lower into a reversing boundary condition, which line 110 then adds to balls. The main simulation loop now consists of nothing more than advancing the balls in each time step, and calling applyBoundaryCondition() to enforce the boundary conditions.
Particles are a fundamental construct in physical calculations. POOMA's Particles class, and the classes that support it, allow programmers to create and manage sets of particles both efficiently and flexibly. While doing this is a multi-step process, the payoff as programs are extended and updated is considerable. The list below summarizes the most important aspects of Particles' interface.
Note: if the patchId given to destroy() or
deferredDestroy() is negative, the domain argument
must specify a global domain (i.e., global numbering). If the argument
is non-negative, the domain is interpreted as being local, i.e., the
index 0 refers to the first particle in that patch.
| [Prev] | [Home] | [Next] |