8  Particles

Particles in IPPL are treated as a particle container, which stores particle attributes (such as charge, position, and velocity). Furthermore, it contains layout information, which tells each rank which particles they own according to domain decomposition.

This is the Lagrangian half of IPPL’s particle-mesh model.

8.1 Main classes

Class Role
ParticleBase Base class for user particle containers.
ParticleAttrib Typed particle attribute stored in a Kokkos view.
ParticleAttribBase Common attribute interface used by the container.
ParticleLayout Base layout interface.
ParticleSpatialLayout Spatial decomposition.
ParticleSpatialOverlapLayout Spatial layout with overlap regions.
ParticleBC Particle boundary condition support.

8.2 Particle anatomy

classDiagram
  class ParticleBase {
    +R
    +ID
    +create(nLocal)
    +globalCreate(nTotal)
    +destroy(invalid, destroyNum)
    +update()
    +addAttribute(attribute)
    +forAllAttributes(functor)
    +getLocalNum()
    +getTotalNum()
  }
  class ParticleAttrib {
    +getView()
    +getHostMirror()
    +create(n)
    +destroy(...)
    +scatter(field, positions, policy, hash)
    +gather(field, positions, addToAttribute)
    +sum()
    +min()
    +max()
    +prod()
  }
  class ParticleLayout {
    +setParticleBC(bc)
    +applyBC(R, region)
  }
  class ParticleSpatialLayout {
    +update(container)
    +updateLayout(fieldLayout, mesh)
    +locateParticles(...)
    +getRegionLayout()
  }

  ParticleSpatialLayout --|> ParticleLayout
  ParticleBase --> ParticleSpatialLayout
  ParticleBase --> ParticleAttrib

8.3 User model

A user defines a particle container by deriving from ParticleBase and adding attributes such as position, velocity, mass, charge, or application-specific quantity. The layout determines which MPI rank owns each particle and how particles move when positions change.

The programming interface clearly defines the role of each attribute and method:

  • R stores particle positions.
  • ID stores particle IDs when IDs are enabled.
  • user attributes are registered with addAttribute.
  • getLocalNum() returns the rank-local particle count.
  • getTotalNum() returns the global particle count.
  • update() sends particles to the correct rank after positions change (which could result in particles leaving the local domain of the owning rank).

ParticleBase automatically registers the inherited position attribute R. If particle IDs are enabled through the IDProperties template parameter, it also registers ID. User attributes must be registered in the derived class constructor with addAttribute.

8.4 Minimal derived container

template <class PLayout>
class ChargedParticles : public ippl::ParticleBase<PLayout> {
public:
    using Base = ippl::ParticleBase<PLayout>;
    using Vector_t = typename PLayout::vector_type;

    ippl::ParticleAttrib<double> qm;
    typename Base::particle_position_type P;
    typename Base::particle_position_type E;

    ChargedParticles(PLayout& layout)
        : Base(layout) {
        this->addAttribute(qm);
        this->addAttribute(P);
        this->addAttribute(E);
        this->setParticleBC(ippl::BC::PERIODIC);
    }
};

This is the same pattern used by test/particle/PICnd.cpp: charge-to-mass ratio, velocity, and gathered electric field are particle attributes; the inherited R position attribute drives migration and interpolation.

8.5 Attribute storage and expressions

ParticleAttrib<T, Properties...> stores one value per local particle in a Kokkos view. Common types include scalar attributes such as charge and vector attributes such as velocity:

ippl::ParticleAttrib<double> q;
typename Base::particle_position_type velocity;

Useful operations:

Operation Example
scalar assignment bunch.qm = totalCharge / totalParticles;
expression assignment bunch.P = bunch.P + dt * bunch.qm * bunch.E;
device view access auto view = bunch.qm.getView();
host mirror auto host = bunch.R.getHostMirror();
reduction double total = bunch.qm.sum();
field deposit (particle-mesh) scatter(bunch.qm, rho, bunch.R);
field gather (particle-mesh) gather(bunch.E, efield, bunch.R);

To initialize an attribute on host (CPU), one should use a host mirror followed by Kokkos::deep_copy to copy back to the device view. Then call update() so the rank ownership of particles according to domain decomposition is enforced.

For example, to initialize particle positions:

auto R_host = bunch.R.getHostMirror();
for (size_t i = 0; i < bunch.getLocalNum(); ++i) {
    R_host(i) = initialPosition(i);
}
Kokkos::deep_copy(bunch.R.getView(), R_host);
bunch.update();

8.6 Creation and initialization

For the particle container creation, a ParticleSpatialLayout is needed, which uses the same domain decomposition description as for the fields. Therefore, a FieldLayout of an IPPL mesh is needed first (to obtain the domain decomposition, which is fully grid-based).

Object Role in particle setup
NDIndex<Dim> Global grid extent (number of cells per axis).
UniformCartesian<T, Dim> (mesh) Physical spacing, origin, and cell geometry.
FieldLayout<Dim> (fieldLayout) MPI decomposition: which rank owns which part of the domain.
ParticleSpatialLayout<T, Dim> Maps particle positions to rank-owned regions using fieldLayout and mesh.

End-to-end setup (same pattern as test/particle/PICnd.cpp):

constexpr unsigned Dim = 3;

// 1) Domain and mesh to obtain FieldLayout
ippl::Vector<int, Dim> nr = {64, 64, 64};
ippl::NDIndex<Dim> domain;
for (unsigned d = 0; d < Dim; ++d) {
    domain[d] = ippl::Index(nr[d]);
}

ippl::Vector<double, Dim> origin = {0.0, 0.0, 0.0};
ippl::Vector<double, Dim> hr = {1.0 / nr[0], 1.0 / nr[1], 1.0 / nr[2]};

using Mesh_t = ippl::UniformCartesian<double, Dim>;
Mesh_t mesh(domain, hr, origin);

// we parallely decompose all dimensions
std::array<bool, Dim> isParallel;
isParallel.fill(true);

using FieldLayout_t = ippl::FieldLayout<Dim>;
FieldLayout_t fieldLayout(MPI_COMM_WORLD, domain, isParallel);

// 2) Particle layout and container
using PLayout_t = ippl::ParticleSpatialLayout<double, Dim>;
PLayout_t particleLayout(fieldLayout, mesh);

ChargedParticles<PLayout_t> bunch(particleLayout);

// 3) Create local particles and set attributes
const size_t totalParticles = 10000;
size_t nloc = totalParticles / ippl::Comm->size();
bunch.create(nloc);
bunch.qm = totalCharge / totalParticles;
bunch.P  = 0.0;

// 4) Enforce ownership after positions are set
bunch.update();

globalCreate(totalParticles) is an alternative to create(nloc), which requires hand-splitting nloc. If using globalCreate(totalParticles), every rank calls it with the same global count, and IPPL distributes particles across ranks.

8.6.1 Collective operations

create, globalCreate, createWithID, destroy, and update are collective in the MPI sense: every rank in the communicator must execute the call, even if local arguments differ (for example each rank passes its own nloc).

Do not guard these calls so that only some ranks run them:

// Wrong: rank 0 skips create(); other ranks hang or corrupt state
if (ippl::Comm->rank() == 0) {
    bunch.create(nloc);
}

// Correct: every rank calls create; nloc may differ per rank
bunch.create(nloc);

The same rule applies after if branches: if one code path calls update(), every rank must take a path that also calls update() (or none must). Avoid return on one rank before a collective while others continue.

8.7 Particle boundary conditions

Particle boundary conditions are selected with the ippl::BC enum:

BC Behavior
PERIODIC Wraps coordinates back into the domain.
REFLECTIVE Mirrors coordinates back into the domain.
SINK Clips coordinates to the boundary.
NO Leaves coordinates unchanged.

Use one boundary condition on all faces:

bunch.setParticleBC(ippl::BC::PERIODIC);

or provide one value per face:

typename Bunch::bc_container_type bcs;
bcs.fill(ippl::BC::PERIODIC);
bunch.setParticleBC(bcs);

The spatial layout applies particle boundary conditions at the start of update(), before migration of particles to appropriate ranks. Unit tests cover upper and lower faces for periodic, reflective, sink, and no-boundary behavior.

8.8 Particle update (spatial migration)

ParticleSpatialLayout maps particle positions to rank-owned spatial regions derived from a FieldLayout and a mesh. The update algorithm is:

  1. Apply particle boundary conditions to R.
  2. Locate each particle in the current rank, neighboring ranks, or the global region list.
  3. Count sends per destination rank.
  4. Pack and send particles that left the current rank.
  5. Delete invalid local particles.
  6. Receive incoming particles and unpack all registered attributes.

flowchart LR
  R["Updated positions R"]
  BC["applyBC"]
  Locate["locateParticles"]
  Pack["pack/send"]
  Delete["internalDestroy"]
  Recv["recv/unpack"]
  Valid["Rank-local valid particles"]

  R --> BC --> Locate --> Pack --> Delete --> Recv --> Valid

The ParticleSendRecv unit test verifies that after migration every remaining local particle has an expected rank equal to ippl::Comm->rank(), and the global particle count is conserved.

8.9 Developer model

Particle data movement is a correctness-critical layer. Changes to particle layouts must be documented with:

  • communication pattern,
  • memory-space assumptions,
  • unit tests covering multiple dimensions and execution spaces.

The test suite instantiates particle tests across floating point precision, dimensions (1D to 3D), and execution spaces. When adding a particle feature, prefer tests that cover both host and device memory-space paths and more than one dimension.