Cuda pr 6 particles on gpu - #651
Merged
Merged
Conversation
The SPH linear smoothing kernels returned a non-zero gradient at zero separation, so every particle pushed on itself. This PR makes all SPH kernel gradients vanish at r=0, which matches the symmetric value the trigonometric and gaussian kernels already give. Closes #437 ### What was wrong - `grad_linear_uni(0, h)` returned `+1/h**2`, the left derivative. That value carries into `grad_linear_1d`, `grad_linear_2d_{1,2}` and `grad_linear_3d_{1,2,3}`, and it applies whenever the component's own coordinate is 0, not only at the origin. - `grad_linear_isotropic_3d_{1,2,3}` returned `-1/h / (C h^3)` at `r == 0`. - As a result the SPH gradient of a lattice-loaded constant density was not zero: it came out at about 3.5 for 1d linear, 3.8 for 3d tensor linear and -0.58 for isotropic in the check below. ### What changed - `src/struphy/pic/sph_smoothing_kernels.py`: at the cusp the gradient now returns 0 (`x == 0` in `grad_linear_uni`, `r == 0` in the isotropic variants). Values away from 0 are unchanged. The trigonometric and gaussian gradients were checked too and already give 0 there. - `src/struphy/pic/tests/test_kernel_setup.py`: - New test `test_sph_kernel_gradients_vanish_at_zero` covers every gradient type in `smoothing_kernel`, at the origin and on the coordinate planes. - `test_sph_tensor_destinations` relied on the old self-gradient `1/h**2`. It now uses a second neighbouring particle so the viscosity-tensor destinations still get non-trivial values. ### Verification - A standalone check calls all 21 gradient kernel types at the origin and on the coordinate planes. Before the fix, 9 types were non-zero at the origin and all 4 plane cases were non-zero. After the fix all of them are 0. At random points away from 0 the values are bit-identical before and after. - The gradient of a lattice-loaded constant density is now ~1e-16 in 1d linear, 3d tensor linear and 3d isotropic. - `pytest src/struphy/pic/tests/test_kernel_setup.py`: 74 passed. The new test fails without the fix. - A subset of `test_sph.py::test_sph_evaluation_1d` (linear_1d, periodic) passes. ### Risks / follow-ups - SPH results that use the linear kernels change slightly, because the self-contribution is gone. That removal is the intended correction. - #436 (viscosity strain without DF^-T) is in the same area and is not addressed here. - This branch includes the commit bumping `feectools` to 0.1.11 so that CI can run. 🤖 Generated with [Claude Code](https://claude.com/claude-code) --------- Co-authored-by: Claude Opus 5.5 <noreply@anthropic.com>
…ld just be a proof of concept
CudaKernel no longer flattens argument objects or casts scalars at call time: the arguments must already be in the format of a cupy.RawKernel (flat tuple of CuPy arrays and NumPy scalars, built once at setup), and the number of threads is passed as n_threads. Add docstrings to the kernel and argument classes. Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
Converting the scalars and joining the .values of the CUDA argument classes costs about 1 us per call (the kernel launch alone about 70 us), prevents silently misaligned arguments when a Python int is passed as a 64-bit value, and keeps the calls the same on both backends. Arrays are still never converted or copied; n_threads stays an explicit argument. Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
…sses CUDA reads each kernel argument with the size declared in the signature, so Python int/float already arrive correctly in int/double parameters; the casts did not change the result, and did not prevent the actual failure cases (e.g. an integer passed to a double parameter). Checking scalars against the kernel signature is noted as a follow-up in CUDA_STRATEGY.md. Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
Add CudaKernel.from_file, which reads the CUDA source from a <name>_cuda.cu file and takes the kernel name from the file name, and ship .cu/.cuh files as package data. Step 2 of CUDA_STRATEGY.md. Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
KernelCatalog.from_package pairs <name>/<name>_kernels.py (pyccel) with <name>/<name>_cuda.cu (CUDA, optional). The CUDA kernel of a Kernel is now optional: on the CuPy backend, a kernel without a CUDA version raises NotImplementedError naming the expected .cu file. Step 3 of CUDA_STRATEGY.md. Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
Pusher now takes a Kernel or, as before, a PyccelKernel (wrapped into a Kernel without CUDA version) and calls get_kernel() in its constructor. On the CuPy backend, a pusher whose kernel has no CUDA version thus fails when it is created instead of in the time loop. No changes to the propagators, no behaviour change on the CPU. Step 4 of CUDA_STRATEGY.md. Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
- Test only analytic mappings on the CuPy backend: spline mappings such as IGAPolarCylinder cannot be created there yet (interp_mapping passes CuPy arrays to scipy.sparse.csc_matrix); correct CUDA_STRATEGY.md. - Build the CUDA domain arguments with cupy instead of xp, so they can be built whichever backend is active (the arrays are on the device); test this. - Move the new tests above the __main__ block of test_domain.py. Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
9 of 20 tasks
max-models
marked this pull request as ready for review
October 1, 2026 22:26
max-models
enabled auto-merge (squash)
October 1, 2026 22:26
spossann
approved these changes
Oct 2, 2026
max-models
added a commit
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>
max-models
added a commit
that referenced
this pull request
Oct 2, 2026
**Solves the following issue(s):** PR 8 of the CUDA strategy (`CUDA_STRATEGY.md`). No issue closed. Stack: #649 → #651 → #652 → this PR. Only the commit `Shared CUDA header for the kernel argument classes` is new relative to #652 (plus a merge of #652). **Core changes:** CUDA kernels used to repeat the full flat signature of the argument classes (26 parameters for markers and domain alone). CuPy does not check kernel signatures, so adding a field to `MarkerArguments` would have shifted all following arguments of every CUDA kernel without an error. Now each argument class is one C struct, defined once. - New `src/struphy/kernel_arguments/pusher_args.cuh` with `struct MarkerArgs`, `DerhamArgs` and `DomainArgs`. Kernels take them by value in place of the pyccel argument classes, e.g. `void push_eta_linear(double dt, int stage, MarkerArgs args_markers, DomainArgs args_domain)`, and read `args_markers.markers`, `args_markers.n_markers`, and so on. - The member names are the attribute names of the pyccel classes. The only CUDA-specific member is `MarkerArgs.n_cols`, since pyccel kernels take it from `markers.shape[1]`. Accordingly, `CudaMarkerArguments.first_pusher_idx` is now `first_init_idx`, as in `MarkerArguments`. - `src/struphy/utils/cuda_arguments.py`: - Each `Argument` subclass lists its members once, e.g. `fields = (("double*", "markers"), ...)`. - `struct_dtype()` builds a NumPy structured dtype with `align=True`, which gives C alignment and padding. Pointer members are `uint64` device addresses (`array.data.ptr`). - The struct is a `numpy.void`, which CuPy passes by value. It is packed once in the constructor, so kernel calls do no extra work. - Scalar members are checked when the struct is packed. A `float` for an `int` member raises `TypeError`, and a value that does not fit raises `OverflowError`. Before, they arrived in the kernel truncated or wrapped around; plain NumPy assignment also truncates `1.7` to `1`. This solves the "scalar types" open question for everything inside the argument classes. Scalars passed directly to a kernel (`dt`, `stage`) are still unchecked. - Because the struct holds device addresses, it is repacked when an argument object is deepcopied or unpickled (`__getstate__`/`__setstate__`). - `CudaKernel` (`src/struphy/utils/kernel_backends.py`) compiles with the folder that contains the `struphy` package on the NVRTC include path, so `.cu` files can `#include "struphy/kernel_arguments/pusher_args.cuh"`. `.cuh` files are already package data. **Why structs, and why `_pack()` and the copy handling are needed:** The underlying problem: CuPy passes kernel arguments positionally and never checks them against the signature. With flat parameter lists, a new field in an argument class shifts every later argument of every CUDA kernel, silently and on the GPU only. Some single, checked definition of each argument layout is therefore necessary. A C struct is the simplest one that also works for device helper functions. What each piece is for: - `struct_dtype()` with `align=True`: CuPy passes a struct by value only as a NumPy structured scalar, and `align=True` gives the C padding (the `int` members are followed by pointers). - Pointers as `uint64` (`data.ptr`): a struct cannot hold a CuPy array, only its device address. - Repacking in `__setstate__`: a direct consequence of storing addresses. The arrays of a deepcopied or unpickled object have new addresses, and without the repack a copy would silently read the original's arrays. - Packing once in the constructor: so kernel calls do no extra work. Packing at every call would also be fine (well under 1 µs next to a launch of about 70 µs). - Scalar checks: the only optional part. They are cheap and turn truncation or wrap-around into an exception. For: - Silent argument shifts cannot happen, and a test that needs no GPU checks the header against `fields`. - CUDA code mirrors the pyccel code: pyccel kernels and helpers take `args_domain` and read `args_domain.kind_map`, and CUDA reads exactly the same. This decides it for PR 10. The mapping evaluation, B-spline and boundary helpers take these argument objects, so a `__device__` helper takes one `const DomainArgs&` instead of 12+ parameters forwarded at every call site. - Short kernel signatures. Adding a field touches `fields` and the header, not every kernel. - No runtime cost. Against: - Raw pointers bypass CuPy's per-launch checks: a host array or an array on another GPU is no longer rejected at launch. `_cupy_array()` checks at construction that the array is a C-contiguous CuPy array of the right dtype, but not which device it is on. That will matter with one GPU per MPI rank; adding a device check is a small follow-up. - If an owner reallocates an array without rebuilding its argument object, the struct points at the old array. The pyccel argument classes have the same issue, since they also hold references, so this is not new. - About 100 lines of infrastructure plus layout tests before any real CUDA kernel exists. But this choice fixes the signatures of all device helpers in PR 10, so changing it later would mean rewriting them. - Alignment is subtle. The host `sizeof`/`offsetof` check and the NVRTC layout test cover it. Alternatives considered: - Flat parameters generated from a C macro (`#define MARKER_ARGS double* markers, ...`), with `get_cuda_args()` returning the arrays in `fields` order. This also prevents silent shifts and keeps CuPy's checks, without pointers or repacking. But every device helper then needs a long parameter list, or a second macro that rebuilds a struct inside the kernel. That is more macro machinery, and the CUDA code no longer mirrors the pyccel helpers. - Plain flat lists without a single definition: the state before this PR, and the dangerous one. Tests (`src/struphy/pic/tests/test_kernel_backends.py`): - Run everywhere: - `test_cuda_argument_structs_match_header`: parses `pusher_args.cuh` and compares member names, C types and order with `fields`. - `test_cuda_struct_members_are_pyccel_attributes`: the member names exist as attributes in `pusher_args_kernels.py`, except `n_cols`. - GPU only: - `test_cuda_struct_layout`: an NVRTC-compiled kernel reports `sizeof` and the member offsets and sizes, which are compared with the NumPy dtype. - `test_cuda_struct_scalars_are_checked` and `test_cuda_struct_follows_copies`. - The demo kernels `push_eta_linear` and `write_scalars` now use the structs, and `write_scalars` also reads a pointer member of each struct. - Fixed: several GPU-only tests here and in `geometry/tests/test_domain.py` still used the `values` attribute, which `get_cuda_args()` replaced. They would have failed on a GPU. The GPU-only tests have not run on a GPU yet. On the host: - Compiling the header with a C++ compiler gives the same `sizeof`/`offsetof` as the NumPy dtypes. x86-64 and arm64 lay these structs out like CUDA. - The tests that do not launch kernels pass with a strict host stand-in for CuPy. - `test_pusher_accepts_kernel` fails with and without this PR (see #652). 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
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.
Particles can now be initialized with the CuPy backend. The new
cuda_args_markersproperty lazily builds and cachesCudaMarkerArgumentsfrom the particle’s device arrays, validity mask, indices, and boundary-condition codes. Existingargs_markersremains available for Pyccel calls.The initialization path now converts Python process-count lists to
cunumpyarrays before reductions. Scalar MPI gathers use small NumPy buffers, then copy the results to the active backend.GPU particle creation and CUDA argument construction still need a GPU smoke test.