Repository navigation
Declare kernel outputs for PyccelKernels - #723
Merged
Merged
Conversation
Each kernel folder's __init__.py declares the arguments its kernel writes to
(OUTPUTS, passed to cunumpy as host_options={"outputs": OUTPUTS}); the
PyccelKernel wraps in feec/mass.py, boundary_mass.py and
basis_projection_ops.py pass outputs= too. The parity tests compare the
declared outputs; new tests check the declarations on the NumPy backend
and in the CPU emulation of the CUDA kernels.
Solves #342.
Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
Stack the PR on #725. Conflicts: feec/mass.py, boundary_mass.py, basis_projection_ops.py (took #725's kernel folders; their outputs moved into the folders), test_cuda_emulation.py (#723's argument_arrays check on top of the tests added below), CUDA_STRATEGY.md (kept both notes). Semantic fixes for the PRs below: - OUTPUTS for the 16 FEEC assembly folders of #725 (same arguments as the former PyccelKernel(outputs=...) wraps, plus kernel_*d_vec and surface_kernel_3d_vec) and for eval_spline_mpi_tensor_product_fixed (#724). - OUTPUTS of eval_spline_mpi_markers/matrix/sparse_meshgrid point at values in the SplineArguments signatures of #721 (3, 5, 5). - argument_arrays (cuda_parity_cases.py) knows CudaSplineArguments (#721). Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
Pulls in the two follow-up fixes of #724 (already resolved identically in the previous merge). Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
max-models
added a commit
that referenced
this pull request
Oct 8, 2026
Stack the PR on #723. Conflicts: eval_spline_mpi_{markers,matrix,sparse_meshgrid}_cuda.cu (SplineArgs of #721 with the qualified evaluation_kernels_3d::eval_spline_mpi), linalg_kernels.cuh, filler_kernels.cuh, particle_to_mat_kernels.cuh (helpers of #705/#713 moved into their module namespaces). Semantic fixes: the struphy_cuda::<pyccel module> rule applied to the CUDA code of the PRs below: linear_vlasov_ampere (#705), vlasov_maxwell (#713), bstar_parallel_3form (#718) and eval_spline_mpi_tensor_product_fixed (#724) use 'using namespace struphy_cuda;' and qualified helper calls; m_v_fill_b_v1_symm calls filler_kernels:: and evaluation_kernels_3d::; the fill_v1_symm test wrapper calls particle_to_mat_kernels::m_v_fill_b_v1_symm. Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
This was referenced Oct 8, 2026
max-models
added this pull request to stack #728
October 8, 2026 05:45
max-models
removed this pull request from stack #728
October 8, 2026 06:00
max-models
added a commit
that referenced
this pull request
Oct 8, 2026
max-models
marked this pull request as draft
October 8, 2026 09:13
max-models
added a commit
that referenced
this pull request
Oct 8, 2026
Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
max-models
added a commit
that referenced
this pull request
Oct 12, 2026
**Stack:** part 1 of 14, based on `devel`, next: #708. Full order: #709 → #708 → #711 → #712 → #705 → #713 → #714 → #718 → #720 → #721 → #724 → #725 → #723 → #722. --- Two small bugs that stop particle runs early on `CUNUMPY_BACKEND=cupy`: - **`SplineFunction.__init__`** (`feec/psydac_derham.py`): the knots of `derham.V0fem` are NumPy arrays, and `xp.ascontiguousarray` (= `cupy.ascontiguousarray`) does not accept them: `TypeError: Argument 'a' has incorrect type (expected cupy._core.core._ndarray_base, got numpy.ndarray)`. They are now made contiguous on the host and copied to the device with `xp.asarray`. - **`SortingBoxes._set_boxes`** (`pic/sorting.py`): `xp.sqrt` of a Python int returned a 0-d CuPy array, which `round()` rejects (`TypeError: type ndarray doesn't define __round__ method`). It is a host scalar, so this now uses `math.sqrt`. The NumPy backend is unaffected. ### Note With these fixes, `examples/VlasovAmpereOneSpecies/strong_Landau_damping` gets further on CuPy but still stops in `initialize_neighbours`: box sorting and the `vlasov_maxwell` accumulation are not on the GPU yet (see `CUDA_STRATEGY.md`, "Next"). This PR does not try to change that.
max-models
added a commit
that referenced
this pull request
Oct 12, 2026
**Stack:** part 2 of 14, based on #709 (merge that first), next: #711. Full order: #709 → #708 → #711 → #712 → #705 → #713 → #714 → #718 → #720 → #721 → #724 → #725 → #723 → #722. --- ## Fuse the pusher preamble into one CUDA kernel On the CuPy backend, the three column-slice assignments at the start of `Pusher._push` ```python markers[:, first_pusher_idx:first_shift_idx] = markers[:, : 3 + vdim] markers[:, first_shift_idx:residual_idx] = 0.0 markers[:, residual_idx:-2] = 0.0 ``` each sweep the whole row-major marker buffer, because the slices are strided. With 4M markers (an 8M-row buffer of 1.66 GB), this took 4.6 ms per pusher call, more than the push kernels themselves. ## Changes - New src/struphy/pic/pushing/prepare_push_cuda.cu: one thread per entry of the contiguous column range [first_pusher_idx, n_cols - 2), so all three assignments are done in a single coalesced pass. - Pusher._push calls it on the CuPy backend; the kernel is compiled once, on first use. The NumPy path is unchanged. This roughly speeds up the preamble by 2x. The next big improvement for the markers would be #707 --------- Co-authored-by: Claude Opus 5.5 <noreply@anthropic.com>
max-models
added a commit
that referenced
this pull request
Oct 12, 2026
…kernels) (#711) **Stack:** part 3 of 14, based on #708 (merge that first), next: #712. Full order: #709 → #708 → #711 → #712 → #705 → #713 → #714 → #718 → #720 → #721 → #724 → #725 → #723 → #722. --- **Solves the following issue(s):** Solves #704, tracked in #650. Until now a CUDA kernel compiled on its first call, which put the compile time into the first time step and its profiling regions. Now `Simulation.run()` calls the new `Simulation.compile_cuda_kernels()` after `allocate()`, so the kernels compile before the time stepping. **Core changes:** - **`Simulation.compile_cuda_kernels()`** (`simulation/sim.py`): - It runs after `allocate()` in `run()`, inside a new profiling region, `setup: compile cuda kernels`. - On the CuPy backend it calls `Kernel.compile()` once for each kernel of `Simulation.kernels()` and returns the compiled kernels. - On the NumPy backend it does nothing and returns `()`. NumPy runs are unaffected. - **Fail fast:** on CuPy, if any kernel the time stepping calls has no CUDA version, the method raises `NotImplementedError` before compiling anything. - The error lists every such kernel and the path where its `.cu` file is expected. For example, `VlasovAmpereOneSpecies` names `vlasov_maxwell`. - Without this check, the run fails at the first call of the first missing kernel, in the time loop. - Kernels with `missing_cuda="fallback"` are not counted as missing. - **`Simulation.kernels()`:** collects the kernels from the objects that keep and call them, through a new `kernels()` method on each owner. Each kernel is listed once (`utils/kernel_compilation.collect_kernels`). The owners: - `Domain`: the four geometry kernels. - `Derham`: the three spline evaluation kernels of `SplineFunction`, and feectools' `stencil_dot_3d`, `stencil_transpose_3d`, `stencil_inner_3d` and `stencil_axpy_3d` (called by the FEEC operators and solvers). - `Particles`: `reflect`, if an axis has reflecting boundary conditions. - `Pusher`: the pusher kernel and the kernels of its init/eval `KernelSetup`s. `KernelSetup`, `Accumulator` and `AccumulatorVector` list their own kernel. - `Propagator.kernels()`: by default it collects from the propagator's attributes. That covers `Kernel`s kept as attributes and objects with a `kernels()` method, also inside lists, tuples and dicts. It does not go into other objects. A propagator that calls a kernel it does not keep can override the method; none needs to today. - **Not listed:** kernels that only diagnostics or the initial solves call. They run during setup ("setup: initial diagnostics", "model post allocate") anyway. **Model-specific changes:** None. I checked the kernel lists against what the time loop actually calls. For each of 19 models (all particle models, plus Maxwell and LinearMHD), I ran the default parameter file for one step on NumPy while recording every `Kernel` call during `model.integrate`. Every kernel called in the time loop is in `Simulation.kernels()` for all 19 models. (This was an ad-hoc script; the committed test does the same check for `VlasovAmpereOneSpecies`.) **Documentation changes:** `CUDA_STRATEGY.md` has a short new section, "Compiling the CUDA kernels at setup (struphy#704)". The PR 16 notes now point to it from the sentence "when to compile them is decided later". **Tests** (new `simulation/tests/test_compile_cuda_kernels.py`): - `collect_kernels` walks owners and containers, lists each kernel once, and `missing_cuda` finds the kernels without a CUDA version. - For a small `VlasovAmpereOneSpecies`, `Simulation.kernels()` is exactly the domain + Derham kernels plus `push_eta_stage`, `vlasov_maxwell` and `push_v_with_efield`, and the only kernel missing CUDA is `vlasov_maxwell`. Reflecting particles add `reflect`. - Every `Kernel` called in the time loop of a one-step run is listed. - On NumPy, `compile_cuda_kernels()` is a no-op (`Kernel.compile` patched to fail). - With the backend switched to CuPy by monkeypatching: every kernel is compiled, and a missing CUDA version raises before anything is compiled. - The same CuPy path runs in a subprocess on cunumpy's fake CuPy (`CUNUMPY_FAKE_CUPY=1`, clean MPI environment), with `CudaKernel.compile` recorded. - GPU test (skipped here): `Vlasov` compiles all its kernels at setup, and `VlasovAmpereOneSpecies` raises naming `vlasov_maxwell`. **Testing** (macOS, no GPU; pyccel kernels compiled with GNU/Fortran; **GPU tests not run**): The suites were `simulation/tests`, `pic/tests/{test_kernel_setup,test_reflect,test_pushers,test_kernel_backends,test_cuda_emulation}.py`, `feec/tests/test_derham_gpu.py` and `geometry/tests/test_domain.py`. | | before (`devel`) | after | | --- | --- | --- | | those suites | 336 passed, 55 skipped | 343 passed, 56 skipped (7 new passed, the new GPU test skipped) | The "before" numbers are the same suites without the new test file; no existing test file was changed. Not verified without a GPU: the real NVRTC compile through `Kernel.compile()` and the GPU test. On the CPU, only the call path up to `CudaKernel.compile` is tested. 🤖 Generated with [Claude Code](https://claude.com/claude-code) --------- Co-authored-by: Claude Opus 5.5 <noreply@anthropic.com>
max-models
added a commit
that referenced
this pull request
Oct 12, 2026
**Stack:** part 4 of 14, based on #711 (merge that first), next: #705. Full order: #709 → #708 → #711 → #712 → #705 → #713 → #714 → #718 → #720 → #721 → #724 → #725 → #723 → #722. --- **Solves the following issue(s):** Solves #698, part of #650. Markers move between ranks without going through the host on the CuPy backend: `Particles.mpi_sort_markers` exchanges the marker rows as device buffers with CUDA-aware MPI. **Core changes:** - **Device exchange of marker rows** (`Particles._sendrecv_markers_device` in `pic/base.py`, called by `_sendrecv_markers` on CuPy): - The rows to send (already gathered on the device per destination by `_sendrecv_get_destinations`) and one device receive buffer go to `Isend`/`Irecv` through `cunumpy.mpi.mpi_buffer`. With CUDA-aware MPI this passes the device arrays directly, after `cunumpy.mpi.synchronize_for_mpi`, as feectools' data exchangers have done since struphy-hub/feectools#86. Without CUDA-aware MPI the same buffers are staged through (pinned) host memory, with a one-time warning, instead of handing device pointers to a host-only MPI. - Rows from rank `i` arrive in a contiguous slice of the receive buffer. The holes are filled in the same order as on NumPy, but by one device scatter instead of one per source rank. - Messages with a count of zero are skipped. Both sides know every count after the `Alltoall`, so they skip the same messages. - Whether MPI is CUDA-aware: the answer already recorded by `xp.mpi.mpi_is_cuda_aware()`/`set_mpi_cuda_aware()`, otherwise a collective probe of `mpi_comm` at the first exchange. This is safe because `mpi_sort_markers` is collective. The probe copies 32 bytes to the host once. - **Counts on the host without copies:** `send_info`/`recv_info` are small NumPy arrays on every backend. The send counts are sizes of `nonzero` results, which are host integers already (CuPy synchronizes to size them). So the `Alltoall` of the counts passes host buffers, and the receive buffer, hole slices and load-imbalance check need no device reads. Before this PR, CuPy used device count arrays: one host write per rank, `Alltoall` on device buffers, and `list(recv_info)`, `recv_info[i]` and `first_hole[i]` read back element by element. - **`alpha`:** `_sendrecv_determine_mtbs` checks a tuple or list `alpha` on the host. Before, it asserted on device comparisons, which made two implicit reads. A uniform `alpha` (the usual case: 1.0 or 0.5 everywhere) becomes a Python scalar, so the 3-vector is no longer copied to the device. The values are the same on NumPy. - **NumPy path unchanged:** the same messages in the same order, and the same hole filling. The only differences are the NumPy count arrays, which `xp` already gave on NumPy, and the scalar `alpha`, which gives the same values. - **Tests:** new `pic/tests/test_sorting_device.py`: - `test_mpi_sort_markers_fake_cupy[bc, cuda_aware]` runs under `mpirun -n N env CUNUMPY_FAKE_CUPY=1 pytest --with-mpi ...` and is skipped unless the fake CuPy is active. It sorts the same markers on NumPy and on the fake CuPy across all ranks, two shifts in a row, so most markers change rank, with periodic and remove boundaries. It compares the whole marker arrays and `valid_mks` exactly, then sorts again with `do_test=True`. - With CUDA-aware MPI (the probe, or forced), it checks that `count_transfers()` is 0 during `mpi_sort_markers` and that no fake CuPy array is read by the host (`int`/`bool`/`float`/`__index__`/iteration/`.get`/`.item`/`.tolist` are patched to count). Before this PR there were 4 such reads per call. - With staging (`set_mpi_cuda_aware(False)`), it checks for exactly one `to_host` per non-empty send buffer, carrying `n_send * n_cols * 8` bytes, and one `to_device` when something is received. - Fake CuPy arrays are host memory, so the non-CUDA-aware Open MPI here takes them through `__cuda_array_interface__`. - `test_mpi_sort_markers_gpu[bc]` (`requires_cupy`, `mpirun -n 2` with one GPU per rank): the same check with the probed MPI, without the host-read patching. - CI: a step in `test-PR-unit-mpi.yml` runs the fake CuPy test with the same `n-procs` matrix (1–4). **Model-specific changes:** None. Not covered by this PR: - The SPH ghost-box exchange (`_communicate_boxes` → `_sendrecv_markers_boxes`) still passes device arrays to MPI without `synchronize_for_mpi`, and it reads `_recv_info_box` element by element. It can use the same pattern in a separate PR. - `apply_kinetic_bc`, which `mpi_sort_markers` calls first, was not changed. With the boundaries tested here, it makes no host/device copy and no host read. - A two-rank GPU run (the last item of #698) has not been done: no GPU here. **Documentation changes:** - `CUDA_STRATEGY.md`: the "MPI + GPUs" open question now refers to the new section "Marker exchange implementation notes (#698)". - Docstrings of `mpi_sort_markers`, `_sendrecv_all_to_all` and the new `_sendrecv_markers_device`/`_mpi_cuda_aware`. **Testing** (macOS, Open MPI 5.0.9 (not CUDA-aware), pyccel kernels compiled with GNU/Fortran; **GPU tests not run**): | | before (`devel`) | after | |---|---|---| | `mpirun -n 2 pytest --with-mpi` `test_sorting.py test_pushers.py test_draw_parallel.py test_neighbor_ranks.py test_accum_vec_H1.py` | 205 passed | 205 passed | | same, `mpirun -n 3` | 205 passed | 205 passed | | `mpirun -n 2 pytest --with-mpi test_sph.py test_kernel_setup.py test_tesselation.py` | – | 743 passed | | `mpirun -n 2 env CUNUMPY_FAKE_CUPY=1 pytest --with-mpi test_sorting_device.py` | 6 failed, 2 skipped | 6 passed, 2 skipped (GPU) | | same with `-n 1`, `-n 3`, `-n 4` | – | 6 passed, 2 skipped | On `devel`, the fake CuPy test fails because of the implicit host reads (`__bool__` ×3, `__iter__`) with CUDA-aware MPI, and because the device buffers are never staged when MPI is not CUDA-aware. 🤖 Generated with [Claude Code](https://claude.com/claude-code) --------- Co-authored-by: Claude Opus 5.5 <noreply@anthropic.com>
max-models
added a commit
that referenced
this pull request
Oct 12, 2026
**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 (#687), then the model end to end (#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](https://claude.com/claude-code) --------- Co-authored-by: Claude Opus 5.5 <noreply@anthropic.com>
max-models
added a commit
that referenced
this pull request
Oct 12, 2026
**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 #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 #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 #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 #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:** #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 (#687), then the models end to end (#689). **Documentation changes:** `CUDA_STRATEGY.md`: - Porting order: #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](https://claude.com/claude-code) --------- Co-authored-by: Claude Opus 5.5 <noreply@anthropic.com>
max-models
added a commit
that referenced
this pull request
Oct 12, 2026
**Stack:** part 7 of 14, based on #713 (merge that first), next: #718. Full order: #709 → #708 → #711 → #712 → #705 → #713 → #714 → #718 → #720 → #721 → #724 → #725 → #723 → #722. --- **Solves the following issue(s):** Solves #696, part of #650 (see `CUDA_STRATEGY.md`). Every equilibrium in `fields_background/equils.py` can now be created and evaluated on the CuPy backend. The time loop does not change: it uses the projected equilibrium. **Core changes:** - **What failed on CuPy.** I checked this with cunumpy's fake CuPy in a clean serial subprocess. Every method models call was evaluated on a meshgrid and at markers. - **Construction failed** for: - `AdhocTorus` with `q_kind` 1 or 2: SciPy `quad`/`UnivariateSpline` got device arrays. - `AdhocTorusQPsi`: `odeint`/`fsolve` got device arrays. - `EQDSKequilibrium`: `RectBivariateSpline` got the device arrays from `xp.linspace`. - **Evaluation failed** for every field of `GVECequilibrium` and `DESCequilibrium`: gvec and DESC/JAX got device arrays. - **Already working:** the analytic ones (`HomogenSlab`, `ShearedSlab`, `ShearFluid`, `ScrewPinch`, `AdhocTorus` with `q_kind=0`, `CircularTokamak`, `ConstantVelocity`, `HomogenSlabITG`, `CurrentSheet`, the generic ones). - **Also found:** `AdhocTorus.psi_r` (`q_kind=0`) did `out *= ...` on a 0-d device array, which fails on CuPy when `psi` is called with arrays. It is now `out = out * ...`. - **New helpers in `fields_background/base.py`:** - `setup_on_host`: a decorator for `__init__`. It runs the setup on the NumPy backend (file reading, ODE solves, SciPy fits), so these equilibria hold only host data on either backend. Used on `AdhocTorus`, `AdhocTorusQPsi` and `EQDSKequilibrium`. - `host_call(fun, *args, **kwargs)`: calls a host-only function, copying device arguments to the host and array results back once per call. Without device arguments it is a plain call, so the NumPy path does not change. Used for every SciPy spline evaluation: `psi_r`/`p_r` of `AdhocTorus`, `psi_r` of `AdhocTorusQPsi`, and `psi`/`q_psi`/`g_psi`/`p_psi` of EQDSK. - `evaluate_on_host`: the same as a method decorator. Used on the GVEC `bv`, `jv`, `p0`, `n0` and the DESC `bv`, `jv`, `p0`, `n0`, `gradB1`. - **Cost of the host evaluation.** Each call costs one host round trip. For example, `gradB_xyz` of EQDSK makes about 10 `psi` calls, so 10 round trips. This is fine because equilibria are only evaluated at setup: initial conditions, projections onto the FEEC spaces and mappings. A device B-spline evaluation of the SciPy knots and coefficients would remove these copies. It is listed as an open question in `CUDA_STRATEGY.md`. - **`Tokamak`:** the workaround from PR 19 is removed. It no longer builds its default `EQDSKequilibrium` inside `use_backend("numpy")`. The field-line tracing (SciPy) still runs on the NumPy backend. - **Tests:** new file `fields_background/tests/test_equils_cupy.py`. - **What it checks:** 19 cases that cover every equilibrium class of `equils`, including EQDSK, GVEC and DESC with the files shipped in the repo or with the packages. `test_all_equilibria_have_cases` checks that no class is missing. Each case is created on NumPy and on CuPy, and 30 logical-domain methods (`absB0`, `p0`, `n0`, `b2`, `unit_b1`, `gradB1`, `curl_unit_b2`, `j2`, …) are evaluated: - on a meshgrid and at markers, - plus `psi`/`g_tor` with all derivatives for the axisymmetric equilibria. Every array result must be a device array and match NumPy to 1e-12. - **`test_equil_fake_cupy`** runs without a GPU, in a subprocess on the fake CuPy. The fake CuPy cannot launch CUDA kernels, so `host_geometry_kernels()` runs the four geometry entry kernels with their pyccel version on the host buffers of the fake arrays. The CUDA geometry kernels have their own parity and emulation tests. The `emulated_launches()` of #705 would also work here, but it compiles every launch, and these cases make hundreds of domain evaluations. - **`test_equil_cupy`** runs the same check on a GPU. **Model-specific changes:** None. **Documentation changes:** `CUDA_STRATEGY.md`: - New section "MHD equilibria on CuPy (#696)". - The open question "Polar splines and MHD equilibria" now covers polar splines only. - New open question: equilibrium splines on the device. **Testing** (macOS, no GPU; pyccel kernels compiled with GNU/Fortran; **GPU tests not run**): | | before (`devel`) | after | |---|---|---| | `fields_background/tests geometry/tests` | 92 passed, 12 skipped | 113 passed, 31 skipped | | `kinetic_background/tests post_processing/tests/test_orbits_tools.py feec/tests/test_field_init.py` | 50 passed | 50 passed | The 21 new passes are the 19 fake-CuPy cases, `test_all_equilibria_have_cases` and `test_host_call_numpy_is_plain_call`. The 19 new skips are the GPU tests (`test_equil_cupy`). The DESC case takes about 35 s, because it loads W7-X from the desc package. 🤖 Generated with [Claude Code](https://claude.com/claude-code) --------- Co-authored-by: Claude Opus 5.5 <noreply@anthropic.com>
max-models
added a commit
that referenced
this pull request
Oct 12, 2026
… generation) (#718) **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. **Decision: every CUDA kernel is written by hand. Kernels are not generated from the pyccel source** (no numba, `cupyx.jit` or other code generation). To decide, one guiding-center kernel, `bstar_parallel_3form`, was ported by hand and through a numba.cuda prototype generator. The prototype has been removed from this PR now that the decision is made. The short version of the reasons is in `CUDA_STRATEGY.md`, section *Hand-written CUDA (decision for #690)*: - **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, and 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. The prototype translated only 31 of the 60 pusher and accumulation entry kernels. - **It does not fit the existing infrastructure.** Structs, `CudaKernel` signature checks, the CPU emulation of the real `.cu` and the parity tests are all built for CUDA C++. numba also adds about 175 MB of dependencies. - **Drift between pyccel and CUDA is caught by the parity cases.** **Core changes:** - **Hand-written port of `bstar_parallel_3form`:** - `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. **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 (decision for #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 | 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](https://claude.com/claude-code) --------- Co-authored-by: Claude Opus 5.5 <noreply@anthropic.com>
max-models
added this pull request to stack #729
October 12, 2026 12:10
spossann
previously approved these changes
Oct 12, 2026
max-models
marked this pull request as ready for review
October 12, 2026 19:05
spossann
approved these changes
Oct 13, 2026
max-models
added a commit
that referenced
this pull request
Oct 13, 2026
**Stack:** part 10 of 14, based on #720 (merge that first), next: #724. Full order: #709 → #708 → #711 → #712 → #705 → #713 → #714 → #718 → #720 → #721 → #724 → #725 → #723 → #722. **Stack merge:** in `SplineFunction.__init__` (`feec/psydac_derham.py`), the `SplineArguments`/`CudaSplineArguments` are now built after #709's degree check and its host-to-device copy of the knots. --- **Solves the following issue(s):** Solves #678 (the remaining point: #667 (comment)). Part of #650. **Core changes:** `SplineFunction` no longer builds `(kind, pn, tn1, tn2, tn3, starts)` tuples for the spline evaluation kernels. It now uses an argument class pair, like `DerhamArguments`/`CudaDerhamArguments`: - New pyccel class `kernel_arguments/spline_args_kernels.SplineArguments(kind, pn, tn1, tn2, tn3, starts)` and its CUDA twin `kernel_arguments/spline_args_cuda.CudaSplineArguments` (a `CudaStructArguments` with struct `SplineArgs`). The header `kernel_arguments/spline_args.cuh` is generated by the new `write_spline_header()`. `SPLINE_STRUCTS` is added to `CUDA_STRUCTS`. - `SplineFunction.__init__` creates one argument object per component: `SplineArguments` on NumPy, `CudaSplineArguments` on CuPy. They are stored in `_args_spline`, which replaces `_args_eval`. The CUDA degree check (degrees 1 to 8) moves from `SplineFunction` into `CudaSplineArguments`, the same as in `CudaDerhamArguments`. - `eval_spline_mpi_markers`, `eval_spline_mpi_matrix` and `eval_spline_mpi_sparse_meshgrid` (pyccel and CUDA, still 1:1) take `args_spline` instead of the six loose arguments. The shared helper `eval_spline_mpi` (pyccel and `__device__`) keeps its flat signature, because other kernels call it too. - `CudaSplineArguments` accepts only C-contiguous knot vectors. `SplineFunction` already passes contiguous ones, so `test_evaluation_parity` now uses contiguous knots instead of strided views. No CUDA kernel other than these three takes the spline arguments, and pushers and accumulation already use `DerhamArguments`. One possible follow-up: `eval_spline_mpi_tensor_product_fixed` is a pyccel-only kernel in `evaluation_kernels_3d.py` that `SplineFunction` still calls with loose `kind, pn, starts` (no knots). Tests: - `test_kernel_backends.py`: the new pair is added to `ARGUMENT_PAIRS`, so the mirror check (and the struct layout check on a GPU) covers it. New test `test_generated_spline_header`. The owner test checks that `SplineFunction` creates `SplineArguments` on NumPy. - `kernel_test_args.spline_evaluation_arguments` returns `(_data, args_spline)`, so the parity and emulation cases now pass the argument object. `test_cuda_emulation.CUDA_CLASSES` includes `CudaSplineArguments`. Local tests (macOS, no GPU), after recompiling the changed kernels: - `bsplines/tests`, `feec/tests/test_eval_field.py`, `feec/tests/test_derham_gpu.py`, `pic/tests/test_cuda_emulation.py`, `pic/tests/test_kernel_backends.py`, `pic/tests/test_cuda_parity.py`, `pic/tests/test_pushers.py`, `pic/tests/test_accum_vec_H1.py`: 258 passed, 315 skipped. On devel, the same files gave 256 passed, 314 skipped; the differences are the 2 new tests and 1 new GPU-only skip. - `mpirun -n 2` on `bsplines/tests/test_eval_spline_mpi.py`, `feec/tests/test_eval_field.py`, `pic/tests/test_pushers.py`, `pic/tests/test_accum_vec_H1.py`: 139 passed (139 on devel). - CPU emulation of the three CUDA kernels against pyccel passes (`test_emulated_parity[eval_spline_mpi_*]`). - With cunumpy's fake CuPy, an ad-hoc check built `CudaSplineArguments` and packed the struct, and the degree check and the rejection of host arrays both work. The fake CuPy cannot construct a full `Derham`, so `SplineFunction` on CuPy was not run end-to-end. **GPU tests were not run** (no GPU available). This affects the struct layout, the parity tests and `test_spline_function_backends[cupy]`. **Model-specific changes:** None. **Documentation changes:** `CUDA_STRATEGY.md`: the new pair is added to the argument class table and the current-state table, plus a short note in the PR 16 implementation notes. 🤖 Generated with [Claude Code](https://claude.com/claude-code) --------- Co-authored-by: Claude Opus 5.5 <noreply@anthropic.com>
max-models
added a commit
that referenced
this pull request
Oct 13, 2026
**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; #705/#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 #716 (struphy-hub/feectools#98) is in. With #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](https://claude.com/claude-code) --------- Co-authored-by: Claude Opus 5.5 <noreply@anthropic.com>
max-models
added a commit
that referenced
this pull request
Oct 13, 2026
**Stack:** part 12 of 14, based on #724 (merge that first), next: #723. Full order: #709 → #708 → #711 → #712 → #705 → #713 → #714 → #718 → #720 → #721 → #724 → #725 → #723 → #722. --- **Solves the following issue(s):** Solves #411, tracked in #650. The FEEC assembly kernels are the last entry kernels that were still called as `PyccelKernel(getattr(module, name))` from shared modules. They now follow the folder rule of `CUDA_STRATEGY.md`, like the pusher, accumulation and geometry kernels. This is a pure restructuring: no numerical change. **Core changes:** - **One folder per kernel:** the 16 entry kernels of `feec/mass_kernels.py` and `feec/basis_projection_kernels.py` each move into `feec/kernels/<name>/`, which holds `<name>_kernels.py` and an `__init__.py` declaring `<name> = Kernel.from_folder(__name__, structs=CUDA_STRUCTS, missing_cuda="fallback")`. There is no `.cu` file, because none of these kernels has a CUDA port yet. Both old modules are removed. - From `mass_kernels.py`: `kernel_{1,2,3}d_mat`, `kernel_{1,2,3}d_vec`, `kernel_{1,2,3}d_eval`, `kernel_3d_matrixfree`, `kernel_3d_diag`, `surface_kernel_3d_mat`, `surface_kernel_3d_vec` - From `basis_projection_kernels.py`: `assemble_dofs_for_weighted_basisfuns_{1,2,3}d` - **The function bodies are unchanged.** I checked this with an AST comparison of all 16 functions against `devel`. Each file imports only what it uses (`numpy as np` and/or `shape`). - **Call sites** in `feec/mass.py`, `feec/boundary_mass.py` and `feec/basis_projection_ops.py` import the kernels. For kernels that depend on the dimension, they pick the right one from a small module-level dict (`_MAT_KERNELS[ldim]`, `_EVAL_KERNELS`, `_MATRIXFREE_KERNELS`, `_DIAG_KERNELS`, `_ASSEMBLY_KERNELS`) instead of `getattr` on a module. The direct pyccel calls of `mass_kernels.kernel_3d_vec` in `L2Projector` also go through the `Kernel` now. - **`missing_cuda="fallback"`:** on the CuPy backend these kernels keep the behaviour of the former `PyccelKernel` calls: they run on host copies of the arrays, with a warning the first time. With the default (`"raise"`), FEEC setups that currently assemble mass matrices on CuPy would start to fail. When a kernel is ported, its port removes the fallback for that kernel. - **`collect_kernel_files()`** picks up the new files with no change, since their names contain `kernels`. It now lists 142 files, 17 of them under `feec/kernels`. - **Tests:** `test_catalog_signatures` expects 17 kernels in `feec.kernels` (was 1). `feec.kernels` was already in `test_cuda_parity.PACKAGES`, so the new folders get the signature check and the folder-declaration check. **Model-specific changes:** None. **Documentation changes:** - `feec/kernels/__init__.py`: the docstring lists the assembly kernels and the fallback. - `CUDA_STRATEGY.md`: - a checklist entry for this PR - the kernel count in "Current state" - a short section "FEEC assembly kernels (#411)" saying that these kernels are now ready for CUDA ports (add `<name>_cuda.cu` and parity cases, then drop the fallback) **Testing** (macOS, no GPU; pyccel kernels compiled in the worktree with GNU/Fortran; **GPU tests not run**; following the agreed scope, I ran only a few directly relevant tests, and CI covers the rest): - `feec/tests/test_mass_reassembly.py`: 2 passed - `feec/tests/test_basis_ops.py::test_some_basis_ops` and `::test_transposed_update_weights_drops_zero_blocks`: passed - `feec/tests/test_boundary_integrals.py -k "unit_cube_constant or hcurl_per_face"`: 12 passed - `feec/tests/test_utilities_kernels.py`, `pic/tests/test_kernel_backends.py::test_catalog_signatures`, `pic/tests/test_cuda_parity.py::test_folders_declare_their_kernels`: passed - An ad-hoc script (not committed) on Colella with 6×5×4 elements and degree (2,3,2): - Matrix-free `M0`–`M3` `dot` agree with the assembled operators to about 1e-15 relative (`kernel_3d_matrixfree` against `kernel_3d_mat`). - Matrix-free and assembled `diagonal()` agree exactly (`kernel_3d_diag`). - `L2Projector.get_dofs` (`kernel_3d_vec`) and `eval_quad` (`kernel_3d_eval`) run. Not run locally: the full test suites, `test_mass_matrices.py` (32³ and 48³ grids, slow) and MPI runs. No MPI code changed. 🤖 Generated with [Claude Code](https://claude.com/claude-code) Co-authored-by: Claude Opus 5.5 <noreply@anthropic.com>
max-models
added a commit
that referenced
this pull request
Oct 13, 2026
**Stack:** part 14 of 14, based on #723 (merge that first), top of the stack. Full order: #709 → #708 → #711 → #712 → #705 → #713 → #714 → #718 → #720 → #721 → #724 → #725 → #723 → #722. **Stack fixes:** applied the `struphy_cuda::<pyccel module>` rule to the CUDA code of the PRs below: `linear_vlasov_ampere` (#705), `vlasov_maxwell` (#713), `bstar_parallel_3form` (#718), `eval_spline_mpi_tensor_product_fixed` (#724) and the `SplineArgs` versions of `eval_spline_mpi_*` (#721). Each of them has `using namespace struphy_cuda;` and qualified helper calls. The helpers `outer`, `fill_mat`, `fill_mat_vec` and `m_v_fill_b_v1_symm` (#705/#713) are inside their module namespaces, and `m_v_fill_b_v1_symm` calls `filler_kernels::`/`evaluation_kernels_3d::`. The `fill_v1_symm` test wrapper is qualified. The "open PRs that need the same treatment" below are therefore covered by this stack. Local re-run on the stacked branch: `test_cuda_parity.py`, `test_kernel_backends.py`, `test_evaluation_cuda.py`, `test_device_helpers.py` (59 passed, 419 skipped), `test_cuda_emulation.py` (74 passed; it compiles every `.cu`/`.cuh` with the host C++ compiler). --- **Solves the following issue(s):** Closes #672, tracked in #650. The CUDA device helpers now live in C++ namespaces that mirror the pyccel modules they port. A pyccel call `linalg_kernels.matrix_inv(...)` becomes `linalg_kernels::matrix_inv(...)` in CUDA, so the 1:1 correspondence is visible at every call site and two helpers with the same name in different modules (`det`/`det_df`, `f`/`df`, `push`/`pull`) can no longer clash. **Core changes:** - **One namespace per pyccel module**, nested in `struphy_cuda` as proposed in the issue (`struphy_cuda::linalg_kernels::matrix_inv`). The namespace is named after the module the header ports: - `bsplines_kernels`, `evaluation_kernels_2d`, `evaluation_kernels_3d` (`bsplines/`) - `evaluation_kernels`, `spline_mappings_kernels`, `transform_kernels` (`geometry/`) - one per mapping, named after its pyccel module: `cuboid_kernels`, `orthogonal_kernels`, `colella_kernels`, `hollow_cylinder_kernels`, `powered_elliptic_cylinder_kernels`, `hollow_torus_kernels`, `shafranov_shift_cylinder_kernels`, `shafranov_sqrt_cylinder_kernels`, `shafranov_dshaped_cylinder_kernels` (`geometry/domains/<name>/<name>_cuda.cuh`) - `linalg_kernels`, `filler_kernels`, `particle_to_mat_kernels`, `pusher_utilities_kernels` - **CUDA-only helpers and constants** go into the namespace of the module whose calculation they extract: `bsplines_kernels::MAX_SPLINE_DEGREE`, `evaluation_kernels_3d::SplineScratch`, `spline_mappings_kernels::spline_index_row`/`first_plane`, `transform_kernels::column_norm`, `pusher_utilities_kernels::reflect_velocity`. - **Call sites:** - Inside a header, calls into another module are qualified (`linalg_kernels::det(...)` in `evaluation_kernels::det_df`). Calls within the same module are not, as in pyccel. - Every `.cu` file writes `using namespace struphy_cuda;` after its includes and qualifies every helper call (`evaluation_kernels_3d::get_spans(...)`, `evaluation_kernels::df(...)`). - The CUDA sources embedded in `test_device_helpers.py` and `geometry/tests/spline_mapping_cases.py` use fully qualified names. - **Left global, on purpose:** - `struphy_cuda::pi` (`geometry/domains/constants_cuda.cuh`). It stands in for `numpy.pi`, not for a struphy module. - The generated argument structs `MarkerArgs`, `DerhamArgs`, `DomainArgs` and `LocalProjectorsArgs` (`kernel_arguments/*.cuh`, written by cunumpy's `write_cuda_header`). cunumpy matches the `extern "C" __global__` signatures against the struct names, and the CPU emulation (`pic/tests/cuda_emulation.py`) rebuilds the structs by name, so the generators and the emulation are unchanged. - **`evaluation_kernels_3d.cuh`:** it had three `namespace struphy_cuda` blocks with an `#include` in between. It now has one block. `eval_spline_mpi` was outside the namespace and is now `evaluation_kernels_3d::eval_spline_mpi`, like its pyccel counterpart. - Pure refactor: no function body, signature or entry kernel changes. Lines changed by the renaming were rewrapped at 120 columns. **Open PRs that need the same treatment after merge:** #705, #713 and #718 add `.cu`/`.cuh` code that still calls `struphy_cuda::<name>`. Whichever of them and this PR merges second needs the rule applied: `using namespace struphy_cuda;` in new `.cu` files, and new helpers in their module's namespace. Their branches were not touched. **Model-specific changes:** None. **Documentation changes:** `CUDA_STRATEGY.md`: a new convention bullet under "Target layout" documents the rule: the namespace is `struphy_cuda::<pyccel module>`, CUDA-only helpers go into their module, cross-module calls are qualified, `.cu` files use `using namespace struphy_cuda;`, and the exceptions are `pi` and the generated structs. **Local tests** (macOS, no GPU; pyccel kernels compiled with GNU/Fortran; **GPU tests not run**, so nothing was compiled with NVRTC). The CPU emulation compiles every `.cu` file and every helper header as C++17 with the host compiler. That is the syntax check for the namespaces; it caught one misplaced brace while I was working on this. - `src/struphy/pic/tests/test_cuda_parity.py` - `src/struphy/pic/tests/test_cuda_emulation.py` - `src/struphy/pic/tests/test_device_helpers.py` - `src/struphy/bsplines/tests/test_evaluation_cuda.py` - `src/struphy/pic/tests/test_kernel_backends.py` Result: 107 passed, 392 skipped (the GPU-only tests). `devel` gives the same counts. One open point for the GPU run: the nested namespace definition (`namespace struphy_cuda::bsplines_kernels {`) needs C++17. NVRTC uses C++17 by default since CUDA 11, and the emulation compiles with `-std=c++17`. 🤖 Generated with [Claude Code](https://claude.com/claude-code) Co-authored-by: Claude Opus 5.5 <noreply@anthropic.com>
This file contains hidden or bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment
Add this suggestion to a batch that can be applied as a single commit.This suggestion is invalid because no changes were made to the code.Suggestions cannot be applied while the pull request is closed.Suggestions cannot be applied while viewing a subset of changes.Only one suggestion per line can be applied in a batch.Add this suggestion to a batch that can be applied as a single commit.Applying suggestions on deleted lines is not supported.You must change the existing code in this line in order to create a valid suggestion.Outdated suggestions cannot be applied.This suggestion has been applied or marked resolved.Suggestions cannot be applied from pending reviews.Suggestions cannot be applied on multi-line comments.Suggestions cannot be applied while the pull request is queued to merge.Suggestion cannot be applied right now. Please check back later.
Stack: part 13 of 14, based on #725 (merge that first), next: #722. Full order: #709 → #708 → #711 → #712 → #705 → #713 → #714 → #718 → #720 → #721 → #724 → #725 → #723 → #722.
Stack fixes: added
OUTPUTSto the 16 FEEC assembly folders of #725, using the same arguments as the formerPyccelKernel(..., outputs=...)wraps inmass.py/boundary_mass.py/basis_projection_ops.py, pluskernel_*d_vecandsurface_kernel_3d_vec. Also addedOUTPUTSforeval_spline_mpi_tensor_product_fixed(#724). The indices ofeval_spline_mpi_{markers,matrix,sparse_meshgrid}were moved to theSplineArgumentssignatures of #721 (3, 5, 5).argument_arraysnow knowsCudaSplineArguments. Thetest_cuda_emulation.pymerge keeps the tests added below (test_emulated_v1_symm_filler). Local re-run on the stacked branch:test_cuda_parity.py,test_kernel_backends.py,test_evaluation_cuda.py(59 passed, 336 skipped),test_cuda_emulation.py(74 passed).Solves the following issue(s):
Solves #342, part of the CUDA work tracked in #650. Struphy's
PyccelKernels now declare which arguments they write to (outputs=), which cunumpy 0.6.1 already supports (PyccelKernel(..., outputs=...),Kernel.from_folder(..., host_options={"outputs": ...}), andassert_kernels_agreeuseskernel.host_kernel.outputsby default). Before this PR, no struphy kernel declared outputs. Without a declaration, the CuPy fallback copies every converted array back to the device after each call, and the parity tests compare every array.Core changes:
Kernel folders (all 90): each
<package>/kernels/<name>/__init__.pydeclaresThe declarations come from a static analysis of the pyccel sources, which follows writes through aliases (
markers = args_markers.markers) and through called kernel helpers. I reviewed the results by hand. The breakdown:args_markers(only itsmarkersfield).mat*/vec*arrays.markersorres.out/values/mat_f/f_eval*argument.The scratch arrays
bn1 … bd3ofDerhamArgumentsare written by some kernels but are not outputs. They are not in the CUDA struct either.Other
PyccelKernel(...)wraps passoutputs=too:WeightedMassOperatorassembly (kernel_{1,2,3}d_mat,-1),eval_quad(kernel_*d_eval,-1), the matrix-freedot(kernel_3d_matrixfree,-2, which writesout._data; the last argument isv._data) anddiagonal(kernel_3d_diag,-1), infeec/mass.pyBoundaryMassOperator(surface_kernel_3d_mat,-1), infeec/boundary_mass.pyBasisProjectionOperator(assemble_dofs_for_weighted_basisfuns_*d,0), infeec/basis_projection_ops.pyNegative indices keep one declaration valid for every dimension. The test-only wraps in
test_mat_vec_filler.pyandtest_kernel_setup.pyare unchanged.Parity tests use the declarations:
test_cuda_parity.py:KernelCatalog.from_packagebuilds its ownKernels and does not see the folders'OUTPUTS. The test therefore takes the kernels from the folder modules (DECLARED), and the GPUtest_parity(assert_kernels_agree) now compares only the declared outputs. A side effect, found by reading the code but not run on a GPU: before,assert_kernels_agreealso collected the hostDerhamArgumentsscratch arrays (argument i.bn1, …). The CUDA struct has no such fields, so the host and CUDA key sets could not match for kernels takingargs_derham.test_kernels_declare_outputs(no GPU needed): every kernel folder declares non-empty, in-rangeOUTPUTS.test_declared_outputs(NumPy backend, every parity case): the pyccel kernel changes at least one declared output and no other argument array. Argument objects are compared through their CUDA struct fields, so scratch is ignored. The helperargument_arrayslives incuda_parity_cases.py.test_cuda_emulation.test_emulated_paritynow compares only the outputs with the pyccel kernel, and also asserts that the emulated CUDA kernel leaves every non-output array unchanged. I checked both new checks by mutation: declaringpush_eta_stagewithOUTPUTS = (4,)makes both tests fail.CUDA_STRATEGY.md: short note onOUTPUTS.Model-specific changes:
None. Behaviour on the NumPy backend is unchanged, because
PyccelKernelusesoutputsonly when it converts CuPy arrays.Documentation changes:
A short note in
CUDA_STRATEGY.md.Local tests run (macOS, no GPU; GitHub CI covers the rest):
pic/tests/test_cuda_parity.py: 18 passed, 213 skipped (the GPU tests)pic/tests/test_cuda_emulation.py -k emulated_parity: all 15 kernels passed. This ran as two invocations: the first hit my 5-minute limit after 13 kernels, and the remaining two (kernel_pullpush_pic,kernel_pullpush) passed separately in 80 s.pic/tests/test_kernel_backends.pyandtest_kernel_setup.py: 89 passed, 22 skippedfeec/tests/test_boundary_integrals.pyandtest_mass_reassembly.py: 25 passedGPU tests were not run. On the GPU,
test_paritynow compares only the declared outputs. Theoutputs=of the FEEC mass/projection wraps only take effect on the CuPy fallback path. No local test exercises them, neither on a GPU nor with the fake CuPy.🤖 Generated with Claude Code