[PATCH] Convert ParticlesDoc.txt to html
Richard Guenther
rguenth at tat.physik.uni-tuebingen.de
Mon Aug 23 11:17:50 UTC 2004
As subject says. Also adds common header to Layout.html.
Ok?
Richard.
2004Aug23 Richard Guenther <richard.guenther at uni-tuebingen.de>
* docs/Layout.html: adjust background color, add head image.
docs/index.html: refer to ParticlesDoc.html.
docs/ParticlesDoc.html: new.
docs/ParticlesDoc.txt: remove.
-------------- next part --------------
Index: index.html
===================================================================
RCS file: /home/pooma/Repository/r2/docs/index.html,v
retrieving revision 1.3
diff -u -u -r1.3 index.html
--- index.html 20 Aug 2004 20:14:19 -0000 1.3
+++ index.html 23 Aug 2004 11:14:54 -0000
@@ -23,7 +23,7 @@
<blockquote>
<h4><a href="parallelism.html">Parallelism Models: Messaging and Threads</a></h4>
<h4><a href="Layout.html">Layouts</a></h4>
- <h4><a href="ParticlesDoc.txt">New description of Particles</a></h4>
+ <h4><a href="ParticlesDoc.html">New description of Particles</a></h4>
<h4><a href="tut-11.html">Text I/O</a></h4>
<h4><a href="tut-12.html">Object I/O</a></h4>
<h4><a href="tut-04.html#tiny">New <tt>Tensor</tt> functionality</a></h4>
--- /dev/null Tue May 18 17:20:27 2004
+++ ParticlesDoc.html Mon Aug 23 13:13:27 2004
@@ -0,0 +1,1520 @@
+<!DOCTYPE HTML PUBLIC "-//W3C//DTD HTML 4.0 Transitional//EN">
+<html>
+<head>
+ <meta http-equiv="Content-Type" content="text/html; charset=iso-8859-1">
+ <meta name="GENERATOR" content="Mozilla/4.72 [en] (X11; U; Linux 2.2.14 i686) [Netscape]">
+ <title>Layout and related classes</title>
+</head>
+<body background="back.gif" LINK="#505062" ALINK="#505062" VLINK="#be7c18">
+
+<CENTER><IMG SRC="banner.gif" ALT="POOMA banner" WIDTH=550 HEIGHT=100
+ALIGN=bottom></CENTER>
+
+
+<h1>POOMA Particles Documentation</h1>
+
+
+<h2>Introduction</h2>
+
+<p>
+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.
+
+<p>
+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 or momentum. Particles
+typically exist in a spatial domain, and they may interact directly
+with one another or with field quantities defined on that domain.
+
+<p>
+This document 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. Later on, we will show
+how particles and fields can interact in a simulation code.
+
+
+<h2>Overview</h2>
+
+<p>
+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 a set of POOMA DynamicArray objects), 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.
+
+<p>
+Each Particles object keeps a list of pointers to its elements'
+attributes. When an application wants to add or delete particles, 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 global particle count and numbering,
+update the data distribution across patches, and apply the particle
+boundary conditions.
+
+<p>
+Applications define a specific 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 derived class should call the method
+Particles::addAttribute() to register each attribute and add
+it to the list maintained by Particles. In this way, the Particles
+class can be extended by the application to accommodate any sort of
+particle description.
+
+<p>
+The distribution of particle data stored in DynamicArray objects
+is directed by a particle layout class. 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 neighboring particles.
+
+
+<h2>Particle Attributes</h2>
+
+<p>
+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
+(ShiftUp). At the same time, DynamicArray objects can be used in
+data-parallel expressions in the same way as ordinary Array objects,
+so that the application can update particle attributes such as
+position and velocity using either a single statement or a loop over
+individual particles.
+
+<p>
+At first glance, it might seem more sensible to have applications
+define a struct 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 makes typical
+computational loops such as
+
+<pre>
+for (int i=0; i<n; ++i)
+{
+ x[i] += dt * vx[i];
+ y[i] += dt * vy[i];
+}
+</pre>
+
+<p>
+run more quickly, as it makes much better use of the data cache.
+
+
+<h2>Particle Layout</h2>
+
+<p>
+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 the method Particles::sync().
+
+
+<h2>Derivation of Particles Subclass</h2>
+
+<p>
+In general, creating a new Particles subclass is a three-step process.
+The first step is to declare a traits class with typedef's specifying
+the type of engine the particle attributes will use and the type of
+particle layout. An example of such a traits class is the following:
+
+<pre>
+struct MyParticleTraits
+{
+ typedef MultiPatch<DynamicTag,Dynamic> AttributeEngineTag_t;
+ typedef UniformLayout ParticleLayout_t;
+};
+</pre>
+
+<p>
+This traits class will be used to specialize the Particles class
+template when an application-specific subclass representing a
+concrete set of particles is derived from it. Particles uses public
+typedef's to give sensible names to these traits parameters, so that
+the derived subclass can access them (as shown below). The traits
+approach is used here to provide flexibility in the Particles design
+for future extensions. In addition to specifying the attribute engine
+and particle layout types, this traits class could also set some
+application-specific parameters.
+
+<p>
+Currently, there is a fairly limited set of valid choices for attribute
+engine type and particle layout strategy. Because we require that
+the particle attributes share a common layout and remain synchronized,
+we must use a MultiPatch engine with a DynamicLayout. The patch
+engines inside the MultiPatch engine must have dynamic capabilities.
+They can be either of type Dynamic or, when running across multiple
+contexts, of type Remote<Dynamic>. As for the particle layout, POOMA
+currently provides only two possible strategies: UniformLayout, which
+just tries to keep a similar number of particles in each patch; and
+SpatialLayout, which organizes the particles into patches based upon
+their current spatial position. For the user'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
+Particles/CommonParticleTraits.h. These define all the combinations
+of multi-patch dynamic and remote dynamic engines with both uniform
+and spatial layouts. Ordinarily, the user can simply choose one of
+these pre-defined traits classes for their Particles subclass.
+
+<p>
+The second step is to actually 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, typedef's may be provided for the instantiated
+parent class and for its layout type. The constructor for the user's
+subclass then usually takes a concrete layout object of the type specified
+in the typedef above as a constructor argument:
+
+<pre>
+template <class PT>
+class MyParticles : public Particles<PT>
+{
+public:
+ // instantiated type of base class
+ typedef Particles<PT> Base_t;
+
+ // type of particle layout (from traits class via base class)
+ typedef typename Base_t::ParticleLayout_t ParticleLayout_t;
+
+ // type of attribute engine tag (from traits class via base 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 passes particle layout to base class
+ MyParticles(const ParticleLayout_t &layout)
+ : Particles<PT>(layout)
+ {
+ // register attributes with base class
+ addAttribute(charge);
+ addAttribute(mass);
+ addAttribute(count);
+ }
+};
+</pre>
+
+<p>
+Note that the attribute elements in this example have different
+element types; i.e., charge and mass are of type double, while
+count is of type int. Attribute elements may in general have any
+type, including any user-defined type.
+
+<p>
+Finally, the application-specific class MyParticles is instantiated
+with a traits class such as MyParticleTraits to create an actual
+set of particles. A particle layout is declared first, and it is
+passed as a constructor argument to the instance of the user's 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. Here is an example of creating a MyParticles object:
+
+<pre>
+const int numPatches = 10;
+MyParticleTraits::ParticleLayout_t layout(numPatches);
+MyParticles<MyParticleTraits> particles(layout);
+</pre>
+
+<p>
+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 a new class 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
+subclasses.
+
+
+<h2>Synchronization and Related Issues</h2>
+
+<p>
+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.
+
+<p>
+Particles::sync() is a member function template taking a single argument.
+This argument must be one of the particle attributes (i.e., a DynamicArray).
+SpatialLayout assumes that the attribute given to sync() is the particle
+positions, and it uses this 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-position
+attribute, such as temperature or pressure, to SpatialLayout via the
+sync() method.
+
+<p>
+UniformLayout, which divides particles as evenly as possible between patches,
+without regard to 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.
+
+<p>
+The use of a parameter in Particles::sync() is one important difference
+between the implementation of particles in POOMA I and POOMA II. In the old
+design, all Particles classes came with a pre-defined attribute named R that
+was the particles' position, and synchronization always referred to this
+attribute. The new scheme allows applications to change the attribute that
+is used to represent the position; e.g., to switch back and forth in a time
+advance algorithm between a "current" position attribute and a "new" position
+attribute. It also allows particles to be weighted according to some attribute,
+so that the distribution scheme load-balances by weight.
+
+<p>
+Of course, before particle data can be distributed, the particles themselves
+must be created. Particles provides two methods for doing this. The first,
+globalCreate(N,renum), creates a specified number of particles N, spread
+as evenly as possible across all patches. The particles are normally renumbered
+after the creation operation, although this can be overridden by passing the
+second parameter (renum) with a value of "false". POOMA automatically uses a
+numbering scheme in which the particles are ordered by patch number and labeled
+consecutively within a patch. For instance, if patch 0 has 6 particles and
+patch 1 has 4 particles, then the particles on patch 0 are labeled 0 through 5,
+and the particles on patch 1 are labeled 6 through 9.
+
+<p>
+Particles::create(N,patchID,renum), on the other hand, creates a specified
+number of particles N on each local context, and adds them to a specific
+patch (or to the last local patch if none is specified). Once again, the
+particles are renumbered after this operation unless renum is false. Used in
+conjunction with the Pooma::context() method, this create() method can be
+utilized to allocate a specific number of particles on each context and in
+each local patch within a context. If a program contains a series of calls
+to the create() method, the user may wish set renum to false to avoid
+renumbering particles until all of the particle creation tasks have been
+completed.
+
+<p>
+After particles have been created (or destroyed), they should be renumbered
+to ensure that each has a unique ID and that the global domain of the
+particle attributes is consistent. This is critical to the proper behavior
+of data-parallel expressions involving attributes. The Particles::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 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 across patches. Note that calls to
+globally synchronizing functions such as renumber() or sync() must be done
+on all contexts in a SPMD fashion.
+
+<p>
+If a program does not (implicitly or explicitly) call renumber() after
+creating or destroying particles, the global domain for the particle
+attributes 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 (incorrect) global domain, in which case the
+layout object for the particle attributes will suffer a run-time assertion
+failure. It is the user's responsibility to ensure that the particle
+attributes are properly synchronized prior to any data-parallel expression.
+
+<p>
+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). This function deletes
+particles from the local patch indicated by patchID, and domain is
+assumed to refer to a subset of the local domain in that patch. If
+patchID is not specified, then domain is assumed to refer to a subset
+of the global domain of the particle attributes, in which case the
+function may delete particles from multiple patches. The renum argument
+indicates whether to renumber the particles after the destroy command
+is performed, and it is true by default.
+
+<p>
+The second destruction method is Particles::deferredDestroy(domain,patchID).
+The meanings of the two arguments are the same as in the destroy() method.
+This is new in POOMA II, and it only does deferred destruction; i.e., 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(). This deferredDestroy() method can be useful when
+there are many separate destroy requests, because it lumps them all together
+and amortizes the expense of having to move particle data around and shrink
+the particle attributes. The performDestroy() method, which causes
+the cached destruction requests to be executed, always performs renumbering.
+The performDestroy() method can be called explicitly by the user in order
+to process and flush any cached destroy requests or implicitly by calling
+the sync() method.
+
+<p>
+All particle destroys are implemented using one of two possible methods:
+BackFill or ShiftUp. With the BackFill method, the "holes" that are left
+behind when particles are deleted are filled with particle data from the
+end of the list for the given patch. The ShiftUp method, on the other hand,
+slides all of the remaining particles up the list in order to fill holes.
+Thus, the ShiftUp destroy method is plainly slower, but it preserves the
+existing order of the particles within a given patch. The user can select
+the preferred destroy method using the setDestroyMethod() function.
+
+
+<h3>Example: Simple Harmonic Oscillator</h3>
+
+<p>
+The example for this tutorial simulates the motion of particles under the
+influence of a simple one-dimensional harmonic oscillator potential. The
+code, a version of which is included in the POOMA II release in the
+examples/Particles/Oscillation directory, is as follows:
+
+<pre>
+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 typedef's 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 (these could be declared
+045 // private and accessed with public methods, but this gives nice syntax)
+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 Dynamic engines, and a DynamicTag, which allows
+052 // the patches to change sizes during the application. This is important
+053 // since we may change the number of particles in each patch.
+054 typedef MultiPatch<DynamicTag,Dynamic> 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 particle 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 parallel 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=1; it<=nIter; ++it)
+127 {
+128 p.x = p.x + dt * p.v;
+129 p.v = p.v - dt * omega * omega * p.x;
+130 out << "Time = " << it*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 }
+</pre>
+
+<p>
+As discussed earlier, the program begins by creating a traits class that
+provides typedef's for 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 the attribute list, while passing the actual particle layout object
+specified by the application up to Particles.
+
+<p>
+Lines 54, 57 and 58 create some convenience typedef's 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.
+
+<p>
+After the POOMA library is initialized (line 71), an Inform object is
+created to manage output. An actual layout is then created (line 76), and
+is used to create a 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 patch 0.
+
+<p>
+The sync() call on line 92 redistributes particles across the available
+patches, using the x coordinate as a template for the current particle
+distribution. As the output from lines 95-100 shows, this distributes the
+particles across patches as evenly as possible.
+
+<p>
+The particle positions are randomized on lines 108-116. (A loop is used
+here rather than a data-parallel expression because parallel 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 and velocities, 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.
+
+
+<h2>Boundary Conditions</h2>
+
+<p>
+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().
+
+<p>
+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
+a specified 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 applyBoundaryConditions() is invoked.
+
+<p>
+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 an expression like 0.5*P.m*P.v*P.v, and the object would be P.
+The object cannot be the expression 0.5*P.m*P.v*P.v because an expression
+contains no actual data and thus cannot be modified.
+
+<p>
+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.
+
+<p>
+POOMA provides the pre-defined boundary condition classes listed in the table below.
+
+<table>
+<tr>
+<td>Class <td>Behavior
+<tr>
+<td>AbsorbBC <td>If attribute is outside specified lower or upper bounds, it is
+ reset to the boundary value.
+<tr>
+<td>KillBC <td>If particles go outside the given bounds, they are destroyed by
+ putting their index into the deferred destroy list.
+<tr>
+<td>PeriodicBC <td>Keeps attributes within a given periodic domain.
+<tr>
+<td>ReflectBC <td>If attribute exceeds a given boundary, its value is reflected
+ back inside that boundary.
+<tr>
+<td>ReverseBC <td>Reflects the value of the subject attribute if it crosses outside
+ the given domain, and reverses (negates) the value of the object
+ attribute.
+</table>
+
+
+<h3>Example: Elastic Collision</h3>
+
+<p>
+As an example of how particle boundary conditions are used, consider a set of
+particles bouncing around in a box in three dimensions. The sample code in
+file examples/Particles/Bounce/Bounce.cpp shows how this can be implemented
+using POOMA for the case of perfectly elastic collisions.
+
+<pre>
+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 // MPDynamicUniform provides MultiPatch Dynamic Engine for
+039 // particle attributes and UniformLayout for particle data.
+040 typedef MPDynamicUniform 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 Particles subclass object
+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 // Simulation timestep loop
+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 }
+</pre>
+
+<p>
+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 storing these attributes is defined in
+terms of the types exported by the traits class with which Balls is instantiated
+(AttributeEngineTag_t, line 19). Meanwhile the type used to represent the points
+is defined in terms of the dimension of the problem (line 22), rather than being
+made dimension-specific. This style of coding makes it much easier to change
+the simulation as the program evolves.
+
+<p>
+Rather than defining a particle traits class explicitly, as the oscillation
+example above did, this program uses one of the pre-defined traits classes
+given in src/Particles/CommonParticleTraits.h. For the purposes of this
+example, a multipatch dynamic 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 one location.
+
+<p>
+Lines 43-56 define the types used in the simulation, and the constants that
+control the simulation's evolution. 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.
+
+<p>
+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()
+renumbers the particles by calling the Particles::renumber() method. Lines 80-93
+then randomize the balls' initial positions and velocities.
+
+<p>
+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
+Particles object. Lines 103-108 defines where particles bounce; again, this
+is done in a dimension-independent fashion in order to make code evolution as
+easy as possible. Line 104 turns lower and upper into a reversing boundary
+condition, which line 105 then adds to balls. The main simulation loop now
+consists of nothing more than advancing the balls in each time step, and
+calling sync() to enforce the boundary conditions.
+
+
+<h2>Summary of Particles Interface</h2>
+
+<p>
+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 the Particles interface.
+
+<ul>
+<li>Particles<PT>(layout):
+Construct the Particles object with the given particle layout. This
+constructor will normally be called by the constructor of the user-defined
+Particles subclass.
+
+<li>initialize(layout):
+Initialize the Particles object with the given particle layout. This is
+used if the Particles object was created with the default constructor.
+
+<li>size():
+Return the current total number of particles since the last renumbering.
+
+<li>domain():
+Return the one-dimensional domain of the particle attributes, represented
+as the Interval<1> object storing the interval 0, 1, 2, ... size()-1.
+
+<li>attributes():
+Return the number of registered attributes.
+
+<li>addAttribute(attrib):
+Add the given attribute (which should be a DynamicArray of the proper
+engine type) to the attribute list stored by Particles.
+
+<li>removeAttribute(attrib):
+Remove the given attribute from the attribute list.
+
+<li>sync(posattrib):
+Apply boundary conditions, carry out cached destroys, swap particles
+across patches as needed, and renumber the particles. If relevant, the
+posattrib is used by the particle layout as a particle position attribute.
+
+<li>swap(posattrib):
+Move particle data between patches as specified by the particle layout
+strategy (uniform or spatial) and renumber the particles. This is more
+efficient than sync() if the user knows that there are no boundary
+conditions or cached destroy requests to carry out.
+
+<li>applyBoundaryConditions():
+Apply the boundary conditions to the current particle attributes, without
+renumbering or destroying any particles.
+
+<li>performDestroy():
+Destroy any particles that were specified in previous deferredDestroy()
+requests. (Note that these requests may be generated by a KillBC.)
+
+<li>renumber():
+Recalculate the per-patch and total domain of the system by inspecting
+the Particles attribute layout and resequencing the particles.
+
+<li>create(N, patchID, renum):
+Create N particles in the specified patch and optionally renumber. If
+the patchID and renum arguments are omitted, this creates the particles
+in the last patch, so as not to disturb the numbering of existing particles.
+
+<li>globalCreate(N, renum):
+Create N particles with a roughly equal number of particles created in
+each patch and optionally renumber.
+
+<li>setDestroyMethod(method):
+Set the preferred destroy method to be either BackFill (default) or ShiftUp.
+
+<li>destroy(domain, patchID, renum):
+Immediately destroy particles in the specified local domain within the
+specified patch and optionally renumber. The domain may be an Interval<1>
+or IndirectionList<int> of particle index numbers. If the patchID is
+omitted, the domain is assumed to be a global domain that may stretch
+across multiple patches.
+
+<li>deferredDestroy(domain, patchID):
+Put the indices of the particles in the given domain in the deferred
+destroy list of the Particles object, so that they will be destroyed by
+the next call to performDestroy(). If patchID is omitted, the domain
+is assumed to be a global domain.
+
+<li>addBoundaryCondition(Subj,Obj,BC) and addBoundaryCondition(Subj,BC):
+Add a new boundary condition BC that depends on the subject Subj and
+affects the object Obj (if different from Subj).
+
+<li>removeBoundaryCondition(i) and removeBoundaryConditions():
+Delete the ith boundary condition, or all boundary conditions.
+</ul>
+
+
+<h2>Particles and Fields</h2>
+
+<p>
+The previous sections have described how POOMA represents a set of
+particles and allows the user to perform typical operations in a
+particle simulation. The remainder of this document shows how POOMA
+Particles and Fields 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
+data layout, and a medium-sized example that illustrates how these
+ideas all fit together. (Note: The current implementation of POOMA
+Particles allows interaction with the original version of the Field
+abstraction created for POOMA II. Particles have not yet been
+modified to work with the new experimental design for POOMA Fields
+that is implemented in the src/Field directory. Thus, all the
+discussion here of POOMA Fields refers to the original implementation.)
+
+
+<h2>Particle/Field Interpolation</h2>
+
+<p>
+POOMA's Particles class is designed to be used in conjunction with the
+Field class. Interpolators are the glue that bind these two together,
+by specifying how to calculate field values at particle (or other)
+locations that lie in between the locations of Field elements. These
+interpolators can also be used to go in the opposite direction,
+acumulating contributions from particles at arbitrary locations into
+the elements of a Field.
+
+<p>
+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 to compute forces on the particles.
+The last code example in this document describes a simulation of this kind.
+
+<p>
+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
+2^D 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.
+
+<p>
+Interpolation is invoked by calling the global functions gather() and
+scatter(), both of which take four arguments:
+
+<ol>
+<li> the particle attribute to be gathered to or scattered from
+(usually a single DynamicArray, although one could scatter an
+expression involving DynamicArray objects as well);
+
+<li> the Field to be gathered from or scattered to;
+
+<li> the particle positions (normally a DynamicArray that is
+a member of a Particles subclass); and
+
+<li> an interpolator tag object of type NGP, CIC or SUDS, indicating
+which interpolation stencil to use.
+</ol>
+
+<p>
+An example of using the gather() function is:
+
+<pre>
+gather(P.efd, Efield, P.pos, CIC());
+</pre>
+
+<p>
+where P is a Particles subclass object whose attributes include efd
+for storing the gathered electric field from the Field Efield and
+pos for the particle positions. The default constructor of the
+interpolator tag CIC is used to create a temporary instance of the
+class to pass to gather(), telling it which interpolation scheme to use.
+
+<p>
+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's
+current 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.
+
+<p>
+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.
+
+<p>
+The other three optimized methods are gatherCache(), scatterCache(), and
+scatterValueCache(). Each of these methods 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 version 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.
+
+
+<h2>Data Layout for Particles and Fields</h2>
+
+<p>
+The use of particles and fields together in a single application brings
+up some issues regarding data layout that do not arise when either is
+used on its own. There are two characteristics of Engine objects that
+must be considered in order to determine whether they can be used for
+attributes in Particles objects:
+
+<ol>
+<li>
+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?
+
+<p>
+If this is the case, then creation, destruction, repartitioning, and
+other operations are done for all the engines sharing the layout.
+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 engine type with this capability in this release of POOMA
+(i.e., the only engine that is usable in Particles attributes) is
+MultiPatch.
+
+<p>
+MultiPatch can use several different types of layouts and single-patch
+engines, and all MultiPatch engines use a shared layout. However, only
+the MultiPatch<DynamicTag,*> specializations of MultiPatch engines are
+useful for Particles attributes, since only that engine type can have
+patches of dynamically varying size.
+
+<li>Can the engine change size dynamically?
+
+<p>
+The engine type used for particle attributes must have dynamic
+capabilities. Thus, we should use dynamic single-patch engines
+inside of MultiPatch. The only engines available in this release
+of POOMA that meet this requirement are Dynamic and Remote<Dynamic>.
+Both of these are inherently one-dimensional and support operations
+such as create(), destroy() and copy(). Remote<Dynamic> is similar
+to Dynamic, but it is context-aware and useful for multi-context codes.
+
+<p>
+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:
+
+<ol>
+<li> the layout for the particle attributes;
+
+<li> 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
+
+<li> the actual SpatialLayout that connects the info about the Field
+layout to the Particles attribute layout.
+</ol>
+
+<p>
+The only thing that needs to match between the attribute and Field
+layouts is the number of patches, which must be exactly the same.
+The engine type (and thus the layout type) of the attributes
+and of the Field do not have to match. Typically, both the attributes
+and the Field will have a MultiPatch engine with the same number of
+patches, but these engines will have different single-patch engine
+types (Dynamic vs. Brick) and use different types of layouts (Dynamic
+vs. Grid).
+
+<p>
+Note once again that in the simple case of a UniformLayout particle
+layout, applications do not need to worry about any Field layout type,
+only the particle attribute layout (which still needs to be shared)
+and the particle layout. This commonly arises during the early
+prototyping (i.e., pre-parallel) stages of application development,
+when you might limit an application to a single patch for simplicity.
+
+</ol>
+
+
+<h3>Example: Particle-in-Cell Simulation</h3>
+
+<p>
+Our third and final example of the Particles 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 (discussed elsewhere)
+with the Particles capabilities we have been describing here.
+
+<p>
+Because this example is longer than the others in this document,
+it will be described in sections. For a unified listing of the source
+code, please see the file examples/Particles/PIC2d/PIC2d.cpp.
+
+The first step is to include all of the usual header files:
+
+<pre>
+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>
+</pre>
+
+<p>
+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 typedef's 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.
+
+<pre>
+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 };
+</pre>
+
+<p>
+The interpolator tag type is included in the traits class because an
+application might want the Particles subclass 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 subclass. 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;
+
+<p>
+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.
+
+<p>
+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
+particle positions R, their velocities V, and their charge/mass
+ratios qm. The class also has a fourth attribute called E, which is
+used to record the electric field at each particle's position in order
+to calculate forces. This calculation will be discussed in greater
+detail below.
+
+<pre>
+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 {
+043 addAttribute(R);
+044 addAttribute(V);
+045 addAttribute(E);
+046 addAttribute(qm);
+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 };
+</pre>
+
+<p>
+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 changed), the type of mesh on which the
+calculations are done, the mesh's size, and so on:
+
+<pre>
+058 // Dimensionality of this problem
+059 static const int PDim = 2;
+060
+061 // Engine tag type for attributes
+062 typedef MultiPatch<DynamicTag,Dynamic> 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;
+</pre>
+
+<p>
+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:
+
+<pre>
+115 int main(int argc, char *argv[])
+116 {
+117 // Initialize POOMA and output stream.
+118 Pooma::initialize(argc,argv);
+119 Inform out(argv[0]);
+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);
+</pre>
+
+<p>
+The program then creates a pair of Field objects. The first, phi,
+is a field of doubles and records the electrostatic potential at
+points in the mesh. The second, EFD, is a field of two-dimensional
+Vectors and stores 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.
+
+<pre>
+133 // Create field layout objects for our electrostatic potential
+134 // and our electric field. Decomposition is 4 x 4 and replicated.
+135 Loc<PDim> blocks(4,4);
+136 FLayout_t flayout(geometry.physicalDomain(),blocks,ReplicatedTag()),
+137 flayoutGL(geometryGL.physicalDomain(),blocks,gl,ReplicatedTag());
+138
+139 // Create and initialize electrostatic potential and electric field.
+140 DField_t phi(geometryGL,flayoutGL);
+141 VecField_t EFD(geometry,flayout);
+</pre>
+
+<p>
+The application now adds periodic boundary conditions to the electrostatic
+field phi, so that particles will not see sharp changes in the 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.
+
+<pre>
+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)
+147 phi.addBoundaryConditions(AllPeriodicFaceBC());
+148 double phi0 = 0.01 * nx;
+149 phi = phi0 * sin(2.0*pi*phi.x().comp(0)/nx)
+150 * cos(4.0*pi*phi.x().comp(1)/ny);
+151 EFD = -grad<Centering_t>(phi);
+</pre>
+
+<p>
+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.
+
+<pre>
+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);
+160 P.addBoundaryCondition(P.R,bc);
+161
+162 // Create an equal number of particles on each processor
+163 // and recompute global domain.
+164 P.globalCreate(NumPart);
+</pre>
+
+<p>
+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 would be:
+
+<pre>
+Particles_t::PointType_t lower, upper;
+for (int d=0; d<PDim; ++d)
+{
+ lower(d) = 0;
+ upper(d) = ng[d];
+}
+</pre>
+
+<p>
+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.
+
+<pre>
+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;
+</pre>
+
+<p>
+The application is now able to enter its main timestep loop.
+In each time step, the particle 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
+particle velocities can be updated:
+
+<pre>
+195 // Begin main timestep loop
+196 for (int it=1; it <= NumSteps; ++it)
+197 {
+198 // Advance particle positions
+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 }
+</pre>
+
+<p>
+Finally, the state of the particles at the end of the simulation is
+printed out, and the simulation is closed down:
+
+<pre>
+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;
+</pre>
+
+
+<h2>Summary</h2>
+
+<p>
+This document 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.
+
+</body>
+</html>
Index: Layout.html
===================================================================
RCS file: /home/pooma/Repository/r2/docs/Layout.html,v
retrieving revision 1.3
diff -u -u -r1.3 Layout.html
--- Layout.html 20 Aug 2004 20:14:18 -0000 1.3
+++ Layout.html 23 Aug 2004 11:16:32 -0000
@@ -5,8 +5,11 @@
<meta name="GENERATOR" content="Mozilla/4.72 [en] (X11; U; Linux 2.2.14 i686) [Netscape]">
<title>Layout and related classes</title>
</head>
-<body text="#000000" bgcolor="#FFFFFF" link="#FF0000" vlink="#800080">
-
+<body background="back.gif" LINK="#505062" ALINK="#505062" VLINK="#be7c18">
+
+<CENTER><IMG SRC="banner.gif" ALT="POOMA banner" WIDTH=550 HEIGHT=100
+ALIGN=bottom></CENTER>
+
<h1>
Layouts and related classes:</h1>
More information about the pooma-dev
mailing list