Cuda pr 8 cuda headers

CUDA_STRATEGY.md

@@ -13,7 +13,7 @@ The work is split into small PRs that can be reviewed and merged one at a time.
 - [x] **PR 5: `Domain` on the GPU** — domain arguments are selected and stored at domain construction; CUDA arguments reference device arrays, and deepcopy/unpickling rebuilds the arguments from the copied or restored arrays.
 - [x] **PR 6: `Particles` on the GPU** — `Particles` can be created on the CuPy backend, and `Particles.args_markers` is selected as the CUDA or Pyccel argument bundle at construction.
 - [x] **PR 7: `Derham` on the GPU** — `Derham` can be created on the CuPy backend, plus `Derham.cuda_args_derham`.
-- [ ] **PR 8: Shared CUDA headers for the argument classes** — one `.cuh` per argument class instead of long flat kernel signatures.
+- [x] **PR 8: Shared CUDA header for the argument classes** — one struct per argument class in `kernel_arguments/pusher_args.cuh`, passed by value, instead of long flat kernel signatures.
 - [ ] **PR 9: One folder per kernel, starting with `pic/pushing`** — pure refactor, no behaviour change.
 - [ ] **PR 10: Device versions of helper kernels** — B-spline evaluation, mapping evaluation (per domain), small linear algebra, as `__device__` functions in `.cuh` headers.
 - [ ] **PR 11: First real CUDA kernel** — `push_eta_stage` with a pyccel/CUDA parity test and an end-to-end run on the GPU.
@@ -41,12 +41,13 @@ CUDA kernels can be added one by one. If the code runs on the GPU and needs a ke
 - **No silent CPU fallback on the GPU.** A kernel without a CUDA version raises an error on the GPU backend. Falling back would mean copying data to the host and back at every call.
 - **Small steps.** Every PR keeps the CPU code path working and tested.
 
-## Current state (PR 7)
+## Current state (PR 8)
 
 | File | Content |
 |---|---|
-| `src/struphy/utils/kernel_backends.py` | `is_cuda_backend()`, `CudaKernel` (wraps a `cupy.RawKernel`, compiled lazily; expands `Argument.get_cuda_args()` and takes `n_threads`), `Kernel` and `KernelCatalog` for backend selection and discovery |
-| `src/struphy/utils/cuda_arguments.py` | `Argument` contract plus `CudaMarkerArguments`, `CudaDerhamArguments` and `CudaDomainArguments`; CUDA arrays are stored individually and returned in signature order by `get_cuda_args()` |
+| `src/struphy/utils/kernel_backends.py` | `is_cuda_backend()`, `CudaKernel` (wraps a `cupy.RawKernel`, compiled lazily with the struphy headers on the include path; expands `Argument.get_cuda_args()` and takes `n_threads`), `Kernel` and `KernelCatalog` for backend selection and discovery |
+| `src/struphy/utils/cuda_arguments.py` | `Argument` base class plus `CudaMarkerArguments`, `CudaDerhamArguments` and `CudaDomainArguments`; each references its device arrays and packs them once into its C struct, which `get_cuda_args()` returns |
+| `src/struphy/kernel_arguments/pusher_args.cuh` | the C structs `MarkerArgs`, `DerhamArgs` and `DomainArgs` that CUDA kernels take in place of the pyccel argument classes |
 | `src/struphy/geometry/base.py` | `Domain.args_domain` is selected once at construction; CUDA domains use device arrays, while direct Pyccel geometry calls retain a host argument bundle |
 | `src/struphy/pic/base.py` | `Particles` arrays and `args_markers` use the backend selected at construction; direct Pyccel methods retain a private host bundle |
 | `src/struphy/pic/tests/test_kernel_backends.py` | the demo kernel pair `push_eta_linear` (pyccel function compiled with `epyccel` at test time, CUDA source string) and tests on both backends |
@@ -57,7 +58,7 @@ Things we learned in the proof of concept:
 - `Particles.args_markers` is the CUDA or Pyccel bundle selected at construction. Existing direct Pyccel methods use a private host bundle. `Derham` still needs its CUDA creation work (PR 7).
 - `Domain.args_domain` returns the argument type selected at domain construction; on CuPy it is a `CudaDomainArguments` object. Direct Pyccel methods on `Domain` use a separate internal host bundle.
 - `cupy.RawKernel` accepts only device arrays (host arrays raise) and does **not** check the kernel signature. Each argument is read with the size declared in the signature, so Python `int`/`float` arrive correctly in `int`/`double` parameters, but a wrongly typed scalar (e.g. an integer for a `double`, or a value that overflows an `int`) gives a wrong value **without an error**. Casting Python scalars in `CudaKernel` does not prevent this, so it is not done; see the follow-up in [Open questions](#open-questions).
-- Flattening the argument classes at each call (joining their `values`) costs well under 1 µs, compared to about 70 µs for launching the kernel.
+- Flattening the argument classes at each call costs well under 1 µs, compared to about 70 µs for launching the kernel. Since PR 8 each class is one struct, packed once when the argument object is created.
 - `struphy compile` compiles every `.py` file whose name contains `kernels`. Non-pyccel modules must not contain `kernels` in their name; `.cu` files are ignored by it.
 - On an H100, the demo kernel pushes 10⁶ markers in about 0.13 ms per step.
 
@@ -100,7 +101,7 @@ kernel = catalog["push_eta_stage"]  # Kernel: pyccel or CUDA depending on the ba
 
 - `CudaKernel.from_file(path)` reads `<name>_cuda.cu`; the kernel name is taken from the file name. The kernel is compiled lazily on first call. CuPy caches compiled kernels on disk (`~/.cupy/kernel_cache`), so the compile cost is paid once per machine.
 - Add `"**/*.cu"` and `"**/*.cuh"` to `[tool.setuptools.package-data]` in `pyproject.toml`.
-- Later (PR 10), when the first shared header is needed: headers are found through NVRTC include paths (`cupy.RawModule(code=..., options=("-I<struphy src>",))`), so a `.cu` file can `#include "struphy/bsplines/bsplines_kernels.cuh"`.
+- Shared headers are found through the NVRTC include path (done in PR 8: `CudaKernel` compiles with `-I<folder containing the struphy package>`), so a `.cu` file can `#include "struphy/kernel_arguments/pusher_args.cuh"` or, later, `"struphy/bsplines/bsplines_kernels.cuh"`.
 
 ### PR 3: Kernel catalog
 
@@ -156,13 +157,16 @@ kernel = catalog["push_eta_stage"]  # Kernel: pyccel or CUDA depending on the ba
   - Field evaluation (`SplineFunction.__call__`, ...) still calls pyccel kernels with the coefficients, which are device arrays on CuPy. It needs CUDA evaluation kernels (PR 10+).
 - Tests: `feec/tests/test_derham_gpu.py`. Without a GPU, a strict host stand-in for CuPy (rejects host/device mixing, cannot run kernels; not part of the repository) was used, on top of struphy-hub/feectools#85. With it, a `Derham` created on the "CuPy" backend matches the NumPy one on 1, 2 and 4 MPI processes.
 
-### PR 8: Argument structs in shared headers
+### PR 8: Argument structs in a shared header (complete)
 
-Today every CUDA kernel repeats the full flat signature (26 parameters for markers and domain alone). CuPy does not check it, so adding a field to `MarkerArguments` would shift all following arguments of all CUDA kernels **without an error**.
+Before, every CUDA kernel repeated the full flat signature (26 parameters for markers and domain alone). CuPy does not check it, so adding a field to `MarkerArguments` would have shifted all following arguments of all CUDA kernels **without an error**.
 
-- Define `struct MarkerArgs { double* markers; bool* valid_mks; int n_markers; ... };` etc. in `kernel_arguments/pusher_args.cuh`, and pass one struct per argument class.
-- To check first: how to pass a struct to a `cupy.RawKernel` (e.g. as a NumPy structured scalar with pointer fields). If this does not work well, keep the flat signature, but generate it from the Python class so that it is defined in one place.
-- A test compares the struct layout (field names, types, order) with the Python class.
+- `kernel_arguments/pusher_args.cuh` defines `struct MarkerArgs`, `DerhamArgs` and `DomainArgs`. A kernel takes them by value in the place of the pyccel argument classes, e.g. `void push_eta_linear(double dt, int stage, MarkerArgs args_markers, DomainArgs args_domain)`, and reads `args_markers.markers`, `args_markers.n_markers`, ... The member names are the attribute names of the pyccel classes; the only CUDA-specific member is `MarkerArgs.n_cols` (pyccel takes `markers.shape[1]`).
+- Passing a struct: CuPy passes a NumPy scalar by value, copying `itemsize` bytes. Each `Argument` subclass lists its members as `fields = (("double*", "markers"), ...)`, from which `struct_dtype()` builds a NumPy structured dtype with `align=True` (C alignment and padding). Pointer members are `uint64` holding `array.data.ptr`. The struct (a `numpy.void`) is packed once, in the constructor, so a kernel call does no extra work.
+- The struct holds device addresses: it is repacked when an argument object is deepcopied or unpickled (`__getstate__`/`__setstate__`), and owners that reallocate an array must rebuild their argument object, as before.
+- Scalar members are checked when packing: a `float` for an `int` member raises `TypeError` and a value that does not fit raises `OverflowError`, instead of arriving truncated or wrapped around (NumPy assignment alone truncates `1.7` to `1`).
+- Tests: the header is parsed and compared with `fields` (names, C types, order) without a GPU; on a GPU, an NVRTC-compiled kernel reports `sizeof` and the member offsets, which are compared with the dtype. The demo kernels use the structs. On the host, `offsetof`/`sizeof` from a C++ compiler agree with the dtypes (x86-64/arm64 lay these structs out like CUDA).
+- Fixed along the way: the GPU-only tests still used the `values` attribute that was replaced by `get_cuda_args()`.
 
 ### PR 9: One folder per kernel
 
@@ -211,7 +215,7 @@ Port the kernels in the order the target models need them, so that complete mode
 
 ## Open questions
 
-- **Scalar types.** Scalars are not checked against the kernel signature (see [Current state](#current-state-pr-1)). `CudaKernel` could read the parameter types from the `extern "C"` signature once, when it is created, and cast each scalar to its declared type or raise if it does not fit (e.g. a Python `float` for an `int`, or an overflowing integer). This could go together with PR 8.
+- **Scalar types.** Since PR 8 the members of the argument structs are checked when they are packed. Scalars passed directly to a kernel (e.g. `dt`, `stage`) are still not checked against the kernel signature (see [Current state](#current-state-pr-8)). `CudaKernel` could read the parameter types from the `extern "C"` signature once, when it is created, and cast each scalar to its declared type or raise if it does not fit (e.g. a Python `float` for an `int`, or an overflowing integer).
 
 - **Marker layout.** The markers array is row-major (`n_markers × n_cols`). With one thread per marker, the memory accesses are strided. This is fine for now (each thread reads a few neighbouring columns), but a column-major or struct-of-arrays layout may be faster later. This would affect the CPU code too, so it is out of scope here.
 - **MPI + GPUs.** One GPU per MPI rank (`cunumpy.set_device(rank % n_gpus)`), and GPU-aware MPI for the marker exchange, so markers do not go through the host.

src/struphy/geometry/tests/test_domain.py

@@ -1135,26 +1135,18 @@ def test_cuda_args_domain(mapping):
         assert isinstance(args, CudaDomainArguments)
         assert domain.args_domain is args  # built once
 
-        kind_map, params, degree, t1, t2, t3, ind1, ind2, ind3, cx, cy, cz = args.values
-        assert int(kind_map) == domain.kind_map
+        assert args.kind_map == domain.kind_map
         # no copies of arrays that already have the right dtype and layout
-        assert t1 is domain.T[0] and ind3 is domain.indN[2]
+        assert args.t1 is domain.T[0] and args.ind3 is domain.indN[2]
 
+        # the struct holds the device addresses of these arrays
+        (struct,) = args.get_cuda_args()
+        assert struct["kind_map"] == domain.kind_map
         host = domain._pyccel_args_domain
-        for dev, ref in (
-            (params, host.params),
-            (degree, host.degree),
-            (t1, host.t1),
-            (t2, host.t2),
-            (t3, host.t3),
-            (ind1, host.ind1),
-            (ind2, host.ind2),
-            (ind3, host.ind3),
-            (cx, host.cx),
-            (cy, host.cy),
-            (cz, host.cz),
-        ):
-            assert (cunumpy.to_numpy(dev) == ref).all()
+        for name in ("params", "degree", "t1", "t2", "t3", "ind1", "ind2", "ind3", "cx", "cy", "cz"):
+            dev = getattr(args, name)
+            assert struct[name] == dev.data.ptr, name
+            assert (cunumpy.to_numpy(dev) == getattr(host, name)).all(), name
 
 
 @requires_cupy
@@ -1172,7 +1164,8 @@ def test_domain_deepcopy_and_pickle_on_cupy(mapping):
             assert (other.args_domain.params == domain.args_domain.params).all()
             other_cuda = other.args_domain
             assert other_cuda is not cuda_args
-            assert other_cuda.values[3] is other.T[0]
+            assert other_cuda.t1 is other.T[0]
+            assert other_cuda.get_cuda_args()[0]["t1"] == other.T[0].data.ptr
 
 
 if __name__ == "__main__":

src/struphy/kernel_arguments/pusher_args.cuh

@@ -0,0 +1,57 @@
+// CUDA versions of the argument classes in pusher_args_kernels.py, one struct per class.
+//
+// CUDA kernels take these structs by value, in the place of the pyccel argument classes, e.g.
+//
+//     #include "struphy/kernel_arguments/pusher_args.cuh"
+//
+//     extern "C" __global__
+//     void push_eta_stage(double dt, int stage, MarkerArgs args_markers, DomainArgs args_domain, ...)
+//
+// The structs are filled on the host by the classes in struphy/utils/cuda_arguments.py, whose `fields` list the
+// members below in the same order and with the same types; test_cuda_argument_structs checks that they agree.
+// The member names are the attribute names of the pyccel classes. Pointers are device pointers.
+#pragma once
+
+// CUDA version of MarkerArguments (struphy.utils.cuda_arguments.CudaMarkerArguments).
+struct MarkerArgs {
+    double* markers;  // (n_markers, n_cols), row-major
+    bool* valid_mks;  // (n_markers,), true for markers that are neither holes nor ghosts
+    int n_markers;
+    int n_cols;
+    int Np;
+    int vdim;
+    int weight_idx;
+    int first_diagnostics_idx;
+    int first_init_idx;
+    int first_shift_idx;
+    int residual_idx;
+    int first_free_idx;
+    int mu_idx;
+    long long* bc_type;  // (3,)
+};
+
+// CUDA version of DerhamArguments (struphy.utils.cuda_arguments.CudaDerhamArguments).
+// The scratch arrays of the pyccel class (bn1, ..., bd3) are local arrays in the kernels.
+struct DerhamArgs {
+    long long* pn;  // (3,)
+    double* tn1;
+    double* tn2;
+    double* tn3;
+    long long* starts;  // (3,)
+};
+
+// CUDA version of DomainArguments (struphy.utils.cuda_arguments.CudaDomainArguments).
+struct DomainArgs {
+    int kind_map;
+    double* params;
+    long long* degree;  // (3,)
+    double* t1;
+    double* t2;
+    double* t3;
+    long long* ind1;  // (number of mapping grid cells, degree + 1)
+    long long* ind2;
+    long long* ind3;
+    double* cx;  // control points
+    double* cy;
+    double* cz;
+};