Repository navigation
feectools: inner products stay on the device - #716
Merged
Merged
Conversation
Moves the feectools submodule to struphy-hub/feectools#96 (branch device-kronecker-solve): KroneckerLinearSolver solves device data with dense inverses of the 1D matrices (one GEMM per direction) instead of copying it to the host for LAPACK/SuperLU, so MassMatrixPreconditioner solves on the CuPy backend make no host/device transfers. Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
Point the feectools submodule at struphy-hub/feectools#97 (branch stencil-6d-views): the 3D stencil dot/transpose kernels take the matrix data as 6D views, and StencilMatrix.dot/vdot/transpose no longer raise for matrices with fewer than 2p+1 diagonals. Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
Move the feectools submodule to struphy-hub/feectools#98: on the CuPy backend inner products return 0-d device arrays and the Krylov solvers copy one scalar per iteration to the host (the convergence test). Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
This was referenced Oct 7, 2026
…6d-views Stack #703 -> #706: feectools submodule points to the head of struphy-hub/feectools#97, which now contains #96. Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
This was referenced Oct 7, 2026
Merged
max-models
marked this pull request as draft
October 8, 2026 09:13
max-models
marked this pull request as ready for review
October 8, 2026 09:14
spossann
previously approved these changes
Oct 8, 2026
spossann
approved these changes
Oct 9, 2026
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
#734) **Solves the following issue(s):** Part of #689, tracked in #650. This PR covers the items of #689 that are left for `VlasovAmpereOneSpecies`: running on 1 and 2 ranks, no host/device transfers in the time loop, and the diagnostics it needs. It does not close #689, because these runs use the CPU emulation; the first run on a real GPU is #687. **Stack:** this is **gpu-1** of a new stack. Merge it after #722 (`cuda-namespaces`), #716, #717, #719 and #765 (MPI sort of the markers after drawing, needed for the 2-rank runs). It merges those branches in, so until they are merged, this diff also shows their commits. **gpu-2** (the benchmark script) is based on this branch. **Core changes:** - **New test `models/tests/test_gpu_VlasovAmpereOneSpecies.py`.** A small Landau-damping run on NumPy and on CuPy: 8 cells, 320 Sobol markers, 3 steps, full-f weights, and an `e1_v1` binning plot so that the diagnostics run on the device too. It compares the energies, the `e_field` and `phi` coefficients, and the markers by ID (tolerances explained in `compare_runs`), and checks the transfer budget of the time loop. - **Shared helpers in `models/tests/gpu_e2e.py`.** The transfer budget, the instrumentation, gathering markers across ranks and the device context, moved out of #719's test. `test_gpu_LinearVlasovAmpereOneSpecies.py` now uses them. - **The GPU tests also run without a GPU.** On cunumpy's fake CuPy (`CUNUMPY_FAKE_CUPY=1`), every CUDA launch is emulated on the CPU and `CudaKernel.compile` is skipped. `Simulation.compile_cuda_kernels` still checks that every kernel of the time loop has a CUDA version. Without a GPU and without the fake CuPy, `test_fake_cupy_matches_numpy` runs the comparison in a serial child process. - **Transfer budget of `model.integrate`.** Only scalar downloads are allowed: the residual norm of each Krylov iteration, which feectools copies for its convergence test (feectools#98 keeps everything else of the solve on the device). No arrays, no host-to-device copies, no host kernels. - **Bug fix (every backend): `VlasovAmpereCoupling` updates the ghost regions of the mid-point field before the push.** The pusher kernel reads `_e_sum._data` including its ghost regions, and `_e_sum` is built with `+=` from the solver output, which does not sync them. Stale ghost regions would give a wrong field near process (and periodic) boundaries. - **`Particles`.** Checking whether MPI is CUDA-aware now happens when the marker array is allocated (in the setup), not on the first marker exchange. Before, that check was the only array transfer inside the time loop on 2 ranks. **Results** (fake CuPy, every CUDA launch emulated on the CPU; macOS; **not run on a GPU**): | test | 1 rank | 2 ranks | |---|---|---| | `test_gpu_VlasovAmpereOneSpecies.py` (NumPy reference + CuPy vs NumPy + transfer budget) | passed | passed | | `test_gpu_LinearVlasovAmpereOneSpecies.py::test_cupy_matches_numpy_without_transfers` | passed | not run | Transfers in the time loop of the CuPy run: - `integrate`: one 8-byte residual norm per pcg iteration. - Diagnostics: one 24-byte copy of the three scalars per update. - Output: one copy per saved dataset. **Diagnostics:** no new CUDA kernels are needed for this model. The scalars use `xp` reductions and feectools' device inner products, and binning uses `xp.histogramdd`. **Found along the way (not changed here):** now #766 (Poisson stabilization) and #767 (energy diagnostic with the control variate). - The initial Poisson solve with the default `stab_eps=1e-14` puts (Monte-Carlo net charge) / 1e-14 into the constant mode of `phi`. That makes `phi` about 1e14, so E = -grad(phi) loses about 5% to cancellation on either backend. The verification tests use this default. The new test uses `stab_eps=1e-6`, as #719 does. - With `control_variate=True`, `kinetic_energy` at t=0 is computed with the full-f weights and afterwards with the delta-f weights, so `total_energy` jumps after the first step. With full-f weights the total is conserved from step to step to 1e-9, but the t=0 value still differs by about 2e-4 (not investigated). **Review fixes** (pushed by @spossann, see the review): - the e2e tests run in CI: unit shard-2 (serial and 2 ranks) and on the fake CuPy on 1 and 2 ranks (`test-PR-unit-mpi.yml`); merged #765, without which they segfault on 2 ranks (#764); - `MassMatrixPreconditioner._core_diagonal` returns a diagonal operator again (with the wrapped band entries); - `test_preconditioner_cupy.py` uses cunumpy's `emulated_launches` and `run_fake_cupy_child`, and its GPU test emulates its launches on the fake CuPy; - `FFTSolver.solve` asserts host input. **Documentation changes:** `CUDA_STRATEGY.md`: the end-to-end entry now covers both models, the CPU emulation and the transfer budget. 🤖 Generated with [Claude Code](https://claude.com/claude-code) --------- Co-authored-by: Claude Opus 5.5 <noreply@anthropic.com> Co-authored-by: Stefan Possanner <stefan.possanner@ipp.mpg.de>
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
2 participants
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 3 of 3 — based on #706 (
feectools-stencil-6d-views, merge that first, after #703), next: none (top of the stack).Corresponding update in feectools: struphy-hub/feectools#98. Mirrors the feectools stack struphy-hub/feectools#96 → struphy-hub/feectools#97 → struphy-hub/feectools#98.
feectools-stencil-6d-viewsis merged into this branch (no rebase), and thefeectoolssubmodule now points at the stacked head of struphy-hub/feectools#98,6d88806, which contains feectools#96 and #97. Against its base (#706) this PR only moves the submoduled16a6ad→6d88806. After each feectools PR of the stack is merged intodevel-tiny, thefeectoolssubmodule of the struphy stack must be moved to the corresponding mergeddevel-tinycommit (thepr-feectools-submodulecheck fails until then).Solves the following issue(s):
Part of #689 (
VlasovAmpereOneSpeciesend to end on the GPU, no host/device transfers in the time loop), tracked in #650. This PR moves thefeectoolssubmodule to struphy-hub/feectools#98 ("Inner products stay on the device"). On the CuPy backend, every CG iteration of a struphy solve copied two scalars from the device to the host; after this change it copies one (the convergence test).The
pr-feectools-submodulecheck fails until struphy-hub/feectools#98 is merged intodevel-tiny. After that, the submodule should point at the merge commit.Core changes:
feectoolssubmodule:d16a6ad(base feectools: stencil 3D device kernels with 6D views #706) →6d88806(branchinner-on-deviceof Inner products stay on the device feectools#98, stacked on feectools#96 and Remove deprecated console commands #97; originallya15a8e8). No struphy code changes.What changes in feectools, and only on the CuPy backend:
StencilVector.inner/BlockVector.inner/dot_innerreturn a 0-d device array, in serial and with MPI.axpyaccepts such a scalar without copying it to the host.On NumPy, results and types are unchanged.
Copies to the host per iteration, counted under cunumpy's fake CuPy:
struphy call sites of
.innerthat run on CuPy now get a 0-d device array. I checked them:models/scalars.pywrite it intoxpbuffers (local_value[0] = ...), which works on the device;linear_mhd.py,shear_alfven.py, ...) return it, and the scalar machinery handles it like the existingdot_innerresults;PolarVector.dotadds a NumPy scalar to it, which works, but polar splines are not on CuPy yet (Polar splines on the CuPy backend #695);.inner. It works with device scalars: the arithmetic stays on the device, and the comparisons with 0 copy to the host implicitly.This removes the first blocker listed in CUDA version of the linear_vlasov_ampere accumulation #705 ("CG inner products ... The Schur solve is outside the transfer guard because of this"). The remaining copy per CG iteration still counts as a transfer for
assert_no_transfers. The transfer guard could then cover the Schur solve with an allowance of one 8-byte copy per iteration.Model-specific changes:
None.
Documentation changes:
None in struphy. The feectools PR documents the change in feectools'
CUDA_STRATEGY.md(section "Inner products on the device").Testing (macOS, no GPU; struphy kernels compiled with GNU/Fortran; GPU tests not run):
feectools (see Inner products stay on the device feectools#98):
ddm/tests/test_cart_*d.py;mpirun -n 2linalg MPI suite: 738 → 739 passed;struphy, solver tests on the NumPy backend, with feectools
devel-tiny(before) and this branch (after):propagators/tests/test_poisson.py -m "not mpi" -k "not multigrid": 40 passed before, 40 passed after;linear_algebra/tests/test_saddlepoint_massmatrices.py: the first case (SaddlePointSolverUzawaNumpy) passed before and after.I stopped the remaining saddle-point and multigrid tests: with the feectools kernels uncompiled in the checkouts, each test took more than 10 minutes. On NumPy, feectools' change is limited to identity helpers and an equivalent loop structure in BiCGStab, and the feectools solver tests check the same iteration counts and results.
🤖 Generated with Claude Code