Repository navigation
GPU stack 6/8: marker rows prepared inside the push kernels - #787
Merged
max-models merged 12 commits intoOct 24, 2026
Merged
Conversation
…the push kernels - Particles allocates the marker buffer column-major on the CuPy backend: kernels and the Python column operations access a few columns of many markers, so with a row-major buffer every access touched whole rows and the threads of a warp read addresses a row apart. The CUDA kernels address the buffer through its strides (Array2D<double> in MarkerArgs); NumPy keeps C order for the pyccel kernels. - push_eta_stage and push_vxb_analytic prepare each marker row themselves at stage 0 (MarkerArgs.prepare_rows), so Pusher._push skips the separate prepare_push pass when no init or evaluation kernel runs before them. prepare_push is launched on a 2D grid. Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
…repare_rows in the CUDA emulation Kernels without a CUDA version run on host copies of the device arrays, which keep the F order (cunumpy.to_numpy uses get(order='A')), and the pyccel kernels reject it: the sorting-box kernels failed on the fake CuPy (test_sorting_boxes_on_cupy_fake_cupy), and the same holds on a GPU for every marker kernel without a CUDA version. The buffer is therefore C-ordered unless the variable is set. The CUDA emulation of the tests reads the struct fields from the pyccel argument objects, which do not have the CUDA-only prepare_rows flag (every emulated marker kernel failed). Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
This was referenced Oct 16, 2026
Merged
max-models
added this pull request to stack #790
October 16, 2026 08:51
…k-5-diffusion-solve) # Conflicts: # src/struphy/kernel_arguments/pusher_args_cuda.py # src/struphy/pic/base.py # src/struphy/pic/pushing/prepare_push_cuda.cu # src/struphy/pic/pushing/pusher.py
…yout (feectools update)
max-models
changed the title
…yout (fused gradient)
…yout (feectools update, CachedLaunch)
…k-5-diffusion-solve)
…yout (pcg_sr in the shear Alfven solves)
max-models
marked this pull request as draft
October 22, 2026 06:39
max-models
marked this pull request as ready for review
October 23, 2026 07:15
spossann
previously approved these changes
Oct 24, 2026
spossann
approved these changes
Oct 24, 2026
max-models
added a commit
that referenced
this pull request
Oct 24, 2026
…ion PCG); pcg_sr solver option (#782) Part 1 of 8 of the split of #772 (`cuda-aware-mpi-fixes`) into stacked PRs. **Requires struphy-hub/feectools#104 to be merged into feectools `devel-tiny` first** (PR 1 bumps the submodule). **Merge struphy-hub/feectools#104 (feectools `fix-cuda-aware-mpi-issues` into `devel-tiny`) first.** This PR points the `feectools` submodule at that branch (commit `1d8d47c`). - Submodule update: device ghost-exchange plans, fused block exchanges, `PConjugateGradientSingleReduction`, faster `stencil_dot_3d` on the GPU (details in the feectools PR). - `pcg_sr` is accepted wherever `pcg` is (`LiteralOptions.OptsSymmSolver`), also with a preconditioner: Poisson, ImplicitDiffusion, GyrokineticPoissonSolve, CurlCurlSolve, VariationalDensityEvolve, VariationalMomentumAdvection. Same iterates as `pcg`, one global reduction per iteration instead of two. ## Tests Local, macOS without a GPU, in a fresh environment with feectools `1d8d47c` and the pinned cunumpy, after merging devel (with #774 and #775). The PR passes `test_cuda_parity`, `test_cuda_emulation`, `test_kernel_backends`, `test_sorting`, `test_accum_matrix_cupy`, `test_preconditioner_cupy`, `test_poisson` and `test_output` on its own. The top of the stack also passes the mass-matrix, preconditioner, Derham-GPU, fused-gradient, Kronecker-block, Poisson, shear Alfvén, post-processing, polar fake-CuPy and all `pic` tests (1477 passed), the examples test for the two profiling cases, the sorting, Poisson and fused-gradient tests under `mpirun -n 2`, and the device MPI sort on the fake CuPy on 2 and 4 ranks. The Kronecker products match the stencil mass matrices to 5e-16 for M0, M1, M2, M3, Mv, M1n and M2n on a Cuboid, periodic and Dirichlet, on 1 and 2 ranks. **Not tested on a GPU here.** The CUDA kernels and the device paths need a run on the cluster. ## Stack (merge in order) 0. struphy-hub/feectools#104 (feectools, into `devel-tiny`) | # | PR | base | |---|---|---| | 1 | #782 update feectools (device ghost exchange, single-reduction PCG); pcg_sr solver option **(this PR)** | `devel` | | 2 | #783 host/device conversions in post-processing; math.sqrt in Species | #782 | | 3 | #784 CUDA mass-matrix assembly kernel kernel_3d_mat | #783 | | 4 | #785 apply separable mass matrices as Kronecker products of 1D banded matrices | #784 | | 5 | #786 Poisson/diffusion solve with fewer passes over the vectors | #785 | | 6 | #787 marker rows prepared inside the push kernels | #786 | | 7 | #788 chunked marker loading, row-based weight update, leaner MPI sorting; static field in PushVxB | #787 | | 8 | #789 profiling cases Poisson cube and Vlasov slab (CPU node vs GPU, GPU scaling) | #788 | 🤖 Generated with [Claude Code](https://claude.com/claude-code) --------- Co-authored-by: Claude Opus 5.5 <noreply@anthropic.com> Co-authored-by: Copilot <223556219+Copilot@users.noreply.github.com>
max-models
added a commit
that referenced
this pull request
Oct 24, 2026
…in Species (#783) Part 2 of 8 of the split of #772 (`cuda-aware-mpi-fixes`) into stacked PRs. **Requires struphy-hub/feectools#104 to be merged into feectools `devel-tiny` first** (PR 1 bumps the submodule). - `post_processing/output.py`: on the CuPy backend the fields are evaluated on the device and written from the host. The evaluation grids, point data and coordinates are converted where they cross, and HDF5 data is moved to the device before it goes into the spline coefficients. - `Species`: the plasma frequency of the unit conversion is a host scalar (`math.sqrt` instead of `xp.sqrt`). ## Tests Local, macOS without a GPU, in a fresh environment with feectools `1d8d47c` and the pinned cunumpy, after merging devel (with #774 and #775). The PR passes `test_cuda_parity`, `test_cuda_emulation`, `test_kernel_backends`, `test_sorting`, `test_accum_matrix_cupy`, `test_preconditioner_cupy`, `test_poisson` and `test_output` on its own. The top of the stack also passes the mass-matrix, preconditioner, Derham-GPU, fused-gradient, Kronecker-block, Poisson, shear Alfvén, post-processing, polar fake-CuPy and all `pic` tests (1477 passed), the examples test for the two profiling cases, the sorting, Poisson and fused-gradient tests under `mpirun -n 2`, and the device MPI sort on the fake CuPy on 2 and 4 ranks. The Kronecker products match the stencil mass matrices to 5e-16 for M0, M1, M2, M3, Mv, M1n and M2n on a Cuboid, periodic and Dirichlet, on 1 and 2 ranks. **Not tested on a GPU here.** The CUDA kernels and the device paths need a run on the cluster. ## Stack (merge in order) 0. struphy-hub/feectools#104 (feectools, into `devel-tiny`) | # | PR | base | |---|---|---| | 1 | #782 update feectools (device ghost exchange, single-reduction PCG); pcg_sr solver option | `devel` | | 2 | #783 host/device conversions in post-processing; math.sqrt in Species **(this PR)** | #782 | | 3 | #784 CUDA mass-matrix assembly kernel kernel_3d_mat | #783 | | 4 | #785 apply separable mass matrices as Kronecker products of 1D banded matrices | #784 | | 5 | #786 Poisson/diffusion solve with fewer passes over the vectors | #785 | | 6 | #787 marker rows prepared inside the push kernels | #786 | | 7 | #788 chunked marker loading, row-based weight update, leaner MPI sorting; static field in PushVxB | #787 | | 8 | #789 profiling cases Poisson cube and Vlasov slab (CPU node vs GPU, GPU scaling) | #788 | 🤖 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 24, 2026
Part 3 of 8 of the split of #772 (`cuda-aware-mpi-fixes`) into stacked PRs. **Requires struphy-hub/feectools#104 to be merged into feectools `devel-tiny` first** (PR 1 bumps the submodule). - `kernel_3d_mat_cuda.cu`: one thread per (element, codomain basis function, domain basis function). It replaces the host fallback on the CuPy backend (`missing_cuda="fallback"` removed, launch size from `n_threads_from`). - `WeightedMassOperator.assemble` reads the quadrature grids (spans, weights, basis values) as backend arrays, copied to the device once and cached. - `Simulation.allocate` compiles the assembly kernel before the first assembly, in its own profiling region (`setup: compile cuda kernels (assembly)`), so the ~1 s compilation is not counted in the assembly. Added while splitting: parity cases for `kernel_3d_mat` (random quadrature data on 3×2×2 elements with mixed degrees, with and without a ghost row). `test_cuda_kernels_have_parity_cases` requires one for every kernel with a CUDA version. The emulation test confirms that the CUDA kernel matches the pyccel kernel without a GPU. ## Tests Local, macOS without a GPU, in a fresh environment with feectools `1d8d47c` and the pinned cunumpy, after merging devel (with #774 and #775). The PR passes `test_cuda_parity`, `test_cuda_emulation`, `test_kernel_backends`, `test_sorting`, `test_accum_matrix_cupy`, `test_preconditioner_cupy`, `test_poisson` and `test_output` on its own. The top of the stack also passes the mass-matrix, preconditioner, Derham-GPU, fused-gradient, Kronecker-block, Poisson, shear Alfvén, post-processing, polar fake-CuPy and all `pic` tests (1477 passed), the examples test for the two profiling cases, the sorting, Poisson and fused-gradient tests under `mpirun -n 2`, and the device MPI sort on the fake CuPy on 2 and 4 ranks. The Kronecker products match the stencil mass matrices to 5e-16 for M0, M1, M2, M3, Mv, M1n and M2n on a Cuboid, periodic and Dirichlet, on 1 and 2 ranks. **Not tested on a GPU here.** The CUDA kernels and the device paths need a run on the cluster. ## Stack (merge in order) 0. struphy-hub/feectools#104 (feectools, into `devel-tiny`) | # | PR | base | |---|---|---| | 1 | #782 update feectools (device ghost exchange, single-reduction PCG); pcg_sr solver option | `devel` | | 2 | #783 host/device conversions in post-processing; math.sqrt in Species | #782 | | 3 | #784 CUDA mass-matrix assembly kernel kernel_3d_mat **(this PR)** | #783 | | 4 | #785 apply separable mass matrices as Kronecker products of 1D banded matrices | #784 | | 5 | #786 Poisson/diffusion solve with fewer passes over the vectors | #785 | | 6 | #787 marker rows prepared inside the push kernels | #786 | | 7 | #788 chunked marker loading, row-based weight update, leaner MPI sorting; static field in PushVxB | #787 | | 8 | #789 profiling cases Poisson cube and Vlasov slab (CPU node vs GPU, GPU scaling) | #788 | 🤖 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 24, 2026
… 1D banded matrices (#785) Part 4 of 8 of the split of #772 (`cuda-aware-mpi-fixes`) into stacked PRs. **Requires struphy-hub/feectools#104 to be merged into feectools `devel-tiny` first** (PR 1 bumps the submodule). - `feec/kronecker_mass.py`: `WeightedMassOperator.dot` applies the blocks of its matrix that are Kronecker products of 1D banded matrices (constant weights, e.g. on a Cuboid) as three banded passes instead of the 3D stencil product. On the GPU this uses the new kernel `kron_band_pass_cuda.cu`. `STRUPHY_KRONECKER_MASS=0` switches it off. - The replacement is built on first use and invalidated by assembly, copies and in-place arithmetic. A matrix changed by arithmetic that is not a Kronecker product (e.g. a Schur complement rebuilt every step) is not checked again. - The boundary and extraction operators around the matrix are skipped when they are identities, instead of being applied as copies. - The band-pass kernel is launched through feectools' `CachedLaunch` (feectools#104): the banded passes run with the same buffers in every solver iteration, so the launch preparation is done once per set of buffers. ## Tests Local, macOS without a GPU, in a fresh environment with feectools `1d8d47c` and the pinned cunumpy, after merging devel (with #774 and #775). The PR passes `test_cuda_parity`, `test_cuda_emulation`, `test_kernel_backends`, `test_sorting`, `test_accum_matrix_cupy`, `test_preconditioner_cupy`, `test_poisson` and `test_output` on its own. The top of the stack also passes the mass-matrix, preconditioner, Derham-GPU, fused-gradient, Kronecker-block, Poisson, shear Alfvén, post-processing, polar fake-CuPy and all `pic` tests (1477 passed), the examples test for the two profiling cases, the sorting, Poisson and fused-gradient tests under `mpirun -n 2`, and the device MPI sort on the fake CuPy on 2 and 4 ranks. The Kronecker products match the stencil mass matrices to 5e-16 for M0, M1, M2, M3, Mv, M1n and M2n on a Cuboid, periodic and Dirichlet, on 1 and 2 ranks. **Not tested on a GPU here.** The CUDA kernels and the device paths need a run on the cluster. ## Stack (merge in order) 0. struphy-hub/feectools#104 (feectools, into `devel-tiny`) | # | PR | base | |---|---|---| | 1 | #782 update feectools (device ghost exchange, single-reduction PCG); pcg_sr solver option | `devel` | | 2 | #783 host/device conversions in post-processing; math.sqrt in Species | #782 | | 3 | #784 CUDA mass-matrix assembly kernel kernel_3d_mat | #783 | | 4 | #785 apply separable mass matrices as Kronecker products of 1D banded matrices **(this PR)** | #784 | | 5 | #786 Poisson/diffusion solve with fewer passes over the vectors | #785 | | 6 | #787 marker rows prepared inside the push kernels | #786 | | 7 | #788 chunked marker loading, row-based weight update, leaner MPI sorting; static field in PushVxB | #787 | | 8 | #789 profiling cases Poisson cube and Vlasov slab (CPU node vs GPU, GPU scaling) | #788 | 🤖 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 24, 2026
…tors (#786) Part 5 of 8 of the split of #772 (`cuda-aware-mpi-fixes`) into stacked PRs. **Requires struphy-hub/feectools#104 to be merged into feectools `devel-tiny` first** (PR 1 bumps the submodule). - `BoundaryOperator.dot` on device Stencil/Block vectors computes `out = v * mask`, one kernel per block (one for all blocks when they are contiguous float64). The masks are found once from `apply_essential_bc_to_array`, so they zero exactly the same entries. This replaces a copy and one fill per boundary and component. - `WeightedMassOperator.without_boundary_ops()` / `MassWithoutBoundaryOps`: the mass operator without its boundary operators. In `grad.T @ M1 @ grad`, `grad` already applies the same 0/1 masks, so this is the same operator with two fewer passes per product. It is not used with the multigrid preconditioner. - `ImplicitDiffusion`: the automatic `1e-14` stabilization stays out of the solver operator when a Dirichlet boundary already makes `grad.T M1 grad` positive definite. There is no `stab_mat` product for the right-hand side when `sigma_2 = 0` (Poisson). - `struphy.feec.fused_derivatives`: `fuse_boundary_gradient(derham.grad)` and `fuse_boundary_gradient(derham.grad.T)` apply `B_out G B_in` with one CUDA kernel each on the GPU (`fused_derivatives_cuda.cu`), instead of a masked copy, three derivatives and a masked copy. The kernel does the same operations in the same order, so the products are bit for bit those of the composed operators. It reads the input with its ghost regions exchanged, times the input mask with exchanged ghost regions, which equals the exchanged masked input in serial and in parallel. Host data and the fake CuPy use the composed operators. `ImplicitDiffusion` uses it for `grad` and `grad.T` (not with the multigrid preconditioner). The original branch imported this module, but it was never committed, so it is written new here. `test_fused_derivatives` compares the kernels with the composed operators on the fake CuPy with emulated launches, and guards that the kernels actually run. Also added: the fused mask kernel is skipped on cunumpy's fake CuPy, which cannot launch raw kernels (`test_mass_preconditioners_fake_cupy` failed). ## Tests Local, macOS without a GPU, in a fresh environment with feectools `1d8d47c` and the pinned cunumpy, after merging devel (with #774 and #775). The PR passes `test_cuda_parity`, `test_cuda_emulation`, `test_kernel_backends`, `test_sorting`, `test_accum_matrix_cupy`, `test_preconditioner_cupy`, `test_poisson` and `test_output` on its own. The top of the stack also passes the mass-matrix, preconditioner, Derham-GPU, fused-gradient, Kronecker-block, Poisson, shear Alfvén, post-processing, polar fake-CuPy and all `pic` tests (1477 passed), the examples test for the two profiling cases, the sorting, Poisson and fused-gradient tests under `mpirun -n 2`, and the device MPI sort on the fake CuPy on 2 and 4 ranks. The Kronecker products match the stencil mass matrices to 5e-16 for M0, M1, M2, M3, Mv, M1n and M2n on a Cuboid, periodic and Dirichlet, on 1 and 2 ranks. **Not tested on a GPU here.** The CUDA kernels and the device paths need a run on the cluster. ## Stack (merge in order) 0. struphy-hub/feectools#104 (feectools, into `devel-tiny`) | # | PR | base | |---|---|---| | 1 | #782 update feectools (device ghost exchange, single-reduction PCG); pcg_sr solver option | `devel` | | 2 | #783 host/device conversions in post-processing; math.sqrt in Species | #782 | | 3 | #784 CUDA mass-matrix assembly kernel kernel_3d_mat | #783 | | 4 | #785 apply separable mass matrices as Kronecker products of 1D banded matrices | #784 | | 5 | #786 Poisson/diffusion solve with fewer passes over the vectors **(this PR)** | #785 | | 6 | #787 marker rows prepared inside the push kernels | #786 | | 7 | #788 chunked marker loading, row-based weight update, leaner MPI sorting; static field in PushVxB | #787 | | 8 | #789 profiling cases Poisson cube and Vlasov slab (CPU node vs GPU, GPU scaling) | #788 | 🤖 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.
Part 6 of 8 of the split of #772 (
cuda-aware-mpi-fixes) into stacked PRs. Requires struphy-hub/feectools#104 to be merged into feectoolsdevel-tinyfirst (PR 1 bumps the submodule).push_eta_stageandpush_vxb_analyticprepare each marker themselves at stage 0 (MarkerArgs.prepare_rows): they copy the phase-space coordinates to the rows fromfirst_pusher_idxand zero the rows after them, asprepare_pushdoes.Pusher._pushthen skips the separateprepare_pushlaunch when no init or evaluation kernel runs before them (those may read the prepared rows). The kernels cover the same markers, one thread per used marker (n_markers_used).prepare_rowsis a CUDA-only field ofMarkerArgs. The CUDA emulation of the tests reads it with default 0.Updated after #774 and #775 were merged into devel: on the original branch this PR also stored the marker buffer column-major on the GPU (
order="F"). devel now stores it as(n_cols, n_markers)on both backends, so that part is gone. The row preparation is ported tomarkers(j, ip)indexing.Tests
Local, macOS without a GPU, in a fresh environment with feectools
1d8d47cand the pinned cunumpy, after merging devel (with #774 and #775). The PR passestest_cuda_parity,test_cuda_emulation,test_kernel_backends,test_sorting,test_accum_matrix_cupy,test_preconditioner_cupy,test_poissonandtest_outputon its own. The top of the stack also passes the mass-matrix, preconditioner, Derham-GPU, fused-gradient, Kronecker-block, Poisson, shear Alfvén, post-processing, polar fake-CuPy and allpictests (1477 passed), the examples test for the two profiling cases, the sorting, Poisson and fused-gradient tests undermpirun -n 2, and the device MPI sort on the fake CuPy on 2 and 4 ranks. The Kronecker products match the stencil mass matrices to 5e-16 for M0, M1, M2, M3, Mv, M1n and M2n on a Cuboid, periodic and Dirichlet, on 1 and 2 ranks.Not tested on a GPU here. The CUDA kernels and the device paths need a run on the cluster.
Stack (merge in order)
devel-tiny)devel🤖 Generated with Claude Code