Store the marker buffer column-major (struct-of-arrays), part 1 of #707

[max-models]

Solves the following issue(s):

Part 1 of #707: the marker buffer is stored as struct-of-arrays (one contiguous block per quantity) on both backends. Part 2 of the issue (tracking the range of marker slots that can hold valid markers, so that launches skip trailing holes) is #775, stacked on this PR.

Core changes:

Particles.markers is now a plain C-ordered array of shape (n_cols, n_markers_max): one column per marker, indexed markers[j, ip] (was (n_rows, n_cols), markers[ip, j]). Every quantity markers[j] is one contiguous block of n_markers_max values.

Naming: n_cols keeps its meaning, the number of quantities per marker. The number of marker slots (holes included) no longer counts rows, so it is renamed after what it counts: Particles.n_rows → n_markers_max, padded_n_rows → padded_n_markers, ROW_ALIGNMENT → MARKER_ALIGNMENT. The quantity indices are still called columns (first_pusher_idx etc.).

An earlier version of this PR kept the shape and allocated the buffer with order="F". That gives the same bytes in memory, but every array passed to a kernel then had to carry the order (pyccel float[:, :](order=F), conversion helpers at every entry point, checks in the CUDA arguments). The buffer is now transposed and the indexing changed instead, so no array needs an order flag.

• struphy.pic.marker_layout:

  • marker_array(n_cols, n_markers_max), plus padded_n_markers/MARKER_ALIGNMENT: n_markers_max is padded to a multiple of 32, so every row of the buffer starts on a 256-byte boundary (point 3 of the issue's second comment).
  • The module docstring documents the convention and its one pitfall. Selecting markers (markers[:, valid_mks]) gives a Fortran-ordered copy in NumPy, so convert it with xp.ascontiguousarray before passing it to a kernel or to MPI.

• pyccel kernels:

  • All kernels index markers[j, ip] and declare a plain float[:,:] (C order).
  • The marker count is shape(markers)[1], and MarkerArguments.n_markers is markers.shape[1].
  • The sorting kernels copy whole markers as columns.

• CUDA:

  • The kernels index markers(j, ip). CudaMarkerArguments takes the C-contiguous (n_cols, n) array again.
  • cunumpy infers the launch size from the first axis of the first array, which is now n_cols. So Pusher and the accumulators pass n_threads=args_markers.n_markers, and the geometry and spline marker kernels pass markers.shape[1].
  • prepare_push (the fused preamble from Fuse the pusher preamble into one CUDA kernel #708) takes an Array2D<double> view (a raw double* would be made C-contiguous, i.e. copied, by the launch). It numbers its threads row by row, so neighbouring threads touch neighbouring markers.

• Particles:

  • The properties follow the layout: markers_wo_holes(_and_ghost) are (n_cols, n), positions is (3, n), velocities is (vdim, n), and phasespace_coords, f_coords and f_jacobian_coords are (k, n). The setters check these shapes.
  • Holes are markers[0] == -1, ghosts are markers[-1] == -2, and whole markers are selected with markers[:, sel].
  • The sorting etas, the bc buffer and _is_on_domain_of_rank work on (3, n).

• MPI exchange:

  • The send lists, the receive buffers and the SPH ghost-box arrays are (n_cols, n).
  • The send lists and ghost-box markers are made C-contiguous. A column selection is Fortran-ordered, and MPI would otherwise send it marker by marker into C-ordered receive buffers.
  • The CuPy exchange receives into one buffer per source rank, because a slice of a common (n_cols, n_recv) buffer is not contiguous.

• Evaluation at markers takes (k, n) arrays, one column per point (eta1..3 in rows 0..2), so the buffer and its selections go in without a transpose:

  • Domain (metric coefficients, pull/push/transform, show), SplineFunction, the GVEC equilibria, the kinetic-background moments and B0 callables, and the canonical Maxwellian (coords[0] is (5, n));
  • user callbacks that are evaluated at markers (the example params files read etas[0][i]).

  Output shapes are unchanged. Domain and SplineFunction make the input C-contiguous (no copy for the buffer).

• Gather only the quantities that are used (point 4 of the issue's second comment): binning, the box count of the sorting boxes, the magnetic energy of DriftKineticElectrostaticAdiabatic, and two emptiness checks in the weight initialization.

• Box sorting from devel (Set up sorting boxes on the CuPy backend #736): devel is merged into this branch. Its box kernels (box_particles, fill_boxes, sort_and_fill_boxes, box_and_sort_particles) call the sorting kernels, and the two it adds there (count_particles_in_boxes, find_holes) now read the box row and row 0 of each marker column, like the other sorting kernels.

• HDF5 output and restart use the buffer layout:

  • restart data restart/<species> is (nt, n_cols, n_markers_max);
  • saved markers kinetic/<species>/markers are (nt, n_cols, n_to_save);
  • external marker files have a markers dataset of shape (n_cols, N).

  Restart and external files written before this change do not load. _load_restart copies the common markers when the padding differs and asserts that any extra ones in the file are holes. The post-processed orbits keep their (n_markers, n_quantities) layout.

• Fix found along the way: the info output of VlasovAmpereCoupling read particles.bufferindex, which does not exist, so it always failed. It now reads the velocities saved by the preamble, gathering 6 rows instead of whole markers.

Tests:

• pic/tests/test_marker_layout.py, on NumPy and on the fake CuPy (emulated CUDA launches), on 1 and 2 ranks. It checks:
  • the layout helpers;
  • that the buffer is C-ordered (n_cols, n_markers_max) with n_markers_max a multiple of 32, and stays the same object through loading, RK4 and v×B pushes, MPI sorting, boundary conditions, evaluation at markers and binning;
  • that marker selections are Fortran-ordered and accepted by MarkerArguments and a kernel after ascontiguousarray;
  • that evaluation on selections gives the same values as on the buffer.
• All tests that build markers or evaluation points by hand use the new layout. The CUDA parity cases pass n_threads explicitly.
• The box checks added by Set up sorting boxes on the CuPy backend #736 in test_sorting.py and test_sorting_device.py read the box row of the markers.
• test_mhd_equils checked its eta1 case with kind="markers" left over from the previous case; fixed.
• Run locally before the full suite was left to CI: the marker-layout, sorting, draw-parallel and load-external tests on 2 ranks, the emulated CUDA prepare_push/parity cases, and the geometry, feec, equilibrium, model-scalar, post-processing and output tests, plus test_verif_LinearVlasovAmpereOneSpecies on 1 and 2 ranks. The MPI tests caught the Fortran-ordered send lists described above. The SPH and pusher tests were stopped partway (no failures until then), and pyccel compilation and the GPU tests run only in CI. Nothing was run locally after the merge of devel (including the new box kernels); CI covers it.

Benchmark (CPU):

pic/tests/bench_marker_layout.py times the layout-dependent operations on the active backend. The numbers below were measured with the order="F" version of this PR. The buffer has the same bytes in memory now, so they are expected to carry over, but they have not been measured again. The setup is 1M markers (1.25M marker slots × 26 quantities), a Cuboid, 8³ cells and cubic splines, with the median of 5 calls on NumPy, 1 rank, no OpenMP, on an M-series MacBook Air. Before: devel at 11adb7b.

Operation	Before (row per marker)	After, run 1	After, run 2	Speedup
preamble (row assignments of Pusher._push)	527 ms	14 ms	34 ms	15–37×
update_holes	53 ms	1.1 ms	0.9 ms	47–61×
apply_kinetic_bc (periodic)	264 ms	70 ms	71 ms	3.7×
mpi_sort_markers (1 rank)	1277 ms	507 ms	519 ms	2.5×
weights of valid markers (particles.weights)	134 ms	40 ms	43 ms	3.1–3.3×
binning e1_v1	2728 ms	2057 ms	2311 ms	1.2–1.3×
push_eta_stage (RK4 pusher)	2685 ms	2021 ms	1985 ms	1.3×
push_vxb_analytic (pusher)	6213 ms	4977 ms	4861 ms	1.25–1.3×
charge_density_0form (accumulation)	4073 ms	3087 ms	3231 ms	1.25–1.3×
all quantities of the valid markers (markers_wo_holes_and_ghost)	225 ms	385 ms	392 ms	0.6×

Caveats:

• The machine was heavily loaded during all runs (load average 210–265 from other jobs), so the absolute numbers are inflated and noisy.
  • The operations on single quantities gain an order of magnitude or more, far beyond the noise.
  • The 1.25–1.35× on the kernel-bound rows (pushers, accumulation) should not be read as a kernel speedup until it is measured on a quiet node. The issue's quiet single-thread microbenchmark found kernel time unchanged within ±3%.
• Gathering all quantities of the valid markers is about 1.7× slower, as expected: a marker now touches 26 separate rows. The hot paths above no longer do it. It is still used by markers_wo_holes_and_ghost/markers_wo_holes and by the saved-markers output.
• With the order="F" version, the valid markers after all timed calls were bit-for-bit equal to devel (--save-markers, which still writes one row per marker).
• Not measured yet:
  • The GPU, where most of the gain is expected (preamble, update_holes, valid_mks masks, diagnostics). It waits for the first H100 run (First GPU run of the CUDA tests (H100) #687): CUNUMPY_BACKEND=cupy python src/struphy/pic/tests/bench_marker_layout.py.
  • CPU with OpenMP on a compute node.
  • SPH neighbour loops, whose per-marker access gets more expensive with this layout (see the issue).

Model-specific changes:

None beyond indexing: every model uses the new layout through Particles. User callbacks evaluated at markers (density, temperature or perturbation functions in params files that accept a single 2D array) now receive (3, n); the example params files are updated.

Documentation changes:

• "Marker array layout" section in the developer docs (dev-pic.rst, automodule of struphy.pic.marker_layout).
• The marker-layout question in CUDA_STRATEGY.md is marked as decided, with the new convention. The kernel and argument docstrings describe the layout.

🤖 Generated with Claude Code

[max-models]

Note this very important caveat: "The machine was heavily loaded during all runs (load average 210–265 from other jobs), so the absolute numbers are inflated and noisy.". So the results in the table should not really be trusted.

[max-models]

Benchmark of this PR against the old row-major layout, on a quiet-ish machine (the table in the description was measured at load average 200+).

Setup: M1 MacBook Air (8 cores), NumPy, 1 rank, OMP_NUM_THREADS=1, bench_marker_layout.py with 1M markers (1,251,000 rows × 26 columns), cubic splines, 8³ cells. Three builds in separate worktrees: devel at 11adb7b (before), this branch at db4a1d9, and #775. 3 rounds, alternating the builds, median of the 3 runs (each run is the median of 5 calls). Three Python processes from another session kept running at 100% CPU throughout, so expect about 10–20% noise; the round 1 (load 14) and round 3 (load 5) speedups agree to within a few percent.

Operation	Before (C)	This PR (F)	Speedup
preamble (column assignments of Pusher._push)	72.6 ms	4.3 ms	16.7×
update_holes	7.2 ms	1.0 ms	7.2×
apply_kinetic_bc (periodic)	31.3 ms	10.1 ms	3.1×
mpi_sort_markers (1 rank)	165 ms	77.5 ms	2.1×
weights of valid markers	16.8 ms	7.4 ms	2.3×
push_eta_stage (RK4)	358 ms	282 ms	1.27×
push_vxb_analytic	806 ms	733 ms	1.1×
binning e1_v1	348 ms	312 ms	1.1×
charge_density_0form	478 ms	481 ms	1.0×
full rows of valid markers (markers_wo_holes_and_ghost)	29 ms	49 ms	0.6×

Summed over all ten operations: 2312 ms → 1957 ms (1.2×).

• The big gains are on the operations that touch single columns. The pushers gain 10–30%, because their kernels spend most of the time in spline evaluation, not in memory traffic on the markers. The accumulation does not change.
• The one regression is gathering whole rows (markers[valid_mks]), as expected for column-major storage. Row gathers that remain are in diagnostics, not in the time loop.

[spossann]

Looks really nice this way! Please update the description before merging.