Hand-written CUDA vs code generation: port bstar_parallel_3form both ways (decision for #690)

[max-models]

Stack: part 8 of 14, based on #714 (merge that first), next: #720. Full order: #709 → #708 → #711 → #712 → #705 → #713 → #714 → #718 → #720 → #721 → #724 → #725 → #723 → #722.

---

Solves the following issue(s):

Solves #690 (the decision part), tracked in #650. Decides whether the guiding-center kernels (step 3 of the porting order in CUDA_STRATEGY.md) are ported by hand or generated from the pyccel source. To decide, one guiding-center kernel, bstar_parallel_3form, was ported both ways.

Decision: hand-written CUDA stays. The comparison table, the reasoning and what would change the decision are in the new CUDA_STRATEGY.md section Hand-written CUDA vs. code generation (decision for #690). Short version:

• The shared code is already ported. Most of what a generator would save is the helper chain: mappings, metric, splines, linear algebra. PRs 10–19 have already ported and tested it by hand. What is left per kernel is its body: here 36 lines of CUDA plus one 5-line helper.
• A generator would be a second compiler to maintain.
  • The numba prototype has 458 lines and covers only part of what the kernels use.
  • On all 60 pusher and accumulation entry kernels it translates 31. It fails on 2D slice assignments in the linear-algebra helpers and on numpy.abs in the SPH helpers.
  • Its mistakes are silent. It changes array lengths when a size is known only at run time, and it cannot tell which writes need atomics, so a translated accumulation kernel would race.
• numba does not fit the existing infrastructure.
  • numba kernels do not take the C structs. The prototype needs 38 flat parameters, plus a second launch path next to cunumpy's CudaKernel.
  • The parity tests and the CPU emulation compile the real .cu. numba's simulator only interprets the Python.
  • numba adds about 175 MB of dependencies (numba, llvmlite, numba-cuda) and supports only a range of NumPy versions.
• The other generators:
  • cupyx.jit would need the same rewriting of the pyccel source. Its reference lists no per-thread arrays and no struct arguments. Not run here (no CuPy).
  • pyccel 2.2.1 has no CUDA backend. The separate pyccel-cuda fork (last commit January 2025) still has open issues for thread indexing and passing arrays.
• What would change the decision:
  • pyccel releases a CUDA backend that accepts our argument classes. It could then generate the .cu files in place, behind the same tests.
  • Porting by hand becomes the bottleneck. Then the next step is a generator that writes .cu text (not numba), reusing the prototype's AST transformations.

Core changes:

• Hand-written port (a real port; it stays):
  • pic/pushing/kernels/bstar_parallel_3form/bstar_parallel_3form_cuda.cu, 1:1 with pyccel: same names and arguments, one thread per marker row, holes skipped as in pyccel (markers[ip, 0] == -1). numpy.mod(eta, 1.0) becomes eta - floor(eta) (Fortran MODULO, which pyccel generates).
  • It reuses df (every mapping), det, get_spans/SplineScratch and eval_spline_mpi_kernel.
  • New device helper eval_0form_spline_mpi in bsplines/evaluation_kernels_3d.cuh.
  • The pyccel kernel and its callers (KernelSetup in the guiding-center propagators) are unchanged.
• Parity cases: 14 cases in pic/tests/cuda_parity_cases.py.
  • One per mapping: the 10 analytic mappings plus IGAPolarCylinder, IGAPolarTorus, Tokamak and a 3d spline.
  • α = 1, 0 and mixed. Shifts make eta + shift wrap around in mod, row 0 has a negative eta_n and row 2 is a hole.
  • Two output columns; launch size args_markers.n_markers.
  • The CPU emulation agrees with pyccel to at most 5.7e-14. Removing the mod from the CUDA kernel makes the emulation test fail.
• Generated prototype (test-only): pic/tests/codegen_spike/ (numba_codegen.py, compare.py, test_numba_codegen.py).
  • It translates the pyccel source of the kernel and every helper it calls into one numba.cuda module (490 lines, 24 device functions).
  • With numba installed, the generated kernel agrees with pyccel on the same 14 cases (max 6.3e-13): in the numba.cuda simulator (NUMBA_ENABLE_CUDASIM=1) and compiled with numba.njit.
  • Not checked: whether numba.cuda compiles it for a GPU (needs NVVM).
  • Skipped without numba, which is not a struphy dependency. There is no kernels in the module names, so struphy compile ignores it.
  • It can be deleted once the decision is accepted.

Model-specific changes:

None. bstar_parallel_3form is one of the init kernels of PushGuidingCenterBxEstar and PushGuidingCenterParallel. Neither propagator runs on the GPU yet: the rest of their chain is not ported (driftkinetic_hamiltonian, unit_b_1form, bstar_2form, the discrete-gradient pushers, gc_density_0form).

Documentation changes:

CUDA_STRATEGY.md:

• New section: Hand-written CUDA vs. code generation (decision for Decide between hand-written CUDA and code generation before the guiding-centre kernels #690).
• The "Next" checklist entry, the step 3 note of the porting order and the open question on single-source alternatives point to the new section.

Testing (macOS, no GPU; pyccel kernels compiled with GNU/Fortran; GPU tests not run, so the new kernel has not run on a GPU):

	before (devel)	after
test_cuda_emulation.py, test_cuda_parity.py, test_kernel_backends.py, test_device_helpers.py	400 collected (83 passed, 317 skipped)	415 collected: 84 passed, 331 skipped (+1 emulation test with 14 cases, +14 GPU parity tests skipped)
mpirun -n 1 pytest test_kernel_setup.py test_pushers.py	181 passed	181 passed
codegen_spike/test_numba_codegen.py	–	2 skipped without numba; 2 passed with numba 0.68 and numba-cuda 0.24 installed

The "before" pass/skip split is derived from the collected tests: this PR adds exactly one emulation test and 14 GPU parity tests.

🤖 Generated with Claude Code