Cuda 1 xp arrays - #85
Merged
Merged
Conversation
Make feectools work when cunumpy's backend is CuPy, without any device kernels yet: kernels that stay on the host still copy their arrays. - Wrap all Pyccel kernels (stencil, B-splines, field evaluation, DOF kernels) in cunumpy.PyccelKernel, so they accept CuPy arrays. - Keep host-only metadata on NumPy: MPI/index bookkeeping in ddm (cart, partition, petsc) and fem.partitioning, Kronecker solver sizes, and index arithmetic with Python ints (compute_diag_len, math.prod). - Stage data for host-only libraries by array, not by global backend: LAPACK/SuperLU direct solvers, SciPy FFT, SciPy sparse products. - Fix calls that ran host kernels on device arrays: the second stencil2coo call in StencilMatrix.tosparse and the conjugate transpose. - Vectorize the construction of the 1D collocation matrices in the global projectors (element-wise indexing was a device round trip per entry: 334 s of a 348 s Derham setup on the GPU). - GMRES: take real scalars from CuPy views before modifying them. - Fix StencilMatrix._update_ghost_regions_serial: the ghost region is pads * shifts wide (wrong whenever shifts > 1, on both backends). - Tests: work with CuPy arrays; skip PETSc tests without petsc4py. Serial tests pass on both backends (core, ddm, fem, linalg). Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
This was referenced Oct 1, 2026
Draft
max-models
marked this pull request as ready for review
October 1, 2026 14:01
9 of 20 tasks
max-models
added a commit
to struphy-hub/struphy
that referenced
this pull request
Oct 2, 2026
**Solves the following issue(s):** PR 7 of the CUDA strategy (`CUDA_STRATEGY.md`). No issue closed. Stack: #649 → #651 → this PR → #653. The commits `Derham on the GPU: create Derham on CuPy, add Derham.cuda_args_derham` and `CUDA_STRATEGY: PR 7 depends on feectools#85` are new relative to #651. Depends on struphy-hub/feectools#85 (feectools on the CuPy backend). Without it, feectools itself fails while the spline spaces are built on CuPy. **Core changes:** `Derham` (`src/struphy/feec/psydac_derham.py`) can be created on the CuPy backend. The rule: data that describes the spline spaces stays on the host on every backend, because it goes to pyccel kernels, SciPy or MPI. Only coefficients and stencil matrices live on the device. - The projection and quadrature grids (`get_pts_and_wts`, `get_pts_and_wts_quasi`, `get_span_and_basis`, the local-projector weights) and `spline_types_pyccel` are NumPy arrays. This removes the ad-hoc `cupy` conversions in them. - `domain_array`, `index_array`, `index_array_N`, `index_array_D` and `neighbours` are gathered with NumPy MPI buffers, so no CUDA-aware MPI is needed. They are then converted with `xp.asarray`, so they are device arrays on CuPy, like `Particles.domain_array`. The neighbour search runs on the host because it builds an object array with `None` entries, which CuPy cannot hold. - `args_derham` is built from the host knots, degrees and starts directly; `_to_numpy_for_kernel` is gone. - New `Derham.cuda_args_derham`: `CudaDerhamArguments` (new in `src/struphy/utils/cuda_arguments.py`) with one device copy of the degrees, knots and start indices, built on first access. On the NumPy backend it raises `TypeError`, since host arrays are never copied to the device. The pyccel scratch arrays (`bn1`, ..., `bd3`) are not part of it; they become per-thread local arrays in CUDA (PR 10). - Local projectors (`DerhamOptions.local_projectors=True`) raise `NotImplementedError` on CuPy when the `Derham` is created, instead of failing later inside `CommutingProjectorLocal`. Not in this PR: - Local projectors on CuPy. `CommutingProjectorLocal` builds its data with `xp` and calls pyccel kernels on it, as `Derham` did. - Polar splines on CuPy. They need a spline mapping, which cannot be created on CuPy yet (see PR 5). - Field evaluation (`SplineFunction.__call__`, ...) on the device. It still calls pyccel kernels with the coefficients, so it needs CUDA evaluation kernels (PR 10+). Tests: - New `src/struphy/feec/tests/test_derham_gpu.py`: - `test_cuda_args_derham_needs_device_arrays` (runs everywhere): raises on the NumPy backend. - `test_derham_on_cupy` (GPU only), with periodic and with Dirichlet/free boundaries: the decomposition tables and kernel arguments agree between the backends, and `cuda_args_derham` is built once and holds device copies of the right values. - `test_local_projectors_not_supported_on_cupy` (GPU only). - The GPU-only tests have not run on a GPU yet. On the CPU they were run with a strict host stand-in for CuPy, which rejects host/device mixing and cannot run kernels, on top of feectools#85. With it, a `Derham` created on the "CuPy" backend matches the NumPy one on 1, 2 and 4 MPI processes. - `test_kernel_backends.py::test_pusher_accepts_kernel` fails with and without this PR (`ValueError: array is too big` in `draw_markers`). It is unrelated. **Model-specific changes:** None 🤖 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
1 participant
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.
Summary
feectools now works when cunumpy's backend is CuPy (
ARRAY_BACKEND=cupy). No device kernels yet: Pyccel kernels still run on the host, and their arrays are copied there and back. This step only makes the CuPy backend correct. Speed comes in #88.Changes
cunumpy.PyccelKernel, so they accept CuPy arrays.ddm(cart,partition,petsc) andfem.partitioning, Kronecker solver sizes, and index arithmetic (compute_diag_len,math.prod).stencil2coocall inStencilMatrix.tosparse, and the conjugate transpose.StencilMatrix._update_ghost_regions_serialnow uses a ghost regionpads * shiftswide. It was wrong whenevershifts > 1.petsc4pyis missing.Testing
Serial tests in
core,ddm,femandlinalgpass on the NumPy and CuPy backends.Stack
This is step 1 of 4 for CUDA support. The PRs are stacked, and each one targets
devel-tiny. This PR is the base, so it contains only commit0a55cbd. Merge them in order:🤖 Generated with Claude Code