CUDA version of the linear_vlasov_ampere accumulation

[max-models]

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

---

Solves the following issue(s):

Part of #688 (blocked matrix accumulations), tracked in #650. First matrix accumulation on the GPU: linear_vlasov_ampere, the accumulation of EfieldWeightsCoupling in LinearVlasovAmpereOneSpecies (step 2 of the porting order in CUDA_STRATEGY.md). It uses the Array6D<double> views of cunumpy 0.6.1.

Core changes:

• CUDA kernel pic/accumulation/kernels/linear_vlasov_ampere/linear_vlasov_ampere_cuda.cu: same arguments in the same order as the pyccel kernel. It launches one thread per marker row and skips holes and boundary particles. Per marker it calls df, matrix_inv, matrix_vector, outer and then m_v_fill_b_v1_symm. The six blocks mat11 … mat33 are Array6D<double> (strided), the vectors are Array3D<double>, and f0_values is const double*, as in push_weights_with_efield_lin_va.

• Device helpers, 1:1 with pyccel (same names, arguments and loop order):

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

  Every addition into matrix or vector data is a cunumpy_atomic_add, as in charge_density_0form. The spline values live in the thread-local SplineScratch.

• Accumulator matrix path on CuPy: no code change was needed, and the new test checks it. On CuPy, StencilMatrix._data/StencilVector._data are device arrays. The accumulator zeroes them in place and passes the same arrays to the kernel. The ghost-region exchange, the transposed blocks of the symmetric matrix (feectools' device transpose), the scaled operator built by EfieldWeightsCoupling for SchurSolver, its dot and the Schur solve all run on them.

• Tests:

  • PARITY_CASES["linear_vlasov_ampere"] covers 129 markers (with a hole and a boundary particle) in Cuboid, Colella, HollowTorus and a 3d spline mapping. The data are shaped like the real accumulation matrices (v1_symm_accumulation_data() in kernel_test_args.py). Tolerances are rtol=1e-12, atol=1e-10: atomics change the summation order and entries reach 1e4. The CPU emulation catches a transposed DF.
  • m_v_fill_b_v1_symm is compared with pyccel through a wrapper kernel: test_v1_symm_filler on a GPU, test_emulated_v1_symm_filler by CPU emulation.
  • New test_accum_matrix_cupy.py accumulates as EfieldWeightsCoupling does (Colella, clamped in eta1) and compares with NumPy:
    • It checks the nine kernel arrays, the dense blocks, the accumulated vector, BC.dot(x) and three CG iterations of the Schur solve.
    • It also checks that the kernel writes into the blocks' own _data (no copies), and that nothing is copied to the host from the accumulation up to BC.dot(x).
    • On a GPU it runs with assert_no_transfers. Without a GPU it runs in a subprocess on cunumpy's fake CuPy (CUNUMPY_FAKE_CUPY=1, clean MPI environment). There, every CudaKernel launch is emulated on the CPU in place on the fake device arrays: the new emulated_launches() in pic/tests/cuda_emulation.py covers struphy's and feectools' kernels, struct arguments included. This test also fails on a transposed DF.

Model-specific changes:

None. LinearVlasovAmpereOneSpecies cannot run end to end on the GPU yet. Still needed:

• CG inner products: feectools' StencilVector.inner returns its result on the host (xp.to_numpy), so each CG inner product copies a scalar. The Schur solve is outside the transfer guard because of this.
• Default preconditioner: MassMatrixPreconditioner (the default of EfieldWeightsCoupling) cannot be created on CuPy, because is_circulant asserts a host array. The device Kronecker solve is a separate PR.
• Parallel feectools PRs: the stencil 6D views and the device Kronecker solve.
• H100 run: the first run on an H100 (First GPU run of the CUDA tests (H100) #687), then the model end to end (Run VlasovAmpereOneSpecies end to end on the GPU #689).

Documentation changes:

CUDA_STRATEGY.md:

• linear_vlasov_ampere is marked as ported in the checklist and the porting order.
• 6D views are no longer a gate.
• New section: implementation notes for linear_vlasov_ampere.

Fix added after review (commit 1979ea9): the test hole was not a hole for accumulation kernels.
kernel_test_args.marker_arguments() marked the hole (row 0) only in column 8 (first_init_idx), which the pushers test.
Accumulation kernels, including linear_vlasov_ampere, test column 0, so this PR's parity case had no real hole.
Row 0 is now -1 in every column, as struphy marks holes. After the fix: test_cuda_emulation.py + test_accum_matrix_cupy.py 72 passed, 1 skipped; test_cuda_parity.py -k linear_vlasov_ampere passed on NumPy.

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, 322 skipped
test_accum_matrix_cupy.py (new)	–	1 passed (fake CuPy), 1 skipped (GPU)
mpirun -n 2: test_accum_vec_H1.py test_pushers.py test_verif_LinearVlasovAmpereOneSpecies.py	110 passed	110 passed

The new skips are GPU tests: 4 parity cases, test_v1_symm_filler and test_accumulator_matrix_on_cupy.

🤖 Generated with Claude Code