ippl_particles.tex 41.4 KB
Newer Older
gsell's avatar
gsell committed
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20
\chapter{Particles}
\label{sec:particles}
This section describes the \ippl framework classes which provide the capability to performing particle-based simulations. We first describe how to design and instantiate \texttt{Particle} classes customized to the needs of a specific application, and then discuss the possible operations and expressions in which a particle object may be employed and end with an ready to use example.

\section{Basic \texttt{Particle} Object Characteristics}

The \ippl framework treats \texttt{Particle} classes as containers which store the characteristic data for \texttt{N} individual particles. Each particle has several attributes, such as position, mass, velocity, etc. Looked at in another way, \texttt{Particle} classes store several attribute containers, where each attribute container holds the value of that attribute for all \texttt{N} particles. \texttt{Particle} objects in \ippl may be thought of as shown in the following diagram: ...

There are two particle attributes predefined, namely R (position) and ID a global unique identifier.

The data type of each attribute, the number of attributes, and the names for these attributes are completely customizable, and are specified in the manner described in the following sections. Any number of different \texttt{Particle} classes may be defined and used within a given simulation. Also, the \texttt{Particle} objects may interact with \ippl \texttt{Field} objects or may be used independently. In addition to the attributes, each \texttt{Particle} object uses a specific layout mechanism, which describes the data of the individual particles is spread across the processors in a parallel environment. The \ippl framework provides several different \texttt{Particle} layout classes, any of which may be selected to partition the particle data among processors. The choice of layout depends on the intended use of the \texttt{Particle} object, as discussed later. Once defined and instantiated, \texttt{Particle} objects in the \ippl framework may be used in many ways, including:
\begin{itemize}
    \item Operations involving all the particles within a \texttt{Particle} object may be specified using simple expressions, in a manner very similar to that used for \texttt{Field} objects. These expressions may involve any of the attributes of the particles as well as other scalar data, and they may use not only the standard mathematical operators +, -, *, /, etc., but also standard mathematical functions such a s cos ( ) , exp ( ) , mod ( ) , etc.
    \item Alternatively, you may set up explicit loops that perform operations involving the attributes of a single particle or a subset of all the particles.
    \item \texttt{Particle}s may be created or destroyed during a simulation.
    \item \texttt{Particle}-to-\texttt{Field} and \texttt{Field}-to-\texttt{Particle} operations may be performed (e.g., a particular \texttt{Particle} attribute may be deposited onto a specified \texttt{Field} using a chosen interpolation method). 
\end{itemize}

\section{Defining a User-Specified \texttt{Particle} Class}

21
There is no specific class within the \ippl framework called \texttt{Particle}. Rather, the first step in deploying particles within a \ippl application is to define a user-specified \texttt{Particle} class, which contains the attributes required for each particle, as well as any, specific methods or data the user may need. To do this, the \texttt{IpplParticleBase} and \texttt{ParticleAttrib} classes are used, along with a selected subclass of the \texttt{ParticleLayout} class. The steps to follow in creating a new \texttt{Particle} class are:
gsell's avatar
gsell committed
22 23 24
\begin{itemize}
    \item Based on the type of interactions which the particles have with each other and with external objects such as a \texttt{Field}, select a method of distributing the particles among the nodes in a parallel machine.
    \item Next, decide what attributes each particle should possess.
25 26
    \item Third, create a subclass of \texttt{IpplParticleBase} which includes these attributes (specified as instances of the \texttt{ParticleAttrib} class template).
    \item Finally, instantiate this user-defined subclass of \texttt{IpplParticleBase} and create and initialize storage for the particles which are to be maintained by this object. 
gsell's avatar
gsell committed
27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83
\end{itemize}

The following sections describe in more detail how to accomplish these steps.

\subsection{Selecting a Layout: \texttt{ParticleLayout} and Derived Classes}

When used in a parallel environment, the \ippl framework partitions the particles in a \texttt{Particle} container among the separate processors and includes tools to spread the work of computing and the results of expressions involving \texttt{Particle} attributes among the processing nodes. There are, however, different ways in which particles may be distributed among the processors, and the method which should be used depends upon how the particles in a \texttt{Particle}
object will interact with each other and with \texttt{Field} objects (if at all). The \ippl framework includes different \texttt{Particle} layout mechanisms, which are all derived from the \texttt{ParticleLayout} class. Each \texttt{Particle} object needs its own layout object; that is, you cannot create a layout object and give it to more than one \texttt{Particle} object. The methods typically used to determine how to assign particles to particular nodes are based on analysis of the
position (\texttt{R} attribute) of each particle. Thus, \texttt{ParticleLayout} and its derived classes have two template parameters: the type and the dimensionality of the particle position attribute (this particle position attribute is discussed in more detail later). The following sections describe the particle layout mechanisms currently available in the \ippl framework.

\subsection{The \texttt{ParticleUniformLayout} Class}

The \texttt{ParticleUniformLayout} class maintains an equal number of particles on each node, with no consideration of particle position. \texttt{ParticleUniformLayout} is useful in those cases where particles do not interact with each other but perhaps with some other external agent, so that no consideration need be made about which particles are located near to others. In that case, with this layout maintains an equal balance of memory usage among the processors and requires relatively small amounts of interprocessor communication. If you require the ability to compute an interaction between a particle and its nearest neighbors, this is not the proper layout to use, in that case, the \texttt{Particle}SpatialLayout class (discussed next) is a better choice. The constructor for \texttt{ParticleUniformLayout} takes no arguments, but does require the two template parameters that specify the type and dimensionality of the particle position attribute. An example of creating a new \texttt{ParticleUniformLayout} instance for a 3D particle simulation is:
\begin{smallcode}
ParticleUniformLayout<double,3> uniformlayout();
\end{smallcode}


\subsection{The \texttt{ParticleSpatialLayout} Class}

\texttt{ParticleSpatialLayout}, in contrast to \texttt{ParticleUniformLayout}, assigns particles to nodes based upon their spatial location relative to a \texttt{FieldLayout}. It is useful when the particles will be interacting with other particles in their neighborhood or with a \texttt{Field} object. \texttt{ParticleSpatialLayout} will keep a particle on the same node as that which contains the section of the \texttt{Field} in which the particle is located. If the particle moves to a new position, this layout will reassign it to a new node when necessary. This will maintain locality between the particles and any \texttt{Field} distributed using this \texttt{FieldLayout}. Further more it will help keep particles which are spatially close to each other local to the same processor as well. As with all the layout classes, \texttt{ParticleSpatialLayout} requires the type and dimensionality of the particle position attribute as template parameters. The constructor for \texttt{ParticleSpatialLayout} takes one argument: a pointer to a \texttt{FieldLayout} object that tells the \texttt{ParticleSpatialLayout} how the \texttt{Field} is allocated among the parallel processors, so that the particles may be maintained local to this \texttt{Field}. Note that you do not, need to create a \texttt{Field} instance itself, you only need to give \texttt{ParticleSpatialLayout} a \texttt{FieldLayout} object. An example of creating an instance of this class is as follows:
\begin{smallcode}
FieldLayout<3> myfieldlayout(Index(l6), Index(16), Index(32)); 
ParticleSpatialLayout<double,3> myparticlelayout(&myfieldlayout);
\end{smallcode}


Note that the dimensionality of the \texttt{FieldLayout} and the \texttt{ParticleSpatialLayout} (in this example, 3) must be the same. You may also create a \texttt{ParticleSpatialLayout} instance without providing a \texttt{FieldLayout}. In this case, particles will remain on the node on which they were created. If at some future time you wish to provide a \texttt{FieldLayout} object to tell the \texttt{ParticleSpatialLayout} where to place the particles, you may do so using the \texttt{setFieldLayout} (\texttt{FieldLayout<Dim>*}) method of \texttt{ParticleSpatialLayout}. This is useful when reading particles in from an external source and the size of the spatial domain containing the particles is not known until all the particles have been read. The following example demonstrates the use of the capability:
\begin{smallcode}
ParticleSpatiaILayout<double,3> myparticleLayout; 
// calculate the size of the domain required to contain all the particles 
// create a new FieldLayout object based on these calculations 
FieldLayout<3> myfieldlayout(Index(minx, maxx), Index (miny, maxy), 
                             Index(minz,maxz);   
myparticlelayout.setFieldLayout(&myfieldlayout);
\end{smallcode}

\texttt{ParticleSpatialLayout} also provides functionality to maintain cached ghost particles from neighboring nodes which might be required for particle - particle interaction. A caching policy can be defined using the fourth template parameter of \texttt{ParticleSpatialLayout}:
\begin{smallcode}
typedef UniformCartesian<Dim, double> Mesh_t;
typedef ParticleSpatialLayout<double,Dim,Mesh_t,
        BoxParticleCachingPolicy<double, Dim, Mesh_t> > playout_t;
\end{smallcode}
The available chaching policies are: \texttt{NoParticleCachingPolicy},\texttt{BoxParticleCachingPolicy} and \texttt{CellParticleCachingPolicy}. With \texttt{NoParticleCachingPolicy} there is no caching whatsoever. \texttt{BoxParticleCachingPolicy} extends the interface of \texttt{ParticleSpatialLayout} by two functions \texttt{void setCacheDimension(int d, T length)} and \texttt{void setAllCacheDimensions(T length)} which are used to set the size of the cached region around the local domain in units of space. \texttt{CellParticleCachingPolicy} extends the interface of \texttt{ParticleSpatialLayout} by two functions \texttt{void setCacheCellRange(int d, int length)} and \texttt{void setCacheCellRanges(int d, int length)} which are used to set the size of the cached region around the local domain in units of grid cells of the mesh. \texttt{BoxParticleCachingPolicy} is the default policy.

The caching can be enabled or disabled by calling the \texttt{enableCaching()} or \texttt{disableCaching()} member functions of \texttt{ParticleSpatialLayout}. Caching is disabled by default.


\subsection{Selecting Particle Attributes: The \texttt{ParticleAttrib} \\Class}

\texttt{ParticleAttrib} is a class template that represents a single attribute of the particles in a \texttt{Particle} object. Each \texttt{ParticleAttrib} contains the data for that attribute for all the particles. Within a user-defined \texttt{Particle} class, you declare an instance of \texttt{ParticleAttrib} for each attribute the particles will possess and assigns to it an arbitrary name. \texttt{ParticleAttrib} requires one template parameter, the type of the data for the attribute. As an example, the statement:
\begin{smallcode}
ParticleAttrib<double> density; 
\end{smallcode}

declares an instance of \texttt{ParticleAttrib} named 'density', which will store a quantity of type double for all the particles of the \texttt{Particle} class that contains this data member.

84
\subsection{Specifying a User-Defined \texttt{Particle} Class: \\The \texttt{IpplParticleBase} Class}
gsell's avatar
gsell committed
85

86
\texttt{IpplParticleBase} is the class that all user-defined \texttt{Particle} classes must specify as their base class. It stores the list of attributes for the particles (which are maintained as instances of \texttt{ParticleAttrib}) and a selected parallel layout mechanism. In addition to providing all the capabilities for performing operations on the particles and their attributes, \texttt{IpplParticleBase} also defines two specific attributes which all user-defined \texttt{Particle} classes inherit:
gsell's avatar
gsell committed
87 88 89 90 91 92
\begin{smallcode}
ParticleAttrib<Vektor<T,Dim>>  R; 
ParticleAttrib<unsigned>      ID;
\end{smallcode}


93 94 95
The first attribute, \texttt{R}, represents the position of each particle. Each position is stored as a \texttt{Vektor<T, Dim>}, which is a \ippl data type representing a dim-dimensional vector with elements of type \texttt{T}. The second attribute, \texttt{ID}, stores a unique unsigned integer value for each particle. The values are not guaranteed to be in any particular order, but they are guaranteed to be unique for each particle. \texttt{IpplParticleBase} has one template parameter, the layout class to be used to assign particles
to processors (e.g., \texttt{ParticleSpatialLayout}). The data type and dimensionality of the particle position attribute (\texttt{R}) will be the same as those used to create the specific \texttt{ParticleLayout} derived class. Each \texttt{IpplParticleBase} contains one instance of the chosen layout class. There are two constructors for \texttt{IpplParticleBase}: a default constructor that creates a new instance of the layout class using the layout's default constructor, and a constructor which takes a pointer to an instance of the
layout class. The second version of the \texttt{IpplParticleBase} constructor is useful when the desired layout class requires arguments to its constructor (e.g., \texttt{ParticleSpatialLayout}, which may be give in a \texttt{FieldLayout} pointer).
gsell's avatar
gsell committed
96

97
Using \texttt{IpplParticleBase}, \texttt{ParticleAttrib}, and a selected class derived from \texttt{ParticleLayout,} you can create a user-defined \texttt{Particle} class using the following code template: \\
gsell's avatar
gsell committed
98 99
\clearpage
\begin{codeln}
100
class Bunch : public IpplParticleBase< ParticleSpatiaILayout<double,3> > 
gsell's avatar
gsell committed
101 102 103 104 105 106 107
{ 
public: 
    // Attributes for this particle class (besides position and ID). 
    ParticleAttrib<double>             qm;      // q/m ratio 
    ParticleAttribs Vektor<double,2> > vel; 	// velocity 

    // constructor 
108
    Bunch(Layout\_t *L) : IpplParticleBase<Layout\_t>(L) { 
gsell's avatar
gsell committed
109 110 111 112 113 114 115 116
        addAttribute(qm); 
        addAttribute(vel); 
    } 
}; 
\end{codeln}

Let us describe this example in detail by discussing the important lines in the order of use.

117
Line 1: You may select whatever name is appropriate for the specialized \texttt{Particle} class, but it must be derived from \texttt{IpplParticleBase}. 
gsell's avatar
gsell committed
118

119
In this case, we explicitly specify the type of layout to use (\texttt{ParticleSpatialLayout}), with particle position attribute type and dimensionality template parameters of double and 3, respectively. Alternatively, \texttt{Bunch} may have been declared as a class template itself and may have passed on the layout template parameters to \texttt{IpplParticleBase}. In that case, the first line would instead look like
gsell's avatar
gsell committed
120 121
\begin{smallcode}
template <class PLayout> 
122
class Bunch : public IpplParticleBase<PLayout> 
gsell's avatar
gsell committed
123 124
\end{smallcode}

125
Lines 5-6: Here is where the attributes for the particles in the \texttt{Particle} object are declared. They may be given any name other than \texttt{R} or \texttt{ID}. Instead of stating the type and dimensionality of this attribute specifically, you may also use one of the following typedefs and constants defined in \texttt{IpplParticleBase}:
gsell's avatar
gsell committed
126 127 128 129 130 131 132 133 134 135
\begin{itemize}
          \item \texttt{Dim} - the dimensionality of the particle position attribute (in this example, 3)
          \item \texttt{Position\_t} - the type of data used to store the position attribute components (here, this type is double)
          \item \texttt{Layout\_t} - a synonym for the specified layout class
          \item \texttt{ParticlePos\_t} - a typedef for the particle position attribute; it is shorthand for \texttt{ParticleAttrib< Vektor<Position\_t,Dim> >}
\end{itemize}
and could have been used to specify the attribute \texttt{vel} in the above example as \texttt{ParticlePos\_t vel};
\begin{itemize}
\item \texttt{ParticleIndex\_t} - a typedef for the particle global \texttt{ID} attribute; it is short for \texttt{ParticleAttrib<unsigned>} 
\end{itemize}
136
The constructor for this user-defined class must initialize \texttt{IpplParticleBase} with a pointer to an instance of the selected layout class. 
gsell's avatar
gsell committed
137 138 139 140 141 142

In this example, the layout class is \texttt{ParticleSpatialLayout}, but using one of the typedefs listed above, we can abbreviate this as \texttt{Layout\_t}. Note that we only define one constructor here, omitting the default constructor. This is done because \texttt{ParticleSpatialLayout} (which we have hard-coded as the layout for this user-defined \texttt{Particle} class) requires an argument to its constructor, and this can only be provided if we use a constructor for our \texttt{Particle} class as shown here. A new instance of this class would be declared in an application as follows:
\begin{smallcode}
Bunch myBunch (new ParticleSpatialLayout<double,3>(myFieldLayout)); 
\end{smallcode}
where my\texttt{FieldLayout} was a \texttt{FieldLayout} object created previously. The only action that is required in the constructor for the derived class is to inform the base class of the declared attributes, 
143
using the \texttt{addAttribute(.)} method of \texttt{IpplParticleBase}, which registers the specified \texttt{ParticleAttrib} instance with the parent class \texttt{IpplParticleBase}. The order in which attributes are registered is not important.
gsell's avatar
gsell committed
144 145

\subsection{Example \texttt{Particle} Classes: The \texttt{Genparticle} and \texttt{GenArrayParticle} Classes}
146 147
The \ippl framework provides two classes which are examples of \texttt{Particle} classes derived from \texttt{IpplParticleBase}: \texttt{Genparticle} and \texttt{GenArrayParticle}. They may be used as samples from which to build new classes, or they may be used to quickly include particle capabilities in an application. \texttt{Genparticle} is a \texttt{Particle} class with three attributes: \texttt{R} and \texttt{ID} inherited from \texttt{IpplParticleBase}, and an attribute named data with elements
of an arbitrary type. \texttt{Genparticle} has two template parameters: the type of particle layout to use and the type \texttt{T} of attribute data. It has two constructors just as \texttt{IpplParticleBase} does: the default and one taking a layout pointer. An example of instantiating a \texttt{GenParticle}object is shown below.
gsell's avatar
gsell committed
148 149 150 151 152 153 154 155 156 157 158 159 160 161 162 163 164 165 166
\begin{smallcode}
GenParticle<ParticleUniforrnLayout<float,3>,UserDefinedType> GP(); 
\end{smallcode}


\texttt{GenArrayParticle} is almost identical to \texttt{GenParticle}, the difference being that \texttt{GenArrayparticle} contains not just one but an array of attributes data \texttt{[0 ... N-l]} of a specified type. The number of elements in the attribute array, \texttt{N}, is given as a third template parameter. The following example shows a \texttt{GenArrayParticle} being created with 5 floats stored in the data array for each particle:
\begin{smallcode}
GenArrayParticle<ParticleSpatialLayout<doulble,3>,float,5> GAP(
        new ParticleSpatialLayout<double,3>(myFieldLayout)); 
\end{smallcode}

It is important to note that the array data in \texttt{GenArrayParticle} contains a set of particle attributes of the same type. In situations where it is necessary to have a variety of particle attribute types, you may use the \texttt{Genparticle} class with the type of data being specified as a user-defined struct containing the various attributes needed.

\subsection{Using \texttt{Particle} Classes in an Application}

After a specific \texttt{Particle} class has been defined and created in a \ippl application, you may create and initialize new particles, delete unwanted particles, and perform computations involving these particles. This section describes how to accomplish these tasks.

\subsection{Creating New \texttt{Particle}s}

167
When a \texttt{Particle} object is created, it is initially empty. Storage for new particles is allocated using the create (unsigned) method of \texttt{IpplParticleBase}. For example, if a \texttt{Particle} object bunch has been created already, the statement
gsell's avatar
gsell committed
168 169 170 171 172 173 174
\begin{smallcode}
bunch.create(100); 
\end{smallcode}


will allocate storage for 100 new particles. All the attributes for the particles in the \texttt{Particle} object will have this new storage allocated. The data is uninitialized, except for the global \texttt{ID} attribute; you must assign the proper values to the position and any other attributes that have been defined. The new storage is appended to the end of any existing storage.

175 176
\texttt{IpplParticleBase} includes two methods that allow you to query how many particles exist. The function \texttt{getTotalNum()} will return the total number of particles being stored among all the processors; the function \texttt{getLocalNum()} will return the number of particles just on the local node. Although the new storage space is allocated on the local processor on which the call to create was executed, the \texttt{Particle} class will not officially add the particles to its
local count (and will not tell any other processors it has created these new particles) until you call the \texttt{update()} method of \texttt{IpplParticleBase}. Thus, a call to \texttt{getLocalNum()} will report the same number just before and just after the call to create. The storage does exist after create is called, but only after the \texttt{update} method (which is discussed in more detail in a later section) has been called will all the processors have correct information on their local and total particle counts.
gsell's avatar
gsell committed
177 178 179 180 181 182 183 184 185 186 187 188 189 190 191 192 193 194 195 196 197 198 199 200

\subsection{Initializing Attribute Data}

After calling create to allocate new storage, you must initialize the data. This should be done after calling create and before calling \texttt{update} for the \texttt{Particle} object. After the data is initialized, the \texttt{update} routine will properly distribute the particles to their correct node based on the layout mechanism chosen for that \texttt{Particle} object and possibly the positions of the particles as set during their initialization. The following example shows one way to initialize the data for newly
created particles when running on a single-processor machine. (This example will be modified in the following section for the case of running in parallel.) \\
\begin{code}
// create and'initialize data for an instance of Bunch 
Bunch myBunch(new Bunch::Layout\_t(myFieldLayout)); 
int currLocalPtcls = myPtcls.getLocalNum(); 
myBunch.create(100); 
for (int i = 0; i < 100; i++) { 
    myBunch.R[currLocalPtcls + i] = Vektor<double,3>(0.0, 1.0, 0.0); 
    myBunch.vel[currLocalPtcls + i] Vektor<double,3>(1.0, 1.0, 1.0); 
} 
myBunch.update(); 
\end{code}


In this example, 100 new particles are created, and the \texttt{R} and \texttt{vel} attributes are initialized to \texttt{Vektor} quantities. Notice that each attribute is accessed simply by specifying it as a data member of the \texttt{myBunch} object. After create was called, even though the 100 particles were not added to the \texttt{Particle} object's count of local particles, the storage was allocated and it was possible to assign values to the new elements in the attribute storage
(accessed simply using the \texttt{[]} indexing operator). Finally, calling \texttt{update} added the new storage to the count, of particles stored in \texttt{myBunch}. Further calls to getLocalNum and getTotalNum would report the proper values.

\subsection{Initializing Attribute Data on Parallel Architectures}

The code shown in the previous example has one problem when used on parallel architectures: the call to create is performed on each processor, so if there were P processors a total of 100*P particles would be created. This may be the desired behavior, if so, the previous example is sufficient. However, if you are reading data on particle positions and other attributes from a file or some other source, you may wish to create particles on a single processor, and then distribute the
201
data to the proper nodes. To do this, you need to call create and assign initial data on only one node but call update on all the processors. The \texttt{singlelnitNode()} method of \texttt{IpplParticleBase} will return a boolean value indicating whether the local processor should be used to create and initialize new particles in this way. The following example demonstrates how to use this method for initializing particles: \\
gsell's avatar
gsell committed
202 203 204 205 206 207 208 209 210 211 212 213 214 215 216 217 218
\begin{code}
// create and'initialize data for an instance of Bunch 
Bunch myBunch(new Bunch::Layout\_t(myFieldLayout)); 
int currLocalPtcls = myPtcls.getLocalNum();
if (myBunch.singleInitNode()) { 
    myBunch.create(100);  
    for (int i = 0; i < 100; i++) { 
        myBunch.R[currLocalPtcls + i] = Vektor<double,3>(0.0, 1.0, 0.0); 
        myBunch.vel[currLocalPtcls + i] Vektor<double,3>(1.0, 1.0, 1.0); 
    }
} 
myBunch.update(); 
\end{code}


\subsection{Deleting \texttt{Particle}s}

219
\texttt{Particle}s may also be deleted during a simulation. The method \texttt{destroy (unsigned M, unsigned I)} of \texttt{IpplParticleBase} will delete \texttt{M} particles, starting with the \texttt{I}th particle. The index \texttt{I} here refers to the local particle index, not the global \texttt{ID} value. Thus \texttt{I = 0} means delete particles starting with the first one on the local processor.
gsell's avatar
gsell committed
220 221 222 223 224 225

Unlike the situation when creating new particles, the storage locations for the deleted particles will not be removed from attribute data storage until \texttt{update} is called. Instead, the requests to delete particles are cached until the update phase, at which time all the deletions are performed. You are allowed to issue multiple delete requests between \texttt{update}s. For example, if there are 100 particles on a local node, and you request to delete particles 0 to 10
and then request to delete particles 60 to 70, nothing will change in the attribute storage until you call \texttt{update}, and no change will occur to the local and total particle counts until \texttt{update()} is complete.

\subsection{Updating \texttt{Particle}s: The \texttt{update()} Method}

226
The \texttt{update()} method of \texttt{IpplParticleBase} is responsible for making sure that all processors have correct information about how many particles exist and where they are located in a parallel machine. As mentioned previously, \texttt{update} must be called by all processors after a sequence of particle creation or deletion operations. The \texttt{update} method is also responsible for maintaining a proper assignment of particles to processors, based on the particular
gsell's avatar
gsell committed
227
\texttt{ParticleLayout} class used to create the
228
\texttt{IpplParticleBase} object. Typically, this layout mechanism depends on the position of particles, so when particles change their position, they may need to be reassigned to a new processor to maintain the proper layout. In this case, the \texttt{update} method should be called whenever a computation is complete which alters the attributes (e.g. position) that a layout depends upon. The following short example demonstrates using \texttt{update} in conjunction with some operation that alters the x-coordinate of a
gsell's avatar
gsell committed
229 230 231 232 233 234 235 236 237 238 239 240 241 242 243 244 245 246 247 248 249 250 251 252 253 254 255 256 257 258 259 260 261 262 263 264 265 266 267 268 269 270 271 272 273 274 275 276 277 278 279 280 281 282 283 284 285 286 287 288 289 290 291 292 293 294 295 296 297 298
set of particles. \\
\begin{code}
// do some computation involving myBunch for several time steps 
while (computation_done == false) { 
    // for each particle, add some constant to the x coordinate 
    myBunch.R(0) += 0.li 
    // update the Particle object; this may move particles between nodes 
    myBunch.update(); 
    // determine if the computation is done, etc.
}
\end{code}

\subsection{Using Particle Attributes in Expressions}

Computations involving particle attributes may be performed in many ways. Data-parallel expressions that involve all particles of a given \texttt{Particle} object may be used, or specific loops may be written that employ attribute iterators or nearest-neighbor pairlist iterators.

\subsubsection{Attribute Expressions}

Just as with the \texttt{Field} class, you may perform data-parallel operations on particle attributes using a simple expression syntax, which the \ippl framework will translate into efficient inlined code. These operations will be performed for every particle. The expressions may include any of the attributes in a \texttt{Particle} object as well as scalar values, may use mathematical operators such as +, -, *, / etc., and may call standard mathematical functions such as \texttt{cos(
)}, \texttt{exp( )}, \texttt{mod( )} , etc. for an attribute value of each particle. Some examples are shown below.
\begin{smallcode}
double dt = 2.0; 
myBunch.R += myBunch.vel* dt; 
myBunch. vel = 1. ¡ - log (1. ¡ + myBunch. R * myBunch. R) ; 
myBunch.update(); 
\end{smallcode}


Attribute expressions will perform their operations on all the particles in the \texttt{Particle} object, including any new particles allocated via a call to create, even before update has been called. This fact is useful when initializing the attributes for newly created particles (e.g., to set the init value for some scalar quantity to zero). Generally, however, unless you are performing an initialization of new particles, you should avoid using particle expressions of this type after calls to create or destroy and before a call to update.

Some attributes, such as \texttt{Vektor}s or \texttt{Tenzor}s, have multiple components, and you may wish to involve only the \texttt{N}th component of the attribute in an expression. To do so, use the \texttt{()} operator to select the \texttt{N}th component of that attribute. For instance, using \texttt{myBunch} from the previous example, you can change just the x-coordinate of the particle position attribute \texttt{R} as follows: 
\begin{smallcode}
myBunch.R(0) = myBunch.R(l) - cos(myBunch.R(2));
\end{smallcode}


For 2D or 3D quantities, use two or three indices. For example, if rho is a 3x3 \texttt{Tenzor} attribute of myBunch, you can do the following:
\begin{smallcode}
myBunch.rho(0,0) = -(myBunch.rho(0,l) + myBunch.rho(0,2)); 
\end{smallcode}

Attribute expressions may also use the where operator in much the same way as for \texttt{Field} expressions. The first argument to where is some expression that results in a \texttt{boolean} value for each particle. The second and third arguments are expressions that will be evaluated for a particle if the first argument is \texttt{true} or \texttt{false}, respectively, for that particle. For example,
\begin{smallcode}
myBunch.vel = where(myBunch.R(0) > 0.0, -2.0 * myBunch.vel, myBunch.vel) 
\end{smallcode}

changes the value of the \texttt{vel} attribute in \texttt{myBunch} when the x-coordinate of the particle position is positive.

\subsection{\texttt{Particle} Iterator Loop}

You also have the capability of performing operations on specific particles using iterators or standard indexing operations. The \texttt{ParticleAttrib} containers in a \texttt{Particle} class may be used just as regular \texttt{STL} containers. The \texttt{begin()} and \texttt{end()} methods of the \texttt{ParticleAttrib} class will return an iterator pointing to the first element and just past the last element, respectively, of the attribute. These iterators may be used in an explicit loop just as if they were pointers into the attribute array.
\begin{smallcode}
ParticleAttrib<unsigned>::iterator idptr, idend = myBunch.ID.end(); 
for (idptr = myBunch.ID.begin(); idptr != idend; ++idptr) 
    cout << "Particle ID value: " << *idptr << endl; 
\end{smallcode}


Iterators are available for all \texttt{ParticleAttribs}. As an alternative, you may simply use the \texttt{[]} operator to access the attribute data of the \texttt{N}th particle on a node, treating \texttt{ParticleAttrib} as a regular array of data.
\begin{smallcode}
int nptcls = myBunch.getLocalNum(); 
for(int i=0; i < nptcls ++i) { 
    cout << "Particle ID value: " << myPtcls.ID[i] << endl; 
}
\end{smallcode}


\section{Nearest-Neighbor Interactions (Jakob)}

%A standard need for many types of particle simulation is the execution of some calculation involving an interaction between one particle and all other particles within a given interaction region or radius. The interaction may involve particles that are on other processors, so it is important to efficiently communicate the necessary data between the nodes and limit the amount of work used in determining the nearest-neighbor interaction pairlists. For this type of computation,
299
%\texttt{IpplParticleBase} provides a specialized iterator, called \texttt{pair\_iterator}, which allows you to iterate over all neighbors of a given particle. The specific layout mechanism used by a \texttt{Particle} object takes care of calculating which particles lie within a selected interaction radius of each other, taking into account nearby particles which may reside on other processors in a parallel machine. The first step in calculating nearest-neighbor interactions is to specify what constitutes the
gsell's avatar
gsell committed
300
%neighborhood. Currently, only a radially symmetric interaction neighborhood is supported in \ippl. The interaction radius may be the same for each particle, or it may be a scalar attribute of the \texttt{Particle} class (of the same type as that used for the particle coordinates, e.g., \texttt{double} or \texttt{float}). Before any nearest-neighbor calculation is performed, you must specify the interaction radius of the particles. This is done using the
301
%\texttt{setInteractionRadius()} method of \texttt{IpplParticleBase}, giving as an argument either a single scalar value that will be used for all the particles, or a \texttt{ParticleAttrib<T>} quantity. In the latter case, each particle can have a different interaction radius. \texttt{IpplParticleBase} will remember this data and use it in the calculation of interaction lists. You may call these methods more than once; the most recent value given for the interaction radius is used to determine the
gsell's avatar
gsell committed
302 303 304 305 306 307 308
%nearest-neighbor pairlists. After changing the interaction radius, you should as always perform an \texttt{update()} to make sure all nodes have correct information on how to find neighboring particles. You may, for example, update the \texttt{ParticleAttrib} used to store the individual particle interaction radii as the last step of some calculation, and then call \texttt{update()}. You can retrieve the interaction radius for any local particle by calling the method
%\texttt{getlnteractionRadius(unsigned)} of the \texttt{Particle} object containing that particle.

%The second step in this process is to retrieve a list of nearest-neighborparticles for a given local particle, and then to perform the desired calculation using those nearby particles. This is done using the 
%\begin{smallcode}
%getpairlist(unsigned, pair_iterator&, pair_iterator&)
%\end{smallcode}
309 310
%method of \texttt{IpplParticleBase}. The first argument is the local index of the central particle under consideration. The second and third arguments are, set equal to begin and end
%\texttt{pair\_iterator} objects for the specified particle. The type \texttt{pair\_iterator} is defined in \texttt{IpplParticleBase}. This set of \texttt{pair\_iterator} objects is then used to retrieve the local indices of all particles within the interaction region of the central particle.
gsell's avatar
gsell committed
311 312 313 314 315 316 317 318 319 320 321 322 323 324 325 326 327

%Here is an example of how to use this mechanism: \\
%\clearpage
%\begin{codeln}
%myBunch.setlnteractionRadius(4.0); //determines interaction region 
%Bunch::pair_iterator neighbor, neighbor_end;
%for(int i=0; i < myBunch.getLocalNum(); ++i;) { 
%    myBunch.getPairlist(i, neighbor, neighbor_end);
%    for ( ; neighbor != neighbor_end; ++neighbor) { 
%        User\texttt{Particle}::pair\_t pairdata = *neighbor; 
%        int n = pairdata.first; 	           // local index of next neighbor 
%        double sep2 = pairdata.second;         // distance^2 between particles 
%        myBunch.R[i] += myBunch.vel[n] * 2.0;  // any calculation c<;>uld be here 
%    }
%}
%myPtcls.update();
%\end{codeln}
328
%On line 2, we declare the iterators used to obtain the beginning and end of the list of particle neighbors. These iterators are set in the call to \texttt{getpairlist} on line 4. Of particular note is how the data for each neighbor is obtained from the iterator, as done in lines 6-8: dereferencing the iterator returns an object of type \texttt{IpplParticleBase:: pair\_t}, which stores a pair of values the local index of the neighboring particle, and the square of the distance
gsell's avatar
gsell committed
329 330 331 332 333 334 335 336 337 338 339 340 341 342 343 344 345 346 347 348 349 350 351 352 353 354 355 356 357 358 359 360 361 362 363 364 365 366 367 368 369 370 371 372 373 374 375 376 377 378 379
%between the particles. The index and separation are obtained by accessing the member data elements \texttt{first} and \texttt{second}, respectively, of the \texttt{pair\_t} object.


\subsection{\texttt{Particle} - \texttt{Particle} Interactions}
\label{ppsection}
Efficient particle - particle interactions can be achieved by use of \texttt{PairBuilder} objects. The basic usage is as follows:

\begin{smallcode}
PairBuilder< Bunch<ParticleLayout_t> > PB(myBunch);
PB.for_each(PairCondition(), PairFunctor());
\end{smallcode}
This will call PairFunctor for each pair of particles that fulfills the PairCondition. There are three PairBuilders available: \texttt{HashPairBuilder}, \texttt{BasicBairBuilder} and \texttt{SortingPairBuilder}. HashPairBuilder should be used in most cases since it has the best time complexity. There are three predefined \texttt{PairConditions}. \texttt{TrueCondition} is always true and can be used to iterate over all pairs, \texttt{RadiusCondition} is true for each pair that is closer than a given interaction radius and \texttt{BoxCondition} is true for each pair where one particle lies inside a bounding box around the other particle. The following code shows how to iterate over all pairs within a given interaction radius:

\begin{smallcode}
struct PairFunctor{
  void operator()(std::size_t i, std::size_t j, Bunch<ParticleLayout_t> &P) const
  {
    //some interaction involving particles i and j
  }
};


HashPairBuilder< Bunch<ParticleLayout_t> > HPB(myBunch);
HPB.for_each(RadiusCondition<double, Dim>(interaction_radius), PairFunctor());
\end{smallcode}

To correctly generate all pairs in a multi process simulation a caching strategy has to be chosen so each process also has the required ghost particles. To achieve this for the example given one would call 

\begin{smallcode}
PL->setAllCacheDimensions(interaction_radius);
PL->enableCaching();
\end{smallcode}
with \texttt{PL} being the \texttt{ParticleSpatialLayout} of the bunch.

It is also possible to write custom pair conditions. These have to provide a \texttt{operator()} that takes two vectors and returns a bool, when the functor is used, it is passed two vectors that represent the particle positions. Additionally pair conditions have to provide a \texttt{getRange} function that takes an integer as input and returns the ``radius'' along that dimensions for which the pair condition can be true. In other words: for two particle positions \texttt{a} and \texttt{b} for which the pair condition returns true, the condition \texttt{|a[i] - b[i]| <= getRange(i)} must hold.

\subsection{\texttt{Particle} - \texttt{Field} Interactions}

%\section{\texttt{Particle} - \texttt{Field} Interactions}

Many particle-based simulation methods, including "particle-in-cell" (PIC) simulations, rely on the ability of particles to interact with field quantities. For instance, in particle-based accelerator (plasma) simulations, you typically track the motions of charged plasma particles in a combination of externally applied and self-generated electromagnetic fields. In a \ippl application, such fields might be stored as \texttt{Field} objects of type \texttt{Vektor} existing on a pre-defined mesh. Particles moving through this mesh must be able to "gather" the current value of a \texttt{Field} to their exact positions. Additionally, in order to compute the values of self-generated fields, the particles must be able to "scatter" the value of an attribute onto nearby mesh points, producing a \texttt{Field}. These gather/scatter operations are done using a set of \ippl interpolation methods.

\ippl provides a hierarchy of interpolation classes, each derived from the base class \texttt{Interpolate} and each containing the basic \texttt{gather} and \texttt{scatter} functions. The \texttt{gather} method allows you to gather one or more specified \texttt{Field}s into an equal number of \texttt{ParticleAttribs}. Similarly, \texttt{scatter} will accumulate one or more \texttt{ParticleAttribs} on to an equal number of \texttt{Field} objects. An example of how to scatter the particle density to a \texttt{Field} is shown below.
\begin{smallcode}
InterpolateNGP<Dim> mylnterpolater(myBunch);           // create NGP interpolater 
Field<double,Dim> ptcl_density(myfieldlayout);         // create density field 
myInterpolator.scatter(myBunch.density,ptcl_density);  // do scattet 
\end{smallcode}
The various classes derived from \texttt{Interpolate} implement these \texttt{gather} and \texttt{scatter} methods using different well-known interpolation schemes, such as nearest grid point (NGP), linear interpolation, and the subtracted-dipole scheme (SUDS). You may use these provided classes as a template for deriving new classes from \texttt{Interpolate} that implement other interpolation schemes of interest.

In case of the CIC Interpolation and non-cyclic boundary condition, care has to be taken to not place particles in the outer half of boundary cells. Otherwise values will be scattered out of the grid and be irretrivable.