Field evaluation and diagnostics on the device

[max-models]

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

Stack fixes (for #721 below): SplineFunction.eval_tp_fixed_loc now takes kind, pn and starts from self._args_spline[i], because the tuple _args_eval no longer exists. The parity case of eval_spline_mpi_tensor_product_fixed unpacks spline_evaluation_arguments the same way. test_catalog_signatures now counts 4 spline evaluation kernels. Follow-up, not done here: switch eval_spline_mpi_tensor_product_fixed to the SplineArguments class of #721.

---

Solves the following issue(s):

Part of #697 (field evaluation and diagnostics on the device), tracked in #650. On the CuPy backend, the diagnostics of each output step now run on the device. Only the reduced results reach the host: the scalars in one copy per update, and every other saved dataset in one counted copy per output step. The marker diagnostics kernels in pic/diagnostics/kernels (used by the 5D and hybrid models) are not ported yet. They are listed as follow-up work, so this PR does not close #697.

Core changes:

I found the problems below without a GPU, on cunumpy's fake CuPy (CUNUMPY_FAKE_CUPY=1, clean environment) with cunumpy.profiling.count_transfers. CUDA launches were emulated on the CPU with emulated_launches() from #705's worktree, in ad-hoc probes that are not committed. In those probes, implicit syncs (float(), bool(), __index__, .get() of a device array) were counted by patching the fake ndarray.

Path	Before (fake CuPy)	After
Scalars.update (FunctionScalarPIC etc.)	one float() per function scalar (implicit sync)	no syncs; one xp.to_numpy of all values (8 bytes per scalar)
DataContainer.save_data	one raw .get() per device dataset, not seen by count_transfers; scalars and time copied one by one	one counted xp.to_numpy per device dataset; scalars and time are already on the host, so they need no copy
Particles.binning	fails: xp.count_nonzero(list)	runs xp.histogramdd on the device, no transfers
update_markers_to_be_saved	xp.count_nonzero used as a slice bound (__index__ sync)	the row count is taken from the gather; no extra sync
SplineFunction.__call__ (meshgrid/point)	six bool() syncs per call (process bounds compared as device scalars)	bounds cached on the host once; flagging uses xp.copyto(..., where=...)
SplineFunction.eval_tp_fixed_loc	fails: pyccel kernel called with device arrays; xp.array(degree) uploaded per call	CUDA kernel, no transfers
Derham.prepare_eval_tp_fixed (setup)	fails on CuPy (pyccel helpers given device buffers)	computed on the host, copied to the device once

• Scalars (models/scalars.py): Scalars keeps one device buffer. Each contained scalar's value is a view into it, so update() ends with to_host(), a single xp.to_numpy. host_value(name) returns a NumPy view of the host copy. FunctionScalarFEEC/PIC/SPH store the function's result through _scalar_value instead of float(), so a device reduction stays on the device. StruphyModel.print_scalar_quantities and the HDF5 setup read the host copy.
• Output (io/output_handling.py): DataContainer._as_numpy_array is now xp.to_numpy. This is the one place where output copies device data.
• Time state (simulation/sim.py): time_state is NumPy on both backends. The loop reads it with float() every step, which was a device sync on CuPy. TimeDependentSource still gets a scalar it can pass to xp.sin/xp.cos.
• Binning (pic/base.py): the valid markers are gathered once, and the Jacobian is evaluated once instead of twice. On NumPy the results are the same values as before.
• Kernel eval_spline_mpi_tensor_product_fixed: moved from bsplines/evaluation_kernels_3d.py into bsplines/kernels/eval_spline_mpi_tensor_product_fixed/ (folder rule). It has a new _cuda.cu (one thread per grid point, reusing eval_spline_mpi_kernel) and parity cases in pic/tests/cuda_parity_cases.py (additive; CUDA version of the linear_vlasov_ampere accumulation #705/CUDA version of the vlasov_maxwell accumulation #713 also add entries there). eval_tp_fixed_loc takes kind/pn/starts from SplineFunction._args_eval instead of building them for every call.
• Not changed, still implicit syncs: boolean-mask gathers (markers[valid_mks], also inside KineticEnergyPIC), one per gather. BilinearEnergyFEEC/VolumeFormEnergyFEEC still copy one scalar per inner-product block inside feectools until feectools: inner products stay on the device #716 (Inner products stay on the device feectools#98) is in. With feectools: inner products stay on the device #716 these become device scalars, and this PR's container already handles them.

Model-specific changes:

None. LinearVlasovAmpereOneSpecies._compute_en_w now stays on the device because FunctionScalarPIC no longer calls float().

Documentation changes:

CUDA_STRATEGY.md: one checklist entry and a short section "Diagnostics on the device (#697)" with the follow-ups: CUDA versions of the marker diagnostics kernels in pic/diagnostics/kernels/, and SPH kernel-density plots (step 6).

Tests:

• New src/struphy/models/tests/test_diagnostics_on_device.py has five checks (scalars, binning, saved markers, output, point flagging) with mock particles. They need no CUDA kernel. Each check runs:

  • on NumPy (test_on_numpy);
  • on fake CuPy in a subprocess (test_on_fake_cupy), compared with NumPy, counting transfers, with implicit syncs turned into errors;
  • on a real GPU (test_on_cupy, skipped here).

  test_eval_tp_fixed_loc_on_cupy (GPU, skipped here) compares the CUDA eval_tp_fixed_loc with NumPy.

• New parity cases eval_spline_mpi_tensor_product_fixed (2 cases), run by CPU emulation in test_cuda_emulation.py.

Local tests run (macOS, no GPU, pyccel kernels compiled with GNU/Fortran; GPU tests not run):

• models/tests/test_diagnostics_on_device.py, io/tests/test_output_handling.py, pic/tests/test_binning.py, models/tests/test_scalars.py: 24 passed, 6 skipped (GPU), 2 failed. The two failures are test_en_fB_ignores_ghost_markers and test_guiding_center_en_fB_does_not_divide_by_Np in test_scalars.py. Their mocks have no save_magnetic_background_energy, which GuidingCenter._compute_en_fB calls, and the failures come from code this PR does not touch, so I think they also fail on devel. I did not confirm this with a separate devel run.
• pic/tests/test_cuda_emulation.py -k tensor_product and the parity-case/signature checks of test_cuda_parity.py: passed (GPU parity skipped).
• models/tests/verification/test_verif_LinearVlasovAmpereOneSpecies.py: 2 passed, serial and with mpirun -n 2 ... --with-mpi.
• Not run locally: the full suite, the other model tests, and test_variational_dissipation.py (the models that use eval_tp_fixed_loc). CI covers these.

🤖 Generated with Claude Code