Skip to content

GPU stack 3/8: CUDA mass-matrix assembly kernel kernel_3d_mat - #784

Merged
max-models merged 2 commits into
gpu-stack-2-cupy-fixesfrom
gpu-stack-3-mass-assembly
Oct 24, 2026
Merged

max-models merged 2 commits into
gpu-stack-2-cupy-fixesfrom
gpu-stack-3-mass-assembly

Conversation

@max-models

@max-models max-models commented Oct 16, 2026 •

Copy link
Copy Markdown
Member

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)

  1. Device ghost-exchange plans, fused block exchanges, single-reduction PCG (CUDA-aware MPI) 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

spossann
spossann previously approved these changes Oct 24, 2026
max-models and others added 2 commits October 24, 2026 13:40
- kernel_3d_mat_cuda.cu: one thread per (element, codomain basis function, domain basis function);
  replaces the host fallback on the CuPy backend.
- 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 it is first called, so the compilation
  is not counted in the assembly's profiling regions.

Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
Random quadrature data on 3 x 2 x 2 elements with mixed degrees, with the first row local or a
ghost row. test_cuda_kernels_have_parity_cases requires a case for every kernel with a CUDA
version; the emulation test compares it with the pyccel kernel without a GPU.

Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
@max-models
max-models merged commit 9912a21 into devel Oct 24, 2026
28 checks passed
@max-models
max-models deleted the gpu-stack-3-mass-assembly branch October 24, 2026 19:03
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
… 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>
max-models added a commit that referenced this pull request Oct 24, 2026
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 feectools `devel-tiny` first** (PR 1 bumps the submodule).

- `push_eta_stage` and `push_vxb_analytic` prepare each marker
themselves at stage 0 (`MarkerArgs.prepare_rows`): they copy the
phase-space coordinates to the rows from `first_pusher_idx` and zero the
rows after them, as `prepare_push` does. `Pusher._push` then skips the
separate `prepare_push` launch 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_rows` is a CUDA-only field of `MarkerArgs`. 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 to
`markers(j, ip)` indexing.

## 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 |
#785 |
| 6 | #787 marker rows prepared inside the push kernels **(this PR)** |
#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>
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment

Labels

None yet

Projects

None yet

Development

Successfully merging this pull request may close these issues.

2 participants