Particles

Parthenon provides a data framework for particle methods that allows for Kokkos accelerated parallel dispatch of compute operations. Particle memory is allocated as a separate pool for each variable in the particle.

Swarms

A Swarm contains all the particle data for all particles of a given species. It owns a set of ParticleVariables, one for each value of each particle. For example, the spatial positions x, y, and z of the particles in a swarm are three separate ParticleVariables. ParticleVariables can be either Real- or int-valued, which is specified by the metadata values Metadata::Real and Metadata::Integer. ParticleVariables should also contain the Metadata::Particle flag. By default, ParticleVariables provide one scalar quantity per particle, but up to 2D data per particle is currently supported, by passing std::vector<int>{N1, N2} as the second argument to the ParticleVariable Metadata. All Swarms by default contain x, y, and z ParticleVariables; additional fields can be added as:

Swarm.Add(name, metadata)

For a given species, each MeshBlock contains its own Swarm that holds the particles of that species that are spatially contained by that MeshBlock. The MeshBlock is pointed to by Swarm::pmy_block.

The Swarm is a host-side object, but some of its data members are required for device- side compution. To access this data, a SwarmDeviceContext object is created via Swarm::GetDeviceContext(). This object can then be passed by copy into Kokkos lambdas. Hereafter we refer to it as swarm_d.

To add particles to a Swarm, one calls

ParArray1D<bool> new_particles_mask = swarm->AddEmptyParticles(num_to_add, new_indices)

This call automatically resizes the memory pools as necessary and returns a ParArray1D<bool> mask indicating which indices in the ParticleVariables are newly available. new_indices is a reference to a ParArrayND<int> of size num_to_add which contains the indices of each newly added particle.

To remove particles from a Swarm, one first calls

swarm_d.MarkParticleForRemoval(index_to_remove)

inside device code. This only indicates that this particle should be removed from the pool, it does not actually update any data. To remove all particles so marked, one then calls

swarm.RemoveMarkedParticles()

in host code. This updates the swarm such that the marked particles are seen as free slots in the memory pool.

Parallel Dispatch

Parallel computations on particle data can be performed with the usual MeshBlock par_for calls. Typically one loops over the entire range of active indices and uses a mask variable to only perform computations on currently active particles:

auto &x = swarm.Get("x").Get();
swarm.pmy_block->par_for("Simple loop", 0, swarm.GetMaxActiveIndex(),
  KOKKOS_LAMBDA(const int n) {
    if (swarm_d.IsActive(n)) {
      x(n) += 1.0;
    }
  });

Sorting

By default, particles are stored in per-meshblock pools of memory. However, one frequently wants convenient access to all the particles in each computational cell separately. To facilitate this, the Swarm provides the method SortParticlesByCell (and the SwarmContainer provides the matching task SortParticlesByCell). Calling this function populates internal data structures that map from per-cell indices to the per-meshblock data array. These are accessed by the SwarmDeviceContext member functions GetParticleCountPerCell and GetFullIndex. See examples/particles for example usage.

Defragmenting

Because one typically loops over particles from 0 to max_active_index, if only a small fraction of particles in that range are active, significant effort will be wasted. To clean up these situations, Swarm provides a Defrag method which, when called, will copy all active particles to be contiguous starting from the 0 index. Defrag is not fully parallelized so should be called only sparingly.

SwarmContainer

A SwarmContainer contains a set of related Swarms, such as the different stages used by a higher order time integrator. This feature is currently not exercised in detail.

particles Example

An example showing how to create a Parthenon application that defines a Swarm and creates, destroys, and transports particles is available in parthenon/examples/particles.

Communication

Communication of particles across MeshBlocks, including across MPI processors, is supported. Particle communication is currently handled via paired asynchronous/synchronous tasking regions on each MPI processor. The asynchronous tasks include transporting particles and SwarmContainer::Send and SwarmContainer::Receive calls. The synchronous task checks every MeshBlock on that MPI processor for whether the Swarms are finished transporting. This set of tasks must be repeated in the driver’s evolution function until all particles are completed. See the particles example for further details. Note that this pattern is blocking, and may be replaced in the future.

AMR is currently not supported, but support will be added in the future.

Variable Packing

Similarly to grid variables, particle swarms support ParticleVariable packing, by the function Swarm::PackVariables. This also supports FlatIdx for indexing; see the particle_leapfrog example for usage.

Boundary conditions

Particle boundary conditions are not applied in separate kernel calls; instead, inherited classes containing boundary condition functions for updating particles or removing them when they are in boundary regions are allocated depending on the boundary flags specified in the input file. Currently, outflow and periodic boundaries are supported natively. User-specified boundary conditions must be set by specifying the “user” flag in the input parameter file and then updating the appropriate Swarm::bounds array entries to factory functions that allocate device-side boundary condition objects. An example is given in the particles example when ix1 and ox1 are set to user in the input parameter file.