 # POOMA Tutorial 10Particles and Fields

## Introduction

The previous tutorials have described how POOMA represents fields and particles. This tutorial shows how the two can be combined to create complete simulations of complex physical systems. The first section describes how POOMA interpolates values when gathering and scattering field and particle data. This is followed by a look at the in's and out's of layout, and a medium-sized example that illustrates how these ideas fit together.

## Particle/Field Interpolation

POOMA's Particles class is designed to be used in conjunction with its Fields. Interpolators are the glue that bind these together, by specifying how to calculate field values at particle (or other) locations that don't happen to lie exactly on mesh points.

Interpolators are used to gather values to specific positions in a field's spatial domain from nearby field elements, or to scatter values from such positions into the field. The interpolation stencil describes how values are translated between field element locations and arbitrary points in space. An example of using this kind of interpolation is particle-in-cell (PIC) simulations, in which charged particles move through a discretized domain. The particle interactions are determined by scattering the particle charge density into a field, solving for the self-consistent electric field, and gathering that field back to the particle positions. The last example in this tutorial describes a simulation of this kind.

POOMA currently offers three types of interpolation stencils: nearest grid point (NGP), cloud-in-cell (CIC), and subtracted dipole scheme (SUDS). NGP is a zeroth-order interpolation that gathers from or scatters to the field element nearest the specified location. CIC is a first-order scheme that performs linear weighting among the 2D field elements nearest the point in D-dimensional space. SUDS is also first-order, but it uses just the nearest field element and its two neighbors along each dimension, so it is only a 7-point stencil in three dimensions. Other types of interpolation schemes can be added in a straightforward manner. Figure 1: Interpolation strategies. Black dots show particle positions, and open circles are the interpolation stencil points.

Interpolation is invoked by calling the global functions gather() and scatter(), both of which take four arguments:

1. the particle attribute to be gathered to or scattered from (usually a single DynamicArray, although one could scatter an expression involving DynamicArrays as well, since the evaluation of this expression just produces a temporary one-dimensional ConstArray);
2. the Field to be gathered from or scattered to;
3. the particle positions (normally a DynamicArray that is a member of a Particles-derived class); and
4. an interpolator tag object of type NGP, CIC or SUDS. These tag objects are defined in the header files InterpolatorNGP.h, InterpolatorCIC.h, and InterpolatorSUDS.h respectively.

An example of this is:

```gather(P.efd, Efield, P.pos, CIC());
```

where P is a Particles subclass object whose attributes are efd for storing the gathered electric field from the Field Efield and pos for the particle positions. The default constructor of the interpolator CIC is used to create a temporary instance of the class to pass to gather(), telling it which interpolation scheme to use.

The particle attribute and position arguments passed to gather() and scatter() should have the same layout, and the positions must refer to the geometry of the Field being used. The interpolator will compute the required interpolated values for the particles on each patch. These functions assume each particle is only interacting with field elements in the Field patch that exactly corresponds to the particle patch. Thus, applications must use the SpatialLayout particle layout strategy and make sure that the Field has enough guard layers to accommodate the interpolation stencil.

In addition to the basic gather() and scatter() functions, POOMA offers some variants that optimize other common operations. The first of these, scatterValue(), scatters a single value into a Field rather than a particle attribute with different values for each particle. Its first argument is a single value with a type that is compatible with the Field element type.

The other three optimized methods are gatherCache(), scatterCache(), and scatterValueCache(). Each of these has two overloaded variants, which allow applications to cache and reuse interpolated data, such as the nearest grid point for each particle and the distance from the particle's position to that grid point. The difference between the elements of each overloaded pair of methods is that one takes both a particle position attribute and a particle interpolator cache attribute among its arguments, while the other takes only the cache attribute. When the first of these is called, it caches position information in the provided cache attribute. When the second is called with that cache attribute as an argument, it re-uses that information. This can speed up computation considerably, but it is important to note that applications can only do this safely when the particle positions are guaranteed not to have changed since the last interpolation.

## Laying Out Particles and Fields

The use of particles and fields together in a single application brings up some issues regarding layout that do not arise when either is used on its own. There are two characteristics of Engines that must be considered in order to determine whether they can be used for attributes in Particles objects:

1. Can the engine use a layout that is "shared" among several engines of the same category, such that the size and layout of the engine is synchronized with the other engines using the layout? If this is the case, then creation, destruction, repartitioning, and other operations are done for all the shared engines. Particles require all their attributes to use a shared layout, so only engines that use a shared layout can be used for particle attributes. The only engines with this capability in this release of POOMA (i.e., the only engines that are usable in Particles attributes) are SharedBrick (using a SharedDomainLayout layout) and some specializations of MultiPatch.

MultiPatch can use several different types of layouts and single-block engines, and all MultiPatchs use a shared layout. However, only the MultiPatch<GridTag,*> types of MultiPatch engines are useful for Particles attributes, since only that engine type can have patches of varying size. Future releases of POOMA will add other layouts, such as TileLayout, that will also be useful for attributes.

Note that MultiPatch<UniformTag,Brick> can not be used for particle attributes, as it uses a UniformGridLayout. While that layout is "shared", UniformGridLayout is not useful for particles because it requires all patches to have the same size, and particle attribute patches change their size dynamically.

2. How many patches can the engine have? A SharedBrick can only have one patch, but a Grid-based MultiPatch can have several patches. Either one can be used in serial or in parallel, but their efficiency will differ. If individual particle attribute expressions will normally be run in parallel, the application should use a MultiPatch. Otherwise, it should use a SharedBrick.

Implicit in the discussion above is the fact that there are actually three different types of layout classes that an application programmer must keep in mind:

1. the layout for the particle attributes;
2. the layout for the Field given to the particle SpatialLayout (which is used to determine the layout of the space in which the particles move around); and
3. the actual SpatialLayout that connects the info about the Field layout to the Particles attribute layout.

The only thing that needs to match between the attribute and Field layouts is the number of patches, which must be the same. The engine type (and thus the layout type) of the attributes and of the field do not have to match. An application could therefore use a SharedBrick engine for particle attributes, and a MultiPatch<UniformTag,Brick> for the Field, as long as the MultiPatch engine uses just one patch (since SharedBrick can only have one patch).

Note once again that in the simple case of a UniformLayout, applications do not need to worry about the Field layout type, only the particle attributes' layout (which still needs to be shared) and the particle layout (in this case, UniformLayout). This commonly arises during the prototyping (i.e., pre-parallel) stages of application development.

## Example: Particle-in-Cell Simulation

Our third and final example of this important class is a particle-in-cell program, which simulates the motion of charged particles in a static sinusoidal electrical field in two dimensions. This example brings together the Field classes of the preceding tutorials with this tutorial's Particles class.

Because this example is longer than the others in these tutorials, it will be described in sections. For a unified listing of the source code, please see the file examples/Particles/PIC2d/PIC2d.cpp in the distribution.

The first step is to include all of the usual header files:

```001  #include "Pooma/Particles.h"
002  #include "Pooma/DynamicArrays.h"
003  #include "Pooma/Fields.h"
004  #include "Utilities/Inform.h"
005  #include <iostream>
006  #include <stdlib.h>
007  #include <math.h>
```

Once this has been done, the application can define a traits class for the Particles object it is going to create. As always, this contains typedefs for AttributeEngineTag_t and ParticleLayout_t. The traits class for this example also includes an application-specific typedef called InterpolatorTag_t, for reasons discussed below.

```008  template <class EngineTag, class Centering, class MeshType, class FL,
009            class InterpolatorTag>
010  struct PTraits
011  {
012    // The type of engine to use in the attributes
013    typedef EngineTag AttributeEngineTag_t;
014
015    // The type of particle layout to use
016    typedef SpatialLayout<DiscreteGeometry<Centering,MeshType>,FL>
017      ParticleLayout_t;
018
019    // The type of interpolator to use
020    typedef InterpolatorTag InterpolatorTag_t;
021  };
```

The interpolator tag type is included in the traits class because an application might want the Particles-derived to provide the type of interpolator to use. One example of this is the case in which a gather() or scatter() call occurs in a subroutine which is passed an object of a Particles-derived type. This subroutine could extract the desired interpolator type from that object using:

```// Particles-derived type Particles_t already defined
typedef typename Particles_t::InterpolatorTag_t InterpolatorTag_t;
```

In this short example, this is not really necessary because InterpolatorTag_t is being defined and then used within the same file scope. Nevertheless, this illustrates a situation in which the user might want to add new traits to their PTraits class beyond the required traits AttributeEngineTag_t and ParticleLayout_t.

We can now also define the class which will represent the charged particles in the simulation. As in other examples, this is derived from Particles, and templated on a traits class so that such things as its layout and evaluation engine can be quickly, easily, and reliably changed. This class has three intrinsic properties: the particles' positions R, their velocities V, and their charge/mass ratios qm. The class also has a fourth property called E, which is used to record the electrical field at each particle's position in order to calculate forces. This calculation will be discussed in greater detail below.

```024  template <class PT>
025  class ChargedParticles : public Particles<PT>
026  {
027  public:
028    // Typedefs
029    typedef Particles<PT>                          Base_t;
030    typedef typename Base_t::AttributeEngineTag_t  AttributeEngineTag_t;
031    typedef typename Base_t::ParticleLayout_t      ParticleLayout_t;
032    typedef typename ParticleLayout_t::AxisType_t  AxisType_t;
033    typedef typename ParticleLayout_t::PointType_t PointType_t;
034    typedef typename PT::InterpolatorTag_t         InterpolatorTag_t;
035
036    // Dimensionality
037    static const int dimensions = ParticleLayout_t::dimensions;
038
039    // Constructor: set up layouts, register attributes
040    ChargedParticles(const ParticleLayout_t &pl)
041    : Particles<PT>(pl)
042    {
047    }
048
049    // Position and velocity attributes (as public members)
050    DynamicArray<PointType_t,AttributeEngineTag_t> R;
051    DynamicArray<PointType_t,AttributeEngineTag_t> V;
052    DynamicArray<PointType_t,AttributeEngineTag_t> E;
053    DynamicArray<double,     AttributeEngineTag_t> qm;
054  };
```

With the two classes that the simulation relies upon defined, the program next defines the dependent types, constants, and other values that the application needs. These include the dimensionality of the problem (which can easily be increased to 3), the type of mesh on which the calculations are done, the mesh's size, and so on:

```058  // Dimensionality of this problem
059  static const int PDim = 2;
060
061  // Engine tag type for attributes
062  typedef MultiPatch<GridTag,Brick> AttrEngineTag_t;
063
064  // Mesh type
065  typedef UniformRectilinearMesh<PDim,Cartesian<PDim>,double> Mesh_t;
066
067  // Centering of Field elements on mesh
068  typedef Cell Centering_t;
069
070  // Geometry type for Fields
071  typedef DiscreteGeometry<Centering_t,Mesh_t> Geometry_t;
072
073  // Field types
074  typedef Field< Geometry_t, double,
075                 MultiPatch<UniformTag,Brick> > DField_t;
076  typedef Field< Geometry_t, Vector<PDim,double>,
077                 MultiPatch<UniformTag,Brick> > VecField_t;
078
079  // Field layout type, derived from Engine type
080  typedef DField_t::Engine_t Engine_t;
081  typedef Engine_t::Layout_t FLayout_t;
082
083  // Type of interpolator
084  typedef NGP InterpolatorTag_t;
085
086  // Particle traits class
087  typedef PTraits<AttrEngineTag_t,Centering_t,Mesh_t,FLayout_t,
088                  InterpolatorTag_t> PTraits_t;
089
090  // Type of particle layout
091  typedef PTraits_t::ParticleLayout_t PLayout_t;
092
093  // Type of actual particles
094  typedef ChargedParticles<PTraits_t> Particles_t;
095
096  // Grid sizes
097  const int nx = 200, ny = 200;
098
099  // Number of particles in simulation
100  const int NumPart = 400;
101
102  // Number of timesteps in simulation
103  const int NumSteps = 20;
104
105  // The value of pi (some compilers don't define M_PI)
106  const double pi = acos(-1.0);
107
108  // Maximum value for particle q/m ratio
109  const double qmmax = 1.0;
110
111  // Timestep
112  const double dt = 1.0;
```

The preparations above might seem overly elaborate, but the payoff comes when the main simulation routine is written. After the usual initialization call, and the creation of an Inform object to handle output, the program defines one geometry object to represent the problem domain, and another that includes a guard layer:

```115  int main(int argc, char *argv[])
116  {
117    // Initialize POOMA and output stream.
118    Pooma::initialize(argc, argv);
119    Inform out(argv);
120
121    out << "Begin PIC2d example code" << std::endl;
122    out << "-------------------------" << std::endl;
123
124    // Create mesh and geometry objects for cell-centered fields.
125    Interval<PDim> meshDomain(nx+1,ny+1);
126    Mesh_t mesh(meshDomain);
127    Geometry_t geometry(mesh);
128
129    // Create a second geometry object that includes a guard layer.
130    GuardLayers<PDim> gl(1);
131    Geometry_t geometryGL(mesh,gl);
```

The program then creates a pair of Field objects. The first, phi, is a field of double values, and records the electrostatic potential at points in the mesh. The second, EFD, is a field of two-dimensional Vectors, and records the electric field at each mesh point. The types used in these definitions were declared on lines 74-77 above. Note how these definitions are made in terms of other defined types, such as Geometry_t, rather than directly in terms of basic types. This allows the application to be modified quickly and reliably with minimal changes to the code.

```133    // Create field layout objects for our electrostatic potential
134    // and our electric field.  Decomposition is 4 x 4.
135    Loc<PDim> blocks(4,4);
136    FLayout_t flayout(geometry.physicalDomain(),blocks);
137    FLayout_t flayoutGL(geometryGL.physicalDomain(),blocks,gl);
138
139    // Create and initialize electrostatic potential and electric field.
140    DField_t phi(geometryGL,flayoutGL);
141    VecField_t EFD(geometry,flayout);
```

The application now adds periodic boundary conditions to the electrostatic field phi, so that particles will not see sharp changes in potential at the edges of the simulation domain. The values of phi and EFD are then set: phi is defined explicitly, while EFD records the gradient of phi.

```144    // potential phi = phi0 * sin(2*pi*x/Lx) * cos(4*pi*y/Ly)
145    // Note that phi is a periodic Field
146    // Electric field EFD = -grad(phi);
148    double phi0 = 0.01 * static_cast<double>(nx);
149    phi = phi0 * sin(2.0*pi*phi.x().comp(0)/nx)
150               * cos(4.0*pi*phi.x().comp(1)/ny);
```

With the fields in place, the application creates the particles whose motions are to be simulated, and adds periodic boundary conditions to this object as well. The globalCreate() call creates the same number of particles on each processor.

```153    // Create a particle layout object for our use
154    PLayout_t layout(geometry,flayout);
155
156    // Create a Particles object and set periodic boundary conditions
157    Particles_t P(layout);
158    Particles_t::PointType_t lower(0.0,0.0), upper(nx,ny);
159    PeriodicBC<Particles_t::PointType_t> bc(lower,upper);
161
162    // Create an equal number of particles on each processor
163    // and recompute global domain.
164    P.globalCreate(NumPart);
```

Note that the definitions of lower and upper could be made dimension-independent by defining them with a loop. If ng is an array of ints of length PDim, then this loop is:

```Particles_t::PointType_t lower, upper;
for (int d=0; d<PDim; ++d)
{
lower(d) = 0;
upper(d) = ng[d];
}
```

The application then randomizes the particles' positions and charge/mass ratios using a sequential loop (since parallel random number generation is not yet in POOMA). Once this has finished, the method swap() is called to redistribute the particles based on their positions, i.e., to move each particle to its home processor. The initial positions, velocities, and charge/mass ratios of the particles are then printed out.

```166    // Random initialization for particle positions in nx by ny domain
167    // Zero initialization for particle velocities
168    // Random intialization for charge-to-mass ratio from -qmmax to qmmax
169    P.V = Particles_t::PointType_t(0.0);
170    srand(12345U);
171    Particles_t::PointType_t initPos;
172    typedef Particle_t::AxisType_t Coordinate_t;
173    Coordinate_t ranmax = static_cast<Coordinate_t>(RAND_MAX);
174    double ranmaxd = static_cast<double>(RAND_MAX);
175    for (int i = 0; i < NumPart; ++i)
176    {
177      initPos(0) = nx * (rand() / ranmax);
178      initPos(1) = ny * (rand() / ranmax);
179      P.R(i) = initPos;
180      P.qm(i) = qmmax * (2 * (rand() / ranmaxd) - 1);
181    }
182
183    // Redistribute particle data based on spatial layout
184    P.swap(P.R);
185
186    out << "PIC2d setup complete." << std::endl;
187    out << "---------------------" << std::endl;
188
189    // Display the initial particle positions, velocities and qm values.
190    out << "Initial particle data:" << std::endl;
191    out << "Particle positions: "  << P.R << std::endl;
192    out << "Particle velocities: " << P.V << std::endl;
193    out << "Particle charge-to-mass ratios: " << P.qm << std::endl;
```

The application is finally able to enter its main timestep loop. In each time step, the particles' positions are updated, and then sync() is called to invoke boundary conditions, swap particles, and then renumber. A call is then made to gather() (line 208) to determine the field at each particle's location. As discussed earlier, this function uses the interpolator to determine values that lie off mesh points. Once the field strength is known, the particles' velocities can be updated:

```195    // Begin main timestep loop
196    for (int it=1; it <= NumSteps; ++it)
197    {
199      out << "Advance particle positions ..." << std::endl;
200      P.R = P.R + dt * P.V;
201
202      // Invoke boundary conditions and update particle distribution
203      out << "Synchronize particles ..." << std::endl;
204      P.sync(P.R);
205
206      // Gather the E field to the particle positions
207      out << "Gather E field ..." << std::endl;
208      gather( P.E, EFD, P.R, Particles_t::InterpolatorTag_t() );
209
210      // Advance the particle velocities
211      out << "Advance particle velocities ..." << std::endl;
212      P.V = P.V + dt * P.qm * P.E;
213    }
```

Finally, the state of the particles at the end of the simulation is printed out, and the simulation is closed down:

```215    // Display the final particle positions, velocities and qm values.
216    out << "PIC2d timestep loop complete!" << std::endl;
217    out << "-----------------------------" << std::endl;
218    out << "Final particle data:" << std::endl;
219    out << "Particle positions: "  << P.R << std::endl;
220    out << "Particle velocities: " << P.V << std::endl;
221    out << "Particle charge-to-mass ratios: " << P.qm << std::endl;
222
223    // Shut down POOMA and exit
224    out << "End PIC2d example code." << std::endl;
225    out << "-----------------------" << std::endl;
226    Pooma::finalize();
227    return 0;
```

## Summary

This tutorial has shown how POOMA's Field and Particles classes can be combined to create complete physical simulations. While more setup code is required than with Fortran-77 or C, the payoff is high-performance programs that are more flexible and easier to maintain.