Particles

chombo-discharge stores computational particles in a Struct-of-Arrays (SoA) layout. The source code for the particle functionality resides in $DISCHARGE_HOME/Source/Particle. Particle support contains the following basic features:

• Particle-mesh operations, i.e., deposition and interpolation of particle variables to/from the mesh.

• Particle distribution and remapping with MPI.

• Rudimentary particle output to H5Part files.

Particle support is templated on a small user-defined payload struct, so that users can define new particle types that contain a desired set of variables. The particle position, weight, and bookkeeping IDs are always owned by the container; the payload only adds the extra per-particle variables.

ParticleSoA

ParticleSoA<P, Traits> is the per-patch Struct-of-Arrays leaf that holds the particles of one grid patch. Rather than storing an array of particle objects, it stores one contiguous column per variable. The columns fall into two groups:

• Container-owned columns, present for every particle type: the position (SpaceDim scalar columns), the weight, and the bookkeeping particleID/rankID. These are not declared by the user.

• Payload columns, supplied by the user as a plain struct P whose data members become extra columns.

The payload is described by a ParticleTraits<P> specialization that lists the payload columns as a tuple of member pointers (and, optionally, an h5Columns subset that restricts which columns are written to HDF5 checkpoints). A representative payload (the tracer-particle velocity + Runge-Kutta scratch) looks like

/**
 * @brief SoA payload for tracer particles: interpolated velocity plus Runge-Kutta stage scratch.
 * @details Per-component scalar columns (one member per spatial dimension; the z-component exists only when
 * CH_SPACEDIM == 3). v_x/v_y/v_z hold the interpolated velocity (required by TracerParticleSolver). xk is the
 * position snapshot x^k; k1/k2/k3 hold the Runge-Kutta stage velocities. Weight and position are
 * container-owned columns, NOT payload members.
 */
struct TracerParticle
{
  ParticleReal v_x = 0.0; ///< Interpolated velocity, x-component.
  ParticleReal v_y = 0.0; ///< Interpolated velocity, y-component.
#if CH_SPACEDIM == 3
  ParticleReal v_z = 0.0; ///< Interpolated velocity, z-component.
#endif

  double xk_x = 0.0; ///< Position snapshot x^k, x-component (double: position-like).
  double xk_y = 0.0; ///< Position snapshot x^k, y-component (double: position-like).
#if CH_SPACEDIM == 3
  double xk_z = 0.0; ///< Position snapshot x^k, z-component (double: position-like).
#endif

  ParticleReal k1_x = 0.0; ///< Runge-Kutta stage k1, x-component.
  ParticleReal k1_y = 0.0; ///< Runge-Kutta stage k1, y-component.
#if CH_SPACEDIM == 3
  ParticleReal k1_z = 0.0; ///< Runge-Kutta stage k1, z-component.
#endif

  ParticleReal k2_x = 0.0; ///< Runge-Kutta stage k2, x-component.
  ParticleReal k2_y = 0.0; ///< Runge-Kutta stage k2, y-component.
#if CH_SPACEDIM == 3
  ParticleReal k2_z = 0.0; ///< Runge-Kutta stage k2, z-component.
#endif

  ParticleReal k3_x = 0.0; ///< Runge-Kutta stage k3, x-component.
  ParticleReal k3_y = 0.0; ///< Runge-Kutta stage k3, y-component.
#if CH_SPACEDIM == 3
  ParticleReal k3_z = 0.0; ///< Runge-Kutta stage k3, z-component.
#endif
};

/**
 * @brief ParticleTraits specialization listing the payload columns of TracerParticle.
 */
template <>
struct ParticleTraits<TracerParticle>
{
  /// @brief Payload column member-pointers, listed one field group at a time (SpaceDim entries per group).
  static constexpr auto columns = std::make_tuple(
    D_DECL(&TracerParticle::v_x, &TracerParticle::v_y, &TracerParticle::v_z),
    D_DECL(&TracerParticle::xk_x, &TracerParticle::xk_y, &TracerParticle::xk_z),
    D_DECL(&TracerParticle::k1_x, &TracerParticle::k1_y, &TracerParticle::k1_z),
    D_DECL(&TracerParticle::k2_x, &TracerParticle::k2_y, &TracerParticle::k2_z),
    D_DECL(&TracerParticle::k3_x, &TracerParticle::k3_y, &TracerParticle::k3_z));

  /**
   * @brief Payload columns written to H5Part visualization files: the interpolated velocity vector (id,
   * position and weight are written automatically; the position snapshot and Runge-Kutta stage
   * columns are internal and not plotted).
   */
  static constexpr auto h5PartColumns = std::make_tuple(
    H5Part::Vector{"v", D_DECL(&TracerParticle::v_x, &TracerParticle::v_y, &TracerParticle::v_z)});
};

Per-component vectors are declared as individual scalar columns (there is no RealVect column type); the D_DECL macro expands to the SpaceDim components. The empty payload NoPayload is provided for particles that only need position and weight.

The most common ParticleSoA<P, Traits> member functions are:

• size() – the number of particles in the leaf.

• position(i) / setPosition(i, x) – get/set the position of particle i.

• weight(i) – reference to the weight of particle i.

• get<&P::member>(i) – reference to a payload column entry of particle i.

• column<&P::member>() – raw base pointer to a whole payload column.

• append(pos, weight) or append(pos, weight, payload) – add a particle.

• remove(i) – delete particle i (O(1) swap-and-pop; does not preserve order).

• gather(i) – assemble particle i’s payload back into a P value.

• swap / catenate – swap arenas with, or move all particles from, another leaf.

A typical loop over the particles in a leaf uses an integer index and column access:

ParticleSoA<MyPayload> leaf;

for (std::size_t i = 0; i < leaf.size(); i++) {
   const RealVect x = leaf.position(i);
   const Real     w = leaf.weight(i);

   // Access a payload column entry:
   Real& vx = leaf.get<&MyPayload::vx>(i);
}

Tip

The ParticleSoA C++ API is found at https://chombo-discharge.github.io/chombo-discharge/doxygen/html/classParticleSoA.html.

Checkpoint and HDF5 export

The container-owned position and weight columns are always written to HDF5 checkpoint files. For particles that do not need to checkpoint all payload columns, the ParticleTraits<P> specialization may declare an h5Columns tuple that lists the payload-column subset to export, which reduces the checkpoint file size.

ParticleContainer

The ParticleContainer<P, Traits> is a template class that

1. Stores computational particles of type P over an AMR hierarchy (one ParticleSoA<P, Traits> leaf per grid patch).

2. Provides infrastructure for remapping particles.

3. Provides functionality for getting the particles within a specified grid patch.

4. Provides functionality that is required during regrids.

5. Other types of functionality, like grouping particles into grid cells, and mask/halo particle extraction.

Data structures

ParticleSoA<P> leaves

At the lowest level the particles in one grid patch are stored in a ParticleSoA<P, Traits> (see ParticleSoA), which holds the particles column-by-column with no ordering unless the leaf has been cell-sorted.

AMRParticlesSoA<P>

On each grid level, ParticleContainer<P, Traits> stores the leaves in a LayoutData<ParticleSoA<P, Traits>> (one leaf per patch). The AMR view AMRParticlesSoA<P, Traits> is a vector of these per-level holders:

template <typename P, typename Traits>
using AMRParticlesSoA = Vector<RefCountedPtr<LayoutData<ParticleSoA<P, Traits>>>>;

Again, the Vector indicates the AMR level and the LayoutData is a distributed data holder that holds the leaves on each AMR level. AMRParticlesSoA<P, Traits> always lives within ParticleContainer<P, Traits> and is the class member that actually holds the particles.

Basic usage

Here, we give some examples of basic usage of ParticleContainer. For the full API, see the ParticleContainer C++ API https://chombo-discharge.github.io/chombo-discharge/doxygen/html/classParticleContainer.html.

Getting the particles

To get the per-level holders from a ParticleContainer<P, Traits> one can call getParticles():

/**
 * @brief The valid particles on all levels.
 * @return Mutable per-level holder vector.
 */
AMRParticlesSoA<P, Traits>&
getParticles()
{
  return m_particles;
}

Alternatively, one can fetch the distributed leaves of a specified grid level with operator[]:

int lvl;
ParticleContainer<P> myParticleContainer;

LayoutData<ParticleSoA<P>>& levelParticles = myParticleContainer[lvl];

Iterating over particles

To do something with the particles in a ParticleContainer<P, Traits>, one iterates over the grid levels and patches, gets the ParticleSoA leaf for each patch, and loops over the particle indices.

The code bit below shows a typical example of how the particles can be moved, and then remapped onto the correct grid patches and ranks if they fall off their original one.

ParticleContainer<P> myParticleContainer;

// Iterate over grid levels
for (int lvl = 0; lvl <= m_amr->getFinestLevel(); lvl++){

   // Get the grid on this level.
   const DisjointBoxLayout& dbl = m_amr->getGrids(myParticleContainer.getRealm())[lvl];

   // Iterate over grid patches on this level
   for (DataIterator dit(dbl); dit.ok(); ++dit){

      // Get the SoA leaf in the current patch.
      ParticleSoA<P>& leaf = myParticleContainer[lvl][dit()];

      // Iterate over the particles in the current patch.
      for (std::size_t i = 0; i < leaf.size(); i++){

         // Move the particle
         leaf.setPosition(i, ...);
      }
   }
}

// Remap particles onto new patches and ranks (they may have moved off their original ones)
myParticleContainer.remap();

Sorting particles

Sorting by cell

The particles in a leaf can be sorted by cell by calling ParticleContainer<P>::organizeParticlesByCell():

ParticleContainer<P> myParticleContainer;

myParticleContainer.organizeParticlesByCell();

Internally this performs a counting sort of each leaf’s columns into Fortran (x-fastest) cell order and builds a compressed-sparse-row (CSR) offset array. After the sort, the particles of cell c occupy the contiguous index range [cellStart(c), cellStart(c+1)). Unlike the AoS-era per-cell containers, the cell-sort does not move the particles into a separate structure – it merely reorders the SoA columns in place.

The relevant ParticleSoA<P, Traits> query functions are isSorted(), numCells(), cellStart(c), particlesInCell(c), and cellRange(c) (which returns the {begin, end} index pair). Iteration over cell-sorted particles visits the cells in Fortran order (matching the CSR cell index) and then the contiguous index range of each cell:

ParticleContainer<P> myParticleContainer;
myParticleContainer.organizeParticlesByCell();

// Iterate over all AMR levels
for (int lvl = 0; lvl <= m_amr->getFinestLevel(); lvl++){

   const DisjointBoxLayout& dbl = m_amr->getGrids(myParticleContainer.getRealm())[lvl];

   // Iterate over grid patches on this level
   for (DataIterator dit(dbl); dit.ok(); ++dit){

      const Box       cellBox = dbl[dit()];
      ParticleSoA<P>& leaf    = myParticleContainer[lvl][dit()];

      // Visit cells in Fortran order, matching the CSR cell index.
      std::size_t cellIndex = 0;
      for (BoxIterator bit(cellBox); bit.ok(); ++bit, ++cellIndex){
         const std::pair<std::size_t, std::size_t> range = leaf.cellRange(cellIndex);

         for (std::size_t i = range.first; i < range.second; i++){
            // Do something with particle i in cell bit().
         }
      }
   }
}

Sorting by patch

To return to patch-ordered particles:

ParticleContainer<P> myParticleContainer;

myParticleContainer.organizeParticlesByPatch();

Important

The CSR cell ranges are only valid while the leaf is sorted. append and remove invalidate the sort, and remap() requires patch-ordered particles.

Allocating particles

AmrMesh has a simple function for allocating a ParticleContainer<P, Traits>:

/**
 * @brief Allocate a struct-of-arrays particle container on the given realm.
 * @param[out] a_container SoA particle container to be allocated.
 * @param[in]  a_realm     Realm on which the particles will be allocated.
 * @tparam P      Payload type of the container.
 * @tparam Traits Column descriptor for the payload.
 */
template <typename P, typename Traits>
void
allocate(ParticleContainer<P, Traits>& a_container, const std::string& a_realm) const;

which will allocate the container on the realm a_realm.

Particle mapping

Particles that move off their original grid patch must be remapped in order to ensure that they are assigned to the correct grid. The remapping function for ParticleContainer<P, Traits> is

/**
 * @brief Redistribute every valid particle to the patch/level/rank that owns its cell.
 * @details Pool model: collect all valid particles, map each to its destination via the
 * finest-level tile that contains it (LevelTiles), append same-rank movers into the destination
 * leaf and scatter cross-rank movers with MPI. Particles whose cell is not owned by any tile on
 * any level (off-domain) are dropped and counted in getNumberOfOutcastParticles*(). Invalidates
 * any cell-sort. Particle ids are preserved; rankID is set to the new owning rank.
 */
void
remap();

This is simply used as follows:

ParticleContainer<P> myParticles;

myParticles.remap();

During remapping, the following steps are performed for each MPI rank:

1. Collect all valid particles from this rank.

2. Map each particle to its destination in the AMR hierarchy (level, grid index, and owning MPI rank).

3. Append the particles that stay on this rank to their destination leaf.

4. Scatter the particles that move to another rank with MPI.

5. Assign the scattered particles on each receiving rank.

Particles whose cell is owned by no patch on any level (off-domain) are dropped and counted.

The point-to-patch mapping in step 2 maps the position to a min_block_size tile by integer division and looks it up in a per-level hash map. This is \(\mathcal{O}(1)\) and works also when the grids contain variable-sized (anisotropic) boxes (see Mesh generation), because every box is registered under each min_block_size tile it covers. The same mapping is exposed on the Realm class through Realm::getLevelAndBox (see Patch lookup (tile hash grid)) for non-particle users; the container aliases the realm’s hash grid rather than building its own.

Regridding

As with mesh data, ParticleContainer<P, Traits> requires storing the old-grid particles before assigning them on the new grids. This is done as follows:

1. Before creating the new grids, each MPI rank caches its current particles by calling

   /**
    * @brief Cache the current valid particles ahead of a regrid.
    * @details Snapshots the current per-level holders (and their layout) so that regrid() can
    * rebuild over a new layout and redistribute the cached particles. Call exactly once before
    * regrid().
    */
   void
   preRegrid();
   

   This snapshots the particles off their current grids.

2. When ParticleContainer<P, Traits> regrids, the cached particles are redistributed onto the new layout by calling the regrid function:

   /**
    * @brief Rebuild over a new AMR layout and redistribute the preRegrid()-cached particles onto it.
    * @details Re-allocates the per-level holders and tile maps over the new grids (all levels), then
    * runs the same pool->map->scatter as remap() but sourced from the cache. Positions are unchanged;
    * each cached particle is routed to the patch/level/rank owning its cell on the NEW layout, with
    * ids preserved. Off-domain particles are counted in getNumberOfOutcastParticles*(). Must follow a
    * preRegrid().
    * @param[in] a_grids          New per-level DisjointBoxLayouts.
    * @param[in] a_domains        New per-level problem domains.
    * @param[in] a_dx             New per-level (isotropic) grid spacing.
    * @param[in] a_refRat         New per-level refinement ratios.
    * @param[in] a_minBlockSize Grid blocking factor (tile size) for the new grids.
    * @param[in] a_levelTiles     The owning Realm's per-level tile->box maps for the new grids, aliased
    *                             (never built locally). Must correspond to a_grids.
    * @param[in] a_newFinestLevel New finest AMR level index.
    */
   void
   regrid(const Vector<DisjointBoxLayout>&         a_grids,
          const Vector<ProblemDomain>&             a_domains,
          const Vector<Real>&                      a_dx,
          const Vector<int>&                       a_refRat,
          const int                                a_minBlockSize,
          const Vector<RefCountedPtr<LevelTiles>>& a_levelTiles,
          const int                                a_newFinestLevel);
   

Warning

One must call preRegrid before the regrid. Failure to do so will lead to loss of all particles.

Masked particles

ParticleContainer<P, Traits> also supports the concept of masked particles, where one can fetch a subset of particles that live only in specified grid cells. Typically, this “specified region” is the refinement boundary, but the functionality is generic and might prove useful also in other cases. This functionality is unlikely to be used directly by users of chombo-discharge, but it is nonetheless fruitful to understand the concept in order to more easily fathom how deposition across refinement boundaries proceeds.

When masked particles are used, the user provides a boolean mask over the AMR hierarchy and obtains the subset of particles that live in regions where the mask evaluates to true. This functionality is for example used for some of the particle deposition methods in chombo-discharge where we deposit particles that live near the refinement boundary with special deposition functions.

To fill the masked particles, ParticleContainer<P, Traits> has member functions for copying the particles into internal data containers which the user can later fetch. The function signature for this is

void
copyMaskParticles(const Vector<RefCountedPtr<LevelData<BaseFab<bool>>>>& a_mask);

The argument a_mask holds a bool at each cell in the AMR hierarchy. Particles that live in cells where a_mask is true will be copied to an internal holder which can be retrieved through

/**
 * @brief The halo/mask particles on all levels.
 * @return Mutable per-level holder vector.
 */
AMRParticlesSoA<P, Traits>&
getMaskParticles()
{
  return m_maskParticles;
}

In the above functions the mask particles are copied, and the original particles are left untouched. After the user is done with the particles, they should be released through

/**
 * @brief Drop all mask/halo particles on every level (keeps capacity).
 */
void
clearMaskParticles()

An example pseudocode for working with masked particles is given below:

AmrMask myMask;
ParticleContainer<P> myParticles;

// Copy mask particles
myParticles.copyMaskParticles(myMask);

// Do something with the mask particles
AMRParticlesSoA<P>& maskParticles = myParticles.getMaskParticles();

// Release the mask particles
myParticles.clearMaskParticles();

Ghost particles

ParticleContainer<P, Traits> can populate each patch with a halo of ghost particles: transient, non-owned copies of particles that live in an adjacent region and overlap the patch’s ghosted/grown box. Ghosts let per-patch code (e.g. short-range interactions or custom deposition) see the relevant neighbour particles in place. A copy can come from three directions relative to the receiving patch, and the direction is recorded on each particle:

• Same level – a neighbouring patch on the same AMR level (GhostType::SameLevel).

• Coarse-to-fine – the next-coarser level (GhostType::Coarse).

• Fine-to-coarse – the next-finer level (GhostType::Fine).

The copies are appended into the same ParticleSoA leaf as the valid particles and marked in place through the mandatory 1-byte GhostType column, queried with ParticleSoA::isGhost(i) (and ghost(i)). Ghosts are transient: they exist only between a fillGhostParticles and the next clearGhostParticles / remap / regrid, and are never routed or counted as owned particles (getNumberOfValidParticlesLocal() excludes them).

Fig. 11 illustrates the three scatter directions for a single fine patch and its coarse/same-level neighbours.

Fig. 11 Ghost-particle scatter directions. A patch collects same-level, coarse-to-fine, and fine-to-coarse ghost copies within a user-selected ghost width.

Ghost masks

The scatter is driven by prebuilt masks rather than a per-fill geometric search. For each source cell a mask lists the destination boxes (grid index and receiving rank) whose ghost region the cell reaches, stored as a compact per-box CSR table (ParticleGhostMask). There is one mask per direction, and the collection type is

using AMRParticleGhostMask = Vector<RefCountedPtr<LayoutData<ParticleGhostMask>>>;

The masks live on the Realm and are rebuilt at every regrid. Downstream code must first register the ghost width(s) it needs; only registered widths are built. The ghost width is a minimum, measured in destination cells:

void
registerParticleGhostMask(const std::string& a_realm, const int a_width);

and the per-direction masks for a registered width are then fetched with

const AMRParticleGhostMask& same = amr->getParticleGhostMask(realm, width);
const AMRParticleGhostMask& c2f  = amr->getParticleGhostMaskCoarToFine(realm, width);
const AMRParticleGhostMask& f2c  = amr->getParticleGhostMaskFineToCoar(realm, width);

For a direction that should not scatter, pass AmrMesh::getTrivialParticleGhostMask() (an empty mask).

Note

Because coarse particles are stored per coarse cell but the coarse-to-fine ghost width is measured in fine cells, a whole boundary coarse cell would over-communicate by up to the refinement ratio. The coarse-to-fine mask therefore stores, per target, an acceptance box (the fine box’s ghost region), and the scatter keeps a coarse particle only if its fine cell lands in that box – pruning the transfer to the exact fine shell. Same-level and fine-to-coarse transfers are already exact and carry no acceptance box.

Filling ghosts

With the three masks in hand, the ghost particles are filled with

void
fillGhostParticles(const AMRParticleGhostMask& a_maskSame,
                   const AMRParticleGhostMask& a_maskCoarToFine,
                   const AMRParticleGhostMask& a_maskFineToCoar);

fillGhostParticles first clears any existing halo, then for every valid particle scatters a copy (same-rank directly, cross-rank via MPI) into each destination leaf listed by the masks, tagging it with the receiver-view GhostType and keeping the source’s particleID. A typical call, having registered width during registerRealms/setup, is

ParticleContainer<P>& particles = ...;

const AMRParticleGhostMask& same = amr->getParticleGhostMask(realm, width);
const AMRParticleGhostMask& c2f  = amr->getParticleGhostMaskCoarToFine(realm, width);
const AMRParticleGhostMask& f2c  = amr->getParticleGhostMaskFineToCoar(realm, width);

particles.fillGhostParticles(same, c2f, f2c);   // pass getTrivialParticleGhostMask() to skip a direction

The halo can be dropped at any time with

void
clearGhostParticles();

Tip

The full ghost/mask C++ API is documented at https://chombo-discharge.github.io/chombo-discharge/doxygen/html/classParticleContainer.html and https://chombo-discharge.github.io/chombo-discharge/doxygen/html/classParticleGhostMask.html.

Boundary interaction

ParticleContainer<P, Traits> is EB-agnostic and has no information about the embedded boundary and only partial information about the domain boundary. This means the following:

1. Particles remap just as if the embedded boundary was not there.

2. Particles that completely fall off the domain are deleted when calling the remapping function.

Interaction with the EB is done via the implicit function or discrete information, as well as modifications in the interpolation and deposition steps.

Signed distance function

When signed distance functions are used, one can always query how far a particle is from a boundary:

ParticleSoA<P>& leaf;
BaseIF          distanceFunction;

for (std::size_t i = 0; i < leaf.size(); i++){
   const RealVect pos = leaf.position(i);

   const Real distanceToBoundary = distanceFunction.value(pos);
}

If the particle is inside the EB then the signed distance function will be positive, and the particle can then be removed from the simulation. The distance function can also be used to detect collisions between particles and the EB. E.g, the intersection point can be computed and the particle can be deposited on the boundary, or bounced off it. See AmrMesh for details on how to obtain the distance function.

Domain edges

By default, the ParticleContainer remapping function will discard particles that fall outside of the domain. The user can also check if this happens by checking if the particle position is outside the computational domain:

const RealVect pos    = leaf.position(i);
const RealVect probLo = m_amr->getProbLo();
const RealVect probHi = m_amr->getProbHi();

bool outside = false;
for (int dir = 0; dir < SpaceDim; dir++) {
   if(pos[dir] < probLo[dir] || pos[dir] > probHi[dir]) {
      outside = true;
   }
}

Particle intersection

It is occasionally useful to catch particles that hit an EB or crossed a domain side. Provided that the payload stores the previous position of the particle (one scalar column per component), one can compute the intersection point between the particle trajectory and the EB or domain sides. Currently, AmrMesh supports two methods for computing this

• Using a bisection algorithm with a user-specified step.

• Using a ray-casting algorithm.

These algorithms differ in the sense that the bisection approach will check for a particle crossing between two positions \(\mathbf{x}_0\) and \(\mathbf{x}_1\) using a pre-defined tolerance. The ray-casting algorithm will check if the particle can move from \(\mathbf{x}_0\) towards \(\mathbf{x}_1\) by using a variable step along the particle trajectory. This step is selected from the signed distance from the particle position to the EB such that it uses a large step if the particle is far away from the EB. Conversely, if the particle is close to the EB a small step will be used.

The algorithms that intersect the particles are part of AmrMesh, and the ray-casting variant is called as follows:

void
intersectParticlesRaycastIF(
  ParticleContainer<P, Traits>&                                    a_activeParticles,
  ParticleContainer<P, Traits>&                                    a_ebParticles,
  ParticleContainer<P, Traits>&                                    a_domainParticles,
  const phase::which_phase                                         a_phase,
  const Real                                                       a_tolerance,
  const bool                                                       a_deleteParticles,
  const std::function<void(ParticleSoA<P, Traits>&, std::size_t)>& a_nonDeletionModifier) const noexcept;

The container-owned position holds the trajectory end point, while the start point is read from the payload columns selected by the OldPosition member-pointer pack (SpaceDim members, one per component). The intersected particles are put into the EB-intersected particles (a_ebParticles) and domain-intersected particles (a_domainParticles). The user can choose whether or not to remove intersected particles from a_activeParticles by adjusting a_deleteParticles. The final argument lets the user supply a callback (ParticleSoA&, std::size_t) that modifies particles that were intersected but not deleted (for example to flag the original particle via a payload column).

Both the bisection and ray-casting algorithm have weaknesses. The bisection algorithm requires a user-supplied step in order to operate efficiently, while the ray-casting algorithm is very slow when the particle is close to the EB and moves tangentially along it. Future versions of chombo-discharge will likely include more sophisticated algorithms.

Tip

AmrMesh also stores the implicit function on the mesh, which could also be used to resolve particle collisions with the EB/domain.

Particle-mesh

Particle-mesh operations are required when particles interact with the mesh and vice-versa. There are two main operations involved:

1. Deposition, where particle properties are transferred to the mesh.

2. Interpolation, where mesh properties are transferred to the particles.

Particle deposition

To deposit the particle weight on the mesh, the user can call AmrMesh::depositWeight:

void
depositWeight(EBAMRCellData&                      a_meshData,
              const std::string&                  a_realm,
              const phase::which_phase&           a_phase,
              const ParticleContainer<P, Traits>& a_particles,
              const DepositionType                a_depositionType,
              const CoarseFineDeposition          a_coarseFineDeposition,
              const bool                          a_forceIrregNGP);

The input arguments are the output mesh data holder (must have exactly one component), the realm and phase where the particles live, the SoA particle container (a_particles), the deposition method, the coarse-fine handling, and a flag that enforces nearest grid-point deposition in cut-cells. The last flag is motivated by the fact that some applications might require hard mass conservation, and the user can then ensure that mass is never deposited into covered grid cells.

To deposit a derived per-particle quantity (e.g. weight times a payload column), use AmrMesh::depositGathered with a gatherer callback, or AmrMesh::depositParticles<Members...> to deposit one or more payload columns directly. Surface (EB) deposition of the weight column onto an EBAMRIVData is available through an overload of AmrMesh::depositParticles:

void
depositParticles(EBAMRIVData&                        a_meshData,
                 const std::string&                  a_realm,
                 const phase::which_phase&           a_phase,
                 const ParticleContainer<P, Traits>& a_particles) const noexcept;

The input argument a_depositionType indicates the deposition method, while a_coarseFineDeposition selects deposition modifications near refinement boundaries. These are discussed below.

Base deposition

The base deposition scheme is specified by an enum DepositionType with valid values:

• DepositionType::NGP (Nearest grid-point).

• DepositionType::CIC (Cloud-In-Cell).

• DepositionType::TSC (Triangle-Shaped Cloud).

chombo-discharge supports all of the above methods, which can be combined with various types of modifications near refinement boundaries.

Coarse-fine deposition

The input argument a_coarseFineDeposition determines how deposition near the coarse-fine boundary is handled. Refinement boundaries introduce additional complications in the deposition scheme due to

1. Fine-grid particles whose deposition clouds hang over the refinement boundary and onto the coarse level.

2. Coarse-grid particles whose deposition clouds stick underneath the fine-level.

In addition, there can be complications near physical boundaries, such as domain or embedded boundaries.

Fig. 12 Sketch of deposition schemes near refinement boundaries and cut-cells.

chombo-discharge supports various ways of handling deposition across the refinement boundary. In all of these methods, the mass on the fine grid particles whose deposition clouds hang over the refinement boundaries is simply added to the coarse grid. The main modifications to the deposition scheme are performed for the coarse-grid particles that live around the refinement boundary (see Fig. 13). For the coarse-grid particles the following processes then occur:

Fig. 13 Example regions containing coarse-grid particles that deposit with custom deposition rules.

The following coarse-fine deposition methods are currently supported:

• CoarseFineDeposition::Interp This method permits the coarse-grid particles to deposit into the region underneath the fine grid. The deposited mass that falls underneath the fine grid is then interpolated from the coarse grid to the fine grid. For example, see the indicated coarse-grid particle cloud in the left panel Fig. 12. While this particle has a width given by the coarse-grid cell size, it will deposit into the coarse grid cells underneath the fine grid. The mass that ends up in these cells is interpolated to the fine grid, which in this case will inject mass into two layers of fine-grid cells.

• CoarseFineDeposition::Halo This method extracts the coarse-grid particles that live on the refinement boundary and deposit them with their original width on both the coarse and fine levels. This is done by first depositing the particles on the coarse level, and then transferring them to the fine level and redepositing them there with the original particle width. Taking the left panel in Fig. 12 as an example, the green particle will then deposit into the coarse-grid cell as well as the first layer of fine-grid cells.

• CoarseFineDeposition::HaloNGP Similar to CoarseFineDeposition::Halo discussed above, this method also extracts the coarse-grid particles on the coarse side of the refinement boundary. However, rather than using the original deposition scheme, these particles are deposited with an NGP scheme.

• CoarseFineDeposition::Transition This is a method that was developed in order to minimize spurious gradients in the density across the EB. This method operates by extracting the coarse-grid particles that live around the refinement zone (within some radius), and depositing them with the fine-grid particle width.

Important

Most coarse-fine particle deposition schemes exhibit some artifacts around the refinement boundary, especially when the particle width exceeds the grid cell size (e.g., for TSC). The CoarseFineDeposition::Transition method is the one that we recommend, especially when used with CIC, as it eliminates spurious gradients across the refinement boundary.

Particle interpolation

To interpolate mesh data onto a payload column, the user can call AmrMesh::interpolateParticles:

void
interpolateParticles(ParticleContainer<P, Traits>& a_particles,
                     const std::string&            a_realm,
                     const phase::which_phase&     a_phase,
                     const EBAMRCellData&          a_meshField,
                     const DepositionType          a_interpType,
                     const bool                    a_forceIrregNGP) const;

The template parameter Members is a pack of payload member pointers selecting the target column(s): a single scalar column, or SpaceDim columns for a vector quantity. For example, to interpolate a vector velocity stored as vx/vy/vz payload columns,

RefCountedPtr<AmrMesh> amr;

amr->interpolateParticles<D_DECL(&MyPayload::vx, &MyPayload::vy, &MyPayload::vz)>(...);

The companion AmrMesh::interpolateWeight interpolates a scalar mesh field onto the container-owned weight column.

Note

If interpolating onto a scalar column, the mesh variable must have exactly one component. Likewise, if interpolating a vector quantity (SpaceDim columns), the mesh variable must have SpaceDim components.

Particle visualization

Note

Particle visualization is currently a work in progress with limited functionality.

Simple particle visualization can be performed by writing H5Part compatible files which can be read by VisIt. This is done through the function writeH5Part in the DischargeIO namespace, with the following signature:

template <typename P, typename Traits>
void
writeH5Part(std::string                         a_filename,
            const ParticleContainer<P, Traits>& a_particles,
            RealVect                            a_shift,
            Real                                a_time) noexcept;

This routine permits particles to be written (in parallel, when using MPI) into a file readable by VisIt. The container-owned id, position, and weight are always written automatically. Which payload columns are exported is declared once, per particle type, through an optional h5PartColumns member of the ParticleTraits specialization. This descriptor co-locates each dataset name with the payload member pointer(s) it selects, so there is no positional coupling between names and accessors, and internal scratch columns are simply left out (opt-in). For example, ParticleTraits<ItoParticle> declares

/**
 * @brief Payload columns written to H5Part visualization files: the transport coefficients, the average
 * energy and the velocity vector (id, position and weight are written automatically; the scratch
 * and previous-position columns are internal and not plotted).
 */
static constexpr auto h5PartColumns = std::make_tuple(
  H5Part::Scalar{"mobility", &ItoParticle::mobility},
  H5Part::Scalar{"diffusion", &ItoParticle::diffusion},
  H5Part::Scalar{"energy", &ItoParticle::energy},
  H5Part::Vector{"v", D_DECL(&ItoParticle::vx, &ItoParticle::vy, &ItoParticle::vz)});

Each H5Part::Scalar descriptor becomes one dataset; each H5Part::Vector descriptor is written as name-x/name-y/name-z. If a particle type declares no h5PartColumns (or an empty one), only id, position, and weight are written. The argument a_shift will simply shift the particle positions in the output HDF5 file.

Merging and splitting particles

Often, merging or splitting of particles is required. chombo-discharge supports several strategies, which fall into two families:

• Per-cell mergers operate on one cell at a time and are implemented as factory functions in ParticleManagement that return a ParticleMerger<P, Traits> functor. Each factory accepts user-supplied lambdas for the particle-type-specific gather, reduce, and scatter steps, so the same algorithm can be reused with any ParticleSoA payload type. These are equal_weight_kd, reinitialize, reinitialize_bvh, and nn_sfc.

• Whole-hierarchy mergers operate collectively across the entire AMR hierarchy rather than one cell at a time. None is a per-cell factory, and all are distributed and MPI-safe. These are nn_pair_tree, nn_pair_onecell, nn_pair_hash, kd_carve, kd_patch, and kd_skin_nn.

The recommended pattern for the per-cell family is to cell-sort the leaf, extract a cell’s particles into a small scratch ParticleSoA, merge/split them, and rebuild the leaf. ParticleSoA<P>::extractCell performs the per-cell extraction

/**
 * @brief Extract cell c's particles (the CSR range [cellStart(c), cellStart(c+1))) into a_out.
 * @details a_out is cleared first, then this cell's particles are appended (all columns, including
 * id/rank). Requires the container to be cell-sorted. This is the building block of the per-cell
 * extract->process->rebuild idiom that replaces the AoS BinFab<List<P>> per-cell pattern: extract a
 * cell into a reusable scratch container, mutate it (append/remove), accumulate it into an output
 * container via append(const ParticleSoA&), and finally swap() the output into the leaf. See
 * ItoSolver::makeSuperparticles for the canonical use.
 * @param[in]  a_cell Cell index.
 * @param[out] a_out  Destination scratch container (cleared on entry; reuse it across cells).
 */
void
extractCell(const std::size_t a_cell, ParticleSoA& a_out) const
{
  CH_assert(m_sorted);

  a_out.clear();

  const std::pair<std::size_t, std::size_t> range = this->cellRange(a_cell);
  for (std::size_t i = range.first; i < range.second; i++) {
    a_out.appendParticle(*this, i);
  }
}

and the merged result is accumulated into an output ParticleSoA (via append) which finally replaces the leaf with swap. A complete worked example is ItoSolver::makeSuperparticles in $DISCHARGE_HOME/Source/ItoDiffusion/CD_ItoSolver.cpp.

Tip

The source code for all per-cell merger factories is in $DISCHARGE_HOME/Source/Particle/CD_ParticleManagement.H.

kd-trees

Several of the mergers below partition particles using kd-trees. kd-trees operate by partitioning a set of input primitives into spatially coherent subsets. At each level in the tree recursion one chooses an axis for partitioning one subset into two new subsets, and the recursion continues until the partitioning is complete. Fig. 14 shows an example where a set of initial particles are partitioned using such a tree.

Fig. 14 Example of a kd-tree partitioning of particles in a single cell.

The particles in each leaf of the tree can then be merged into new particles, one per leaf. What distinguishes the kd-based mergers from one another is the rule used to choose each split plane, the rule used to stop the recursion, and whether the tree is built per cell or per patch:

• equal_weight_kd and reinitialize_bvh build one tree per cell, splitting so that the two halves carry as nearly equal weight as possible.

• kd_carve and kd_patch build one tree per patch, splitting at the count or weight median and stopping on a live per-cell quota.

Per-cell kd-merges

The per-cell partitioner buildEqualWeightKDLeaves operates on a lightweight AoS particle type (ParticleManagement::MergeParticle) carrying the position and weight as data members, plus an opaque payload holding any quantities to be preserved across a merge. It recursively bisects the input particles into spatially coherent leaves whose weights are as equal as possible – at each bisection the two halves differ by at most one physical particle. It returns the leaf particle ranges directly; the tree is built in a flat, reusable scratch buffer rather than as linked node objects.

Note

ParticleManagement::detail::buildEqualWeightKDLeaves (in $DISCHARGE_HOME/Source/Particle/CD_ParticleManagement.H) lives in the internal detail namespace – it implements the partitioning but is not part of the public interface and is not meant to be called directly. The public entry point is the makeEqualWeightKDMerger factory (see below), which wraps it.

Warning

buildEqualWeightKDLeaves will usually split particles to ensure that the weight in the two subsets are the same (thus creating new particles). In this case any other members in the particle type are copied over into the new particles.

Equal-weight merging (equal_weight_kd)

ParticleManagement::makeEqualWeightKDMerger wraps buildEqualWeightKDLeaves into a ParticleMerger functor. The caller provides three lambdas:

• A gather function that packs one SoA slot into a MergeParticle.

• A reconcile function (BinaryParticleReconcile) that propagates payload fields to both daughter particles when the median particle is split across a kd boundary.

• A scatter-leaf function that receives the raw [first, last) pointer range of one leaf and appends exactly one merged particle to the SoA.

Since the weight in the leaves of the tree differ by at most one, the resulting computational particles also have weights that differ by at most one, as shown in Fig. 15.

Fig. 15 kd-tree partitioning of particles into new particles whose weight differ by at most one. Left: Original particles with weights between 1 and 100. Right: Merged particles.

In the weighted-centroid variant (equal_weight_kd in ItoSolver), the scatter-leaf computes the weight-averaged position and energy over all particles in the leaf. Particle weights need not be integers, but buildEqualWeightKDLeaves may create new particles at the kd boundaries (see warning above), so the total computational-particle count may exceed the target by a small amount during the build before being reduced.

Reinitialization algorithms

Reinitialization algorithms discard some or all of the spatial information within a cell and redraw particle positions rather than collapsing each group onto its centroid. This avoids accumulating merged particles at a cluster of centroid positions, at the cost of discarding fine-scale spatial information.

Reinitialization (reinitialize)

ParticleManagement::makeReinitializeMerger discards all spatial information and rebuilds the cell from scratch.

template <class Context, class P, class Traits = ParticleTraits<P>>
inline ParticleMerger<P, Traits>
makeReinitializeMerger(std::function<std::pair<long long, Context>(const ParticleSoA<P, Traits>&)> a_aggregate,
                       std::function<void(ParticleSoA<P, Traits>&, const RealVect&, long long, const Context&)> a_emit,
                       std::function<RealVect()> a_probLo) noexcept;

The returned functor proceeds as follows:

1. Calls aggregate once on the input SoA to obtain the total physical-particle count and a caller-defined context (e.g. the weight-averaged energy).

2. Passes the physical count to partitionParticleWeights, which divides it into at most ppc integer weights differing by at most one.

3. For each weight, draws a random position in the cell via Random::randomPosition (cut-cell aware) and calls emit to append the new particle.

All output particles share the same aggregated context, so per-particle information (e.g. individual energies) is lost. This method requires that particle weights are (close to) integers.

Tip

makeReinitializeMerger captures probLo at parse time. The cell-centre position is computed internally as probLo + dx * (gridIndex + 0.5), so no grid pointer needs to be retained in the returned functor.

kd-tree reinitialization (reinitialize_bvh)

The reinitialize_bvh variant uses the same makeEqualWeightKDMerger factory and the same per-cell kd partition as equal_weight_kd, but replaces the centroid scatter with a position-reinitialising scatter:

• Cut-cells (volFrac < 1): the weighted centroid is used to keep the merged particle inside the embedded boundary.

• Full cells: a random point is drawn uniformly from the bounding box of the leaf, reinitialising the spatial distribution within each kd partition rather than collapsing it to a single point.

Unlike reinitialize, the spatial information is discarded only within each leaf, so the coarse-grained distribution across the cell is preserved.

Note

In the full-cell branch, energy is accumulated over the leaf but is not normalised by weight – the stored value is the total (not average) energy of the leaf. This is intentional and matches the original ItoSolver behaviour.

Nearest-neighbour algorithms

Nearest-neighbour algorithms merge each particle with a spatially close partner rather than with everything in its cell. nn_sfc does this per cell, using a space-filling-curve ordering.

The three nn_pair_* variants are distributed, MPI-safe merges that operate collectively across the whole AMR hierarchy. All three share the same propose/judge/verdict protocol and differ only in how merge candidates are found. Each round proceeds as follows:

1. Ghost particles are refilled (fresh, exactly once per round) so that a particle’s nearest neighbour may be one owned by another patch or rank.

2. Every particle’s nearest neighbour is located, and pairs lying entirely within one patch are merged immediately (the trivial tier).

3. Pairs that straddle a patch or rank boundary are resolved through a single cross-patch propose/judge/verdict exchange, so both owners agree on exactly one merge and no particle is merged twice.

Because a round merges pairs, one round only roughly halves an over-full cell’s surplus above the target, so a call repeats the round until no round merges anything – or until nn_pair_max_rounds rounds have run, whichever comes first. A cell far above the target is therefore not necessarily drained by a single call, and in a running simulation the merge relies on being invoked once per time step – and on particle motion between invocations – for further convergence. Merged particles need globally unique ids that cannot collide across ranks or rounds, so the caller supplies an id allocator (a rank-namespaced counter suffices).

Important

All three nn_pair_* methods require a particle ghost mask (width-1 for nn_pair_onecell; width equal to the configured merge distance, or 1, for nn_pair_tree/nn_pair_hash), which is only built during a regrid. The mask must therefore be registered before the grids are (re)built – registering it late leaves it empty and the neighbour search will not see cross-patch particles. ItoSolver handles this automatically when merge_algorithm selects any of the three, and ItoKMCStepper does the same when its regrid-time merge is set to any of the three.

In ItoSolver the behaviour is tuned through nn_pair_iterate (repeat the local trivial-tier merges within a round until no further local pairs remain), nn_pair_fallback (how many additional candidate neighbours to consider when the nearest is unavailable), nn_pair_max_rounds (cap the number of rounds a single call may run, so the cost stays bounded rather than running to full convergence), and, for nn_pair_tree/nn_pair_hash only, nn_pair_max_cell_dist (cap the neighbour search radius in cells; nn_pair_onecell’s search radius is fixed at 1 and does not read this).

SFC pair merging (nn_sfc)

ParticleManagement::makeSfcNearestNeighborMerger sorts particles along a Hilbert space-filling curve and merges adjacent pairs until the count is at most ppc.

template <class Packed,
          Real Packed::*packWeight,
          RealVect Packed::*packPosition,
          class P,
          class Traits = ParticleTraits<P>>
inline ParticleMerger<P, Traits>
makeSfcNearestNeighborMerger(std::function<Packed(const ParticleSoA<P, Traits>&, std::size_t)> a_gather,
                             std::function<void(Packed&, const Packed&)>                       a_combine,
                             std::function<void(ParticleSoA<P, Traits>&, const Packed&)>       a_scatter) noexcept;

The caller provides three lambdas:

• A gather function that packs one SoA slot into a MergeParticle.

• A combine function that merges two adjacent intermediates in place (typically a weighted average of position and energy).

• A scatter function that unpacks one merged intermediate back into the SoA.

Unlike the kd-tree methods, SFC merging does not require integer weights. Particle counts below ppc are handled by splitting the heaviest particle: its weight is halved and a copy is appended, repeating until the target is reached (only if the heaviest particle has weight \(\geq 2\)). The Hilbert ordering ensures that merged pairs are spatially close, which better preserves spatial correlations than random pairing and typically produces smoother merged distributions than the kd centroid.

Whole-patch BVH search (nn_pair_tree)

nn_pair_tree is implemented as ParticleManagement::mergeNearestNeighborsTree in $DISCHARGE_HOME/Source/Particle/CD_NearestNeighborParticleMerge.H, and searches for candidates via one whole-patch PointCloudBVH per patch.

Per-cell BVH search (nn_pair_onecell)

nn_pair_onecell is implemented as ParticleManagement::mergeNearestNeighborsOneCell in the same file, and instead builds one PointCloudBVH per occupied grid cell. A query only ever searches its own cell and its Moore-adjacent neighbours, so its merge distance is structurally fixed at Chebyshev cell distance 1.

Whole-patch hash-grid search (nn_pair_hash)

nn_pair_hash is implemented as ParticleManagement::mergeNearestNeighborsHash in the same file – identical algorithm and tunables to nn_pair_tree, but searches via one whole-patch PointCloudHashGrid (a uniform spatial hash grid) per patch instead of a PointCloudBVH. Both point-cloud types expose an identical query interface, and the two entry points share one generic implementation selected via a template parameter over the point-cloud type.

Whole-patch kd-merges

kd_carve and kd_patch are not built around pairwise nearest-neighbour matching at all. Both are distributed, MPI-safe whole-container merges declared in $DISCHARGE_HOME/Source/Particle/CD_KDParticleMerge.H, and both use the same tree build and the same per-cell quota; they differ only in how patch boundaries are treated.

Each patch’s particles are split purely by position into a kd-tree, never snapped to the grid, since particles near a cell face must be able to merge across it. The longest axis is bisected at the count median – the plane putting half the node’s particles on either side – while the node is still larger than kd_split_weight_leaf_dx cell widths, and at the weight median – half the node’s particle weight on either side – once it is smaller. Above that scale a split’s job is to apportion leaves between cells, which is a question about counts; below it there is little left to apportion and the weight median drives the resulting super-particles toward equal weight.

Splitting is governed by a live per-cell quota. Every group becomes exactly one super-particle, placed at its weighted centroid, so the number of groups centred in a cell is that cell’s post-merge population; that count is tracked while splitting, and a split that would push a cell past the target is refused. Groups are split heaviest-first, because the quota is claimed first-come-first-served: ordering by particle count instead lets light groups take a cell’s slots before a heavy group is considered, stranding the heavy one whole as a single very heavy super-particle.

The quota constrains only those groups that can actually merge. A group whose largest per-axis extent exceeds one cell width is left untouched, so its members survive individually; it therefore costs its cell its full member count rather than one particle, and the quota does not block splitting it further. A cell can consequently finish slightly above target when a group is both too wide to merge and in a cell already at quota, since either choice overshoots. This extent bound is applied only as a filter on the finished partition; it never drives the splitting itself.

Particles in cells at or below target are excluded from the partition and emitted as singletons. A singleton never merges (the merge threshold is two members), so an under-populated cell cannot be drained by a neighbour’s merge.

Unlike the nn_pair_* family there is no iteration, fallback or max-cell-distance equivalent, since there is no drain loop and the per-cell quota alone governs how far splitting goes. The maximum extent of a mergeable group is fixed at one cell width rather than exposed as an option: it is a physical bound rather than a dial, and it matches the hardcoded width-1 ghost halo. One option applies to both algorithms:

• kd_split_weight_leaf_dx – the group size, in cell widths, at or below which the split plane switches from the count median to the weight median. 0 disables the weight median entirely.

Carved patch boundaries (kd_carve)

kd_carve is implemented as ParticleManagement::mergeKDCarve. It builds the tree over each patch’s local plus ghost particles, and then classifies every group in one of three ways:

• Interior – none of its members’ cells are exposed to any neighbouring patch or level, and it holds no ghost member. Such a group merges immediately, with zero communication.

• Unmergeable – the group’s largest per-axis extent exceeds one cell width, so merging it would combine particles more than a cell apart. Left untouched.

• Boundary – otherwise. Such a group becomes a candidate “box” for the carve step.

The boundary/carve step (“z-buffer carve”) resolves every contested particle in a single, fixed, non-iterative protocol. Each box is keyed by (AABB volume, anchor particle id), compared so that the tightest box wins a contested particle (the anchor id is a deterministic tiebreak, since volume alone is not a strict total order). Because a box’s membership can be whittled down by losing individual particles to other, unrelated competing boxes, a box that nominally wins one particle is not guaranteed to end up with enough survivors to actually commit a merge – so the protocol runs in three rounds rather than two:

1. Claims are routed to each particle’s owner.

2. Each owner replies with the nominal winner.

3. Once a box’s final, post-contest membership and validity are known, it tells each surviving foreign member’s owner to actually delete it.

Only this last step authorizes a deletion; a box’s participants are never removed on the strength of the nominal result in step 2 alone, since that could otherwise delete a particle whose “winning” box turns out not to commit.

A patch only ever sees its own particles plus a width-1 ghost halo, so a group can only ever list members drawn from that halo, and claims therefore only ever travel to a patch’s direct neighbours – the exchange is bounded and fixed-shape, not a drain loop, regardless of the third round. This also matches the leaf-extent bound exactly: a mergeable group is at most one cell across, and the halo supplies precisely that reach.

Important

kd_carve requires a width-1 particle ghost mask, registered unconditionally by ItoSolver (the same requirement and registration timing as the nn_pair_* family – see above).

Patch-local merging (kd_patch)

kd_patch is implemented as ParticleManagement::mergeKDPatch. It selects the same tree build and per-cell quota without the boundary tier: no ghost halo is filled, and every group within the extent bound merges immediately regardless of boundary exposure. No particle is ever contested, so there is no claim/verdict/commit exchange and no communication at all.

Each patch therefore reduces only the particles it owns. This is cheaper than kd_carve and free of communication, at the cost of no coordination across a patch boundary: a cell straddled by two patches is reduced by each of them independently.

Warning

The caller must not fill a ghost halo before calling kd_patch. A ghost particle reaching the merge would be merged locally while its true owner merges it too, double-counting its weight.

Nearest-neighbour skin (kd_skin_nn)

kd_skin_nn keeps the tree build and per-cell quota but replaces the carve with the nearest-neighbour pair merge, splitting the particles across two containers so that neither tier can merge what the other already did.

First, ParticleManagement::mergeKDInterior commits every group that holds no ghost into a separate container, and removes its members from the working one. Holding no ghost is the whole safety condition, and it costs no communication to check: every member of such a group is resident in this patch, and any other patch that draws one of them into a group of its own necessarily sees it as a ghost, which disqualifies that group there by the same test. No two patches can therefore commit the same particle. Note this is a weaker condition than kd_carve’s – boundary exposure is not consulted at all, so an exposed but uncontested group merges here rather than going to arbitration, which is what makes the skin smaller.

What remains in the working container is the skin: groups holding a ghost, plus groups too wide to collapse. It is drained by ParticleManagement::mergeNearestNeighborsOneCell (see above), whose search radius is fixed at Chebyshev cell distance 1 and so matches the width-1 halo and the one-cell extent bound exactly. The skin tier drains each cell not to the full target but to what is left of it: every interior super-particle reserves a slot in the cell it landed in, so the two tiers share one budget. Reservations are counted over ghosts as well as local particles, because the nearest-neighbour tier consults the occupancy of cells owned by a neighbouring patch when it judges a cross-patch pair.

Finally the interior container is folded back in, and under-full cells are brought up to the target in the usual way.

Important

kd_skin_nn requires a width-1 particle ghost mask, exactly as kd_carve does. It is tuned by nn_pair_iterate, nn_pair_fallback and nn_pair_max_rounds, shared with the nn_pair_* family; nn_pair_max_cell_dist does not apply, since the skin tier’s radius is structurally fixed at 1.

Because a skin particle and an interior super-particle live in different containers, they can never merge with each other. A cell holding interior results plus a single leftover skin particle therefore finishes one particle above the target – the same class of residual as kd_carve’s, and not a conservation error.