Add transient subcell (TAS) collision-partner selection method (CPU + KOKKOS)

[stanmoore1]

Purpose

Adds Bird's transient adaptive subcell (TAS) collision-partner selection method as an alternative to the existing nearcp near-neighbor option, enabled with a new collide_modify subcell yes keyword. Each timestep the particles in a collision cell are binned into a transient grid of subcells (~1 particle/subcell), and a collision partner is drawn from the same subcell, or from expanding shells of neighbor subcells if the home subcell is empty. Like nearcp, this keeps the mean collision separation (MCS) small even when cells are much larger than the local mean free path, which allows accurate simulations on coarser grids with fewer particles.

The method was prototyped on the long-lived subcell branch. This PR brings that branch up to date with master, rewrites the CPU implementation on top of the current Collide structure (fixing a number of bugs found in the prototype), documents and exercises it with examples, and adds a validated KOKKOS port so it runs on CPU threads and GPUs.

Author(s)

Stan Moore, Sandia National Laboratories, with an original CPU prototype from the subcell branch. Implementation completed, corrected, KOKKOS-ported, and verified with assistance from Claude Code.

Backward Compatibility

Fully backward compatible. The feature is opt-in via the new collide_modify subcell keyword, which defaults to no. With subcell no, all code paths (standard random, nearcp, ambipolar, both CPU and KOKKOS) are byte-for-byte unchanged — verified against the existing blessed example logs and against master behavior. No input syntax changes for existing decks.

subcell is mutually exclusive with nearcp, and (for now) errors out cleanly if combined with the ambipolar approximation or with multiple collision groups.

Implementation Notes

Algorithm (src/collide.cpp, Collide::collisions_one_subcell). Templated on dimension. Per cell: Nsub = floor(N^(1/dim)) subcells per direction so the subcell count does not exceed the particle count; particles are binned into a head-insertion linked list (subcell_first/subcell_next/subcell_count, flattened subcell_id); the first particle is random over the cell, the partner is random within its subcell, else random over an expanding, grid-clipped shell of subcells (subcell_ring). A previous-collision-partner exclusion (reusing the existing nn_last_partner machinery) prevents an unphysical immediate re-collision of the same pair, choosing a different/next-nearest partner instead, per Bird. The transient binning is discarded after each cell, so memory overhead is negligible, and it is rebuilt after any reaction that adds/removes a particle so the subcell lists stay consistent with plist.

KOKKOS port (src/KOKKOS/collide_vss_kokkos.cpp). Mirrors the existing nearcp Kokkos path: the parallel_for/parallel_reduce runs one grid cell per thread, and the subcell scratch is kept thread-private by indexing it as (icell, *) view rows — exactly like d_nn_last_partner. Five 2D device views hold the per-cell subcell state, sized (nglocal, maxcellcount) and reallocated alongside d_plist including in the react-retry path. Partner selection (find_nn_subcell) and rebinning (rebin_subcell) are DIM-templated KOKKOS_INLINE_FUNCTIONs that reproduce the CPU shell search and exclusion exactly. The existing react retry/backup machinery is reused unchanged; the subcell scratch needs no backup since it is rebuilt per cell and after each reaction.

Bugs fixed relative to the prototype: infinite loop on an unknown collide_modify keyword (missing final else error); a destructor memory leak (recomb_ijflag); a heap overflow in the particle-creation growth path (grew the wrong arrays); stale subcell lists after reactions could select deleted or uninitialized particles; pow(N,1/3) integer-truncation of the subcell count (now cbrt/sqrt with a tolerance); an unused npmax*npmax allocation (hundreds of MB for large cells) and an unused vmax array; and out-of-range subcell indices for particles on a cell's upper face (now clamped). Development scratch files and stray debug code from the prototype were removed.

Verification.

• Regression: with subcell no, examples/collide/in.collideInterspecies reproduces the repository's blessed logs column-for-column; the standard and nearcp KOKKOS paths still match the CPU.
• KOKKOS correctness: built with -DSPARTA_KOKKOS_EXACT=ON, the subcell method reproduces the non-KOKKOS results bit-for-bit (Natt, Ncoll, temperature identical every stats step) for the 2d and 3d subcell examples and for a dissociating-gas case that exercises particle creation and the retry path (Np grows 10000→10751, Nreact and T all match). Multithread and production-RNG runs are physically consistent.
• Physics: equilibrium Ar box (273 K), collision rate matches VHS theory to within statistics — standard 0.9999, nearcp 0.9991, subcell 0.9984 of the analytic rate (the small deficits are the expected effect of suppressing repeat collisions); equilibrium temperature exactly conserved.
• Memory safety: dissociation (particle creation), high-rate recombination (deletion), 3000-particle cells (array growth), and 2d all run clean under AddressSanitizer.
• Mean collision separation: for random selection MCS = 0.662·Δx exactly as theory predicts and scales with cell size; subcell MCS is grid-independent (2.355e-6 m on both 1e-5 and 2e-5 m cells) where nearcp's MCS grows with the cell.

Performance vs near-neighbor.

Selection cost — collide-section seconds for 300 steps, 80k-particle 3d box, serial:

grid	particles/cell	std	nearcp 10	nearcp 30	subcell
20³	10	1.05	1.45	1.35	1.53
10³	80	1.07	1.52	1.78	1.69
5³	640	1.00	1.37	1.71	1.63
3³	2963	1.02	1.47	1.69	1.59

Subcell costs about the same as nearcp 30 and, unlike nearcp, is O(1) per selection independent of cell population (no Nlimit to tune) while its partner quality does not degrade as cells coarsen.

Coarse-grid accuracy — Fourier heat conduction between 223/323 K diffuse walls (Kn≈0.024), wall heat-flux error vs a fine-grid (cell≈λ/3) reference (±2–3% statistical noise):

cell size	std	nearcp 30	subcell
λ/3 (ref)	−0.3%	+0.5%	−0.2%
λ	+0.3%	−4.0%	−2.5%
4.2 λ	+26.7%	+1.0%	+2.2%

On cells ~13× coarser than the λ/3 rule, standard selection is off by 27% while subcell stays within noise of the reference. Because it holds accuracy on coarse grids, subcell lets the grid and particle count shrink together: a coarse-grid subcell run with 12.8× fewer particles matched the fine reference to +2.6%, ran ~15× faster per step, and used half the particle memory — consistent with the ~2× overall efficiency reported by Gallis & Torczynski for the sophisticated-DSMC scheme.

Post Submission Checklist

• The feature or features in this pull request is complete
• Suitable new documentation files and/or updates to the existing docs are included (doc/collide_modify.txt)
• One or more example input decks are included (examples/collide/in.collide_{2D,3D}{,_nn,_sub} with reference logs for 1 and 4 MPI ranks)
• The source code follows the SPARTA formatting guidelines

Further Information, Files, and Links

• The KOKKOS subcell path currently supports single-group, non-ambipolar collisions (same scope as its nearcp counterpart); ambipolar and multi-group support is left as follow-up.
• src/MAKE/Makefile.asan is included as a convenience AddressSanitizer build config used during verification (analogous to the existing Makefile.serial_debug); reviewers who prefer not to ship it can drop that one file with no effect on the feature.
• References for the method: G. A. Bird, M. A. Gallis, J. R. Torczynski, D. J. Rader, "Accuracy and efficiency of the sophisticated direct simulation Monte Carlo algorithm for simulating noncontinuum gas flows," Phys. Fluids 21, 017103 (2009), DOI:10.1063/1.3067865; M. A. Gallis, J. R. Torczynski, "Effect of Collision-Partner Selection Schemes on the Accuracy and Efficiency of the DSMC Method," Int. J. Numer. Meth. Fluids 67(8):1057-1072 (2011), DOI:10.1002/fld.2409.

🤖 Generated with Claude Code

---

Generated by Claude Code