[max-models]

Stack: part 2 of 2, based on #774 (merge that first). This PR's diff includes the commits of #774 until it is merged. #774 changed the marker buffer to one column per marker ((n_cols, n_markers_max), markers[j, ip]), and this branch is merged with it, so "marker slot" below means a column of the buffer. Following #774, the names count markers, not rows: the bound is n_markers_used (was n_rows_used), counted by count_markers_used and tracked by _track_markers. Both branches are also merged with the current devel, including the kernel-based box sorting of #736.

Solves the following issue(s):

Closes #707 together with #774. This is part 2 of the issue: launches and per-quantity operations cover only the marker slots that can hold markers, not the trailing holes of the buffer.

Core changes:

The marker buffer has room for (1 + bufsize) * Np markers (n_markers_max slots), and every kernel and per-quantity operation ran over all of them. Particles.n_markers_used is now a bound with the invariant that every marker slot from it on is a hole, and the per-marker work stops there:

• Kernels:
  • args_markers.n_markers holds the bound, so the pyccel kernels loop over it and the CUDA kernels return after it.
  • The 15 accumulation kernels that looped over the full marker count from the array shape now use args_markers.n_markers.
  • Pusher and the accumulators launch n_markers threads.
• Python operations on markers[:, :n_markers_used]:
  • the pusher preamble (prepare_push gets the sliced Array2D view on the GPU);
  • the residual check of the iterative pushers;
  • update_holes;
  • the boundary conditions (_find_outside_particles and the snapshot of the positions);
  • the alpha-weighted sorting positions of mpi_sort_markers (the first n_markers_used columns of the (3, n_markers_max) helper arrays).

Keeping the bound without reading the device:

• New markers always fill the lowest holes (nonzero of the hole mask, ascending). So k incoming markers end below (occupied slots) + k, where _n_occupied_bound is an upper bound on the occupied slots.
• The sent, received and removed counts are already host integers: send_inds.size, the NumPy recv_info, and the lost-marker counts of the boundary conditions. So _track_markers updates the bound on the host. On the CuPy backend the time loop gets no new device-to-host transfer, which the transfer budget of the GPU end-to-end tests checks.
• Where the bound is updated:
  • counted exactly with count_markers_used (reads row 0; one device read on CuPy): after loading (fresh, external or restart; restart files can have markers in any slot) and after box sorting (do_sort, also with the box kernels of Set up sorting boxes on the CuPy backend #736), which moves holes anywhere. Code that writes markers into holes itself, outside Particles, must call it too; it is public for that reason.
  • tracked from host counts: the MPI exchange, removal at the boundary, and the SPH ghost-box exchange (own copy and receives);
  • unchanged on purpose: holes created without a count (removed ghost particles) leave the bound too high, which is safe.

Fix included (also on #774): test_verif_LinearVlasovAmpereOneSpecies::test_delta_f_initial_condition checked the mean of the delta-f weights per rank, assuming every rank holds a random subset of the markers.

• The markers are sorted into the domain of their rank, so on four ranks each rank holds a quarter of the cos(kx) perturbation: mean/rms = (2/π)/(1/√2) ≈ 0.9, exactly the measured 1.25e-7/1.39e-7.
• The sums are now reduced over the clone's communicator.
• The test fails on four ranks without any of these changes, but testmon had not selected it on recent PRs (0 of 3184 verification tests selected on 726-sorting-boxes-cupy). Touching pic/base.py in Store the marker buffer column-major (struct-of-arrays), part 1 of #707 #774 selected it in the domain cloning in Container job.

Tests:

• test_markers_used (in test_marker_layout.py), run on NumPy and on the fake CuPy, on 1 and 2 ranks. It uses removing boundaries in eta1, pushes that move markers out, and MPI sorting. It checks:
  • after every step, that every marker slot after n_markers_used is a hole and that args_markers.n_markers matches;
  • that the final markers (by ID) and the number of lost markers are bit-for-bit those of the same run with every slot counted as used.
• bench_marker_layout.py (from Store the marker buffer column-major (struct-of-arrays), part 1 of #707 #774): before the merge with the new layout, the valid markers after all timed operations were bit-for-bit those of devel. With 1M markers there are 1,251,008 slots, of which n_markers_used = 1,000,000.
• Run locally: the new test and the verification test above, and the marker-layout tests after the renaming. The merges with Store the marker buffer column-major (struct-of-arrays), part 1 of #707 #774 and devel were resolved without running the suite locally; CI ran on the merged state, where test_markers_used called the wrong device-context helper after the merge (fixed).

Benchmark:

• The work of every restricted operation scales with n_markers_used / n_markers_max. That is 0.8 for the default bufsize = 0.25 right after loading, and about 0.5 for bufsize = 1 (the 8M-slot buffer for 4M markers in the issue's trace).
• I could not measure it here. bench_marker_layout.py --all-markers-used compares the same build with and without the bound. I ran it twice, alternating, but the machine had a load average of about 500 from other jobs, and operations this PR does not touch varied by up to 15× between runs. So I don't claim a number.
• It needs a quiet node or the GPU (First GPU run of the CUDA tests (H100) #687):
  • bench_marker_layout.py --all-markers-used --json all.json
  • bench_marker_layout.py --json used.json
  • bench_marker_layout.py --compare all.json used.json

Model-specific changes:

None.

Documentation changes:

Particles.n_markers_used and the tracking (_track_markers, count_markers_used) are documented in their docstrings.

🤖 Generated with Claude Code

[max-models]

Benchmark of this PR, on top of #774 (setup as in my comment there: M1 MacBook Air, NumPy, 1 rank, one thread, bench_marker_layout.py, 1M markers, 3 alternating rounds, median; another session was using three cores throughout, so about 10–20% noise).

The rows used line of the script confirms the bound: 1,000,000 of 1,251,008 rows (0.8, the default bufsize = 0.25). I compared the same build with and without the bound (--all-rows), and against #774:

Operation	#774	#775	Speedup	all rows → bound (same build)
mpi_sort_markers	77.5 ms	63.8 ms	1.22×	1.19×
update_holes	1.0 ms	0.8 ms	1.31×	1.18×
apply_kinetic_bc	10.1 ms	9.0 ms	1.11×	1.10×
preamble	4.3 ms	4.5 ms	0.97×	0.96×
push_eta_stage	282 ms	280 ms	1.01×	1.00×
push_vxb_analytic	733 ms	707 ms	1.04×	0.99×
charge_density_0form	481 ms	481 ms	1.00×	0.99×
binning	312 ms	307 ms	1.01×	0.98×

• The gain is the expected one for the three column operations, bounded by 1 − 0.8 = 20% at this bufsize.
• The pushers and the accumulation do not speed up. Their pyccel loops already skip a hole with one branch, so the 20% fewer rows cost almost nothing to visit.
• The preamble differs by 0.03×, which is inside the noise (run-to-run spread 19%).
• With bufsize = 1 the bound would be about 0.5. I did not measure that.

Against the old row-major layout (devel) the two PRs together give 16× on the preamble, 9.5× on update_holes, 3.5× on apply_kinetic_bc, 2.6× on mpi_sort_markers, and 1.1–1.3× on the pushers. Total of all ten operations: 2312 ms → 1909 ms (1.2×).