CUDA version of the vlasov_maxwell accumulation

[max-models]

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

Stack fix: simulation/tests/test_compile_cuda_kernels.py (from #711) now expects VlasovAmpereOneSpecies to have no kernel without a CUDA version (vlasov_maxwell has one now). The fail-fast check uses push_bxu_Hdiv instead. The CUDA_STRATEGY.md example was updated to match. The linear_vlasov_ampere (#705) and vlasov_maxwell parity cases and imports are both kept.

---

Solves the following issue(s):

Part of #688 (blocked matrix accumulations), tracked in #650. Ports vlasov_maxwell to CUDA. It is the accumulation of VlasovAmpereCoupling in VlasovAmpereOneSpecies and VlasovMaxwellOneSpecies (#3 in the porting order of CUDA_STRATEGY.md). It uses the Array6D<double> views of cunumpy 0.6.1.

This PR is based on devel and does not depend on #705 (linear_vlasov_ampere). The two PRs need the same device helpers. Here they are copied byte-identically from #705 (origin/cuda-linear-vlasov-ampere at 1ad4cd6), so filler_kernels.cuh, particle_to_mat_kernels.cuh, kernel_test_args.py, test_device_helpers.py and test_cuda_emulation.py merge without conflicts. Whichever PR merges second gets small conflicts in three files: cuda_parity_cases.py, test_cuda_parity.py and CUDA_STRATEGY.md. In all three, both sides add lines in the same place, so keep both (checked with git merge-tree).

Core changes:

• CUDA kernel pic/accumulation/kernels/vlasov_maxwell/vlasov_maxwell_cuda.cu:

  • Same arguments in the same order as the pyccel kernel, one thread per marker row.
  • Per marker it calls df, matrix_inv, transpose, matrix_matrix and matrix_vector. Then m_v_fill_b_v1_symm adds A_p = w_p DF^{-1} DF^{-T} and B_p = w_p DF^{-1} v_p.
  • The six blocks mat11 … mat33 are Array6D<double> (any strides) and the vectors are Array3D<double>.
  • As in pyccel, only holes are skipped. Boundary particles (markers[ip, -1] == -2) are accumulated, unlike in linear_vlasov_ampere.

• Device helpers, 1:1 with pyccel, byte-identical to CUDA version of the linear_vlasov_ampere accumulation #705:

  • fill_mat and fill_mat_vec in filler_kernels.cuh
  • m_v_fill_b_v1_symm in particle_to_mat_kernels.cuh

  Every addition into matrix or vector data is a cunumpy_atomic_add. outer from CUDA version of the linear_vlasov_ampere accumulation #705 is not needed here and is not copied.

• Tests:

  • PARITY_CASES["vlasov_maxwell"] covers 129 markers in five mappings: Cuboid, Colella, HollowTorus, ShafranovDshapedCylinder and a 3d spline mapping (geometry_domain 0, 2, 5, 9, 12).
    • Row 0 is a real hole (markers[0, 0] = -1) and row 1 is a boundary particle.
    • The data are shaped like the real accumulation matrices (v1_symm_accumulation_data() from CUDA version of the linear_vlasov_ampere accumulation #705).
  • The CPU emulation agrees with pyccel to 7e-12 or better. Entries reach 6e5 in the spline mapping, because of G^{-1}.
  • Tolerances are rtol=1e-12, atol=1e-8, for the atomic summation order on the GPU. atol is scaled up from linear_vlasov_ampere's 1e-10 because the entries are about 60 times larger.
  • Mutation checks: the emulated parity test fails for a transposed DF, for skipped boundary particles and for holes that are not skipped.
  • vlasov_maxwell is added to test_vlasov_kernel_coverage.
  • From CUDA version of the linear_vlasov_ampere accumulation #705 (identical): test_v1_symm_filler (GPU) and test_emulated_v1_symm_filler compare m_v_fill_b_v1_symm with pyccel.

• Not in this PR: CUDA version of the linear_vlasov_ampere accumulation #705's test_accum_matrix_cupy.py and emulated_launches(). They test the matrix path of Accumulator on CuPy and are not copied, to keep this PR small. I ran an ad-hoc copy of that test with vlasov_maxwell in place of linear_vlasov_ampere, not committed. It runs on the fake CuPy with emulated launches and matches NumPy: the nine kernel arrays, the dense blocks, the vector, BC.dot(x) and three Schur CG iterations, with no host copies from the accumulation to BC.dot(x). After both PRs are merged, that test could get a vlasov_maxwell case.

Model-specific changes:

None. Status of VlasovAmpereOneSpecies and VlasovMaxwellOneSpecies on the GPU:

Part	Kernels / operations	GPU status
PushEta	push_eta_stage, reflect	ported
PushVxB (analytic and implicit)	push_vxb_analytic, push_vxb_implicit	ported
VlasovAmpereCoupling	vlasov_maxwell	this PR
	push_v_with_efield	ported
	Schur solve	needs the FEEC parts below
Initial PoissonSolve	charge_density_0form	ported
	weights through the control variate (update_weights, Maxwellian on xp, jacobian_det through the geometry kernels)	ported
	stiffness solve	needs the FEEC parts below
MaxwellWeakAmpere (VlasovMaxwell only)	FEEC only (curl, M1, M2, Schur solve)	needs the FEEC parts below
Scalars	KineticEnergyPIC (xp.sum)	runs on the device
	BilinearEnergyFEEC	goes through StencilVector.inner, so a host scalar
	gauss_error diagnostic	uses toarray(), a host copy

Still missing for an end-to-end GPU run:

• Preconditioner: MassMatrixPreconditioner, the default of VlasovAmpereCoupling, MaxwellWeakAmpere and PoissonSolve, cannot be created on CuPy (is_circulant asserts a host array). The device Kronecker solve is a separate PR.
• CG inner products: feectools' StencilVector.inner returns a host scalar, so each CG iteration and the FEEC energy scalars copy a scalar from the device.
• FEEC operators: the derham operators and the mass matrices on the device need the feectools CUDA stack (stencil 6D views, device Kronecker solve, PR 17).
• Marker sorting: multi-rank marker sorting still goes through the host.
• H100 run: the first run on an H100 (First GPU run of the CUDA tests (H100) #687), then the models end to end (Run VlasovAmpereOneSpecies end to end on the GPU #689).

Documentation changes:

CUDA_STRATEGY.md:

• Porting order: Add Pylint workflow for Python code analysis #3 vlasov_maxwell is marked ✓.
• New short section: vlasov_maxwell implementation notes.

The checklist and the "6D views are a gate" paragraph are left to #705, to keep the conflict small.

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

	before (devel)	after
test_kernel_backends.py test_cuda_parity.py test_device_helpers.py test_cuda_emulation.py	83 passed, 317 skipped	85 passed, 323 skipped
mpirun -n 2: test_verif_VlasovMaxwellOneSpecies.py test_verif_VlasovAmpereOneSpecies.py test_accum_vec_H1.py	–	7 passed

• The 2 new passes are test_emulated_parity[vlasov_maxwell] and test_emulated_v1_symm_filler.
• The 6 new skips are GPU tests: 5 vlasov_maxwell parity cases and test_v1_symm_filler.
• The MPI tests only use the pyccel kernels, which this PR does not change. I ran them only on this branch.

🤖 Generated with Claude Code