MHD equilibria on the CuPy backend by max-models · Pull Request #714 · struphy-hub/struphy

max-models · 2026-10-07T19:12:00Z

Stack: part 7 of 14, based on #713 (merge that first), next: #718. Full order: #709 → #708 → #711 → #712 → #705 → #713 → #714 → #718 → #720 → #721 → #724 → #725 → #723 → #722.

Solves the following issue(s):

Solves #696, part of #650 (see CUDA_STRATEGY.md). Every equilibrium in fields_background/equils.py can now be created and evaluated on the CuPy backend. The time loop does not change: it uses the projected equilibrium.

Core changes:

What failed on CuPy. I checked this with cunumpy's fake CuPy in a clean serial subprocess. Every method models call was evaluated on a meshgrid and at markers.
- Construction failed for:
  - AdhocTorus with q_kind 1 or 2: SciPy quad/UnivariateSpline got device arrays.
  - AdhocTorusQPsi: odeint/fsolve got device arrays.
  - EQDSKequilibrium: RectBivariateSpline got the device arrays from xp.linspace.
- Evaluation failed for every field of GVECequilibrium and DESCequilibrium: gvec and DESC/JAX got device arrays.
- Already working: the analytic ones (HomogenSlab, ShearedSlab, ShearFluid, ScrewPinch, AdhocTorus with q_kind=0, CircularTokamak, ConstantVelocity, HomogenSlabITG, CurrentSheet, the generic ones).
- Also found: AdhocTorus.psi_r (q_kind=0) did out *= ... on a 0-d device array, which fails on CuPy when psi is called with arrays. It is now out = out * ....
New helpers in fields_background/base.py:
- setup_on_host: a decorator for __init__. It runs the setup on the NumPy backend (file reading, ODE solves, SciPy fits), so these equilibria hold only host data on either backend. Used on AdhocTorus, AdhocTorusQPsi and EQDSKequilibrium.
- host_call(fun, *args, **kwargs): calls a host-only function, copying device arguments to the host and array results back once per call. Without device arguments it is a plain call, so the NumPy path does not change. Used for every SciPy spline evaluation: psi_r/p_r of AdhocTorus, psi_r of AdhocTorusQPsi, and psi/q_psi/g_psi/p_psi of EQDSK.
- evaluate_on_host: the same as a method decorator. Used on the GVEC bv, jv, p0, n0 and the DESC bv, jv, p0, n0, gradB1.
Cost of the host evaluation. Each call costs one host round trip. For example, gradB_xyz of EQDSK makes about 10 psi calls, so 10 round trips. This is fine because equilibria are only evaluated at setup: initial conditions, projections onto the FEEC spaces and mappings. A device B-spline evaluation of the SciPy knots and coefficients would remove these copies. It is listed as an open question in CUDA_STRATEGY.md.
Tokamak: the workaround from PR 19 is removed. It no longer builds its default EQDSKequilibrium inside use_backend("numpy"). The field-line tracing (SciPy) still runs on the NumPy backend.
Tests: new file fields_background/tests/test_equils_cupy.py.
- What it checks: 19 cases that cover every equilibrium class of equils, including EQDSK, GVEC and DESC with the files shipped in the repo or with the packages. test_all_equilibria_have_cases checks that no class is missing. Each case is created on NumPy and on CuPy, and 30 logical-domain methods (absB0, p0, n0, b2, unit_b1, gradB1, curl_unit_b2, j2, …) are evaluated:
  - on a meshgrid and at markers,
  - plus psi/g_tor with all derivatives for the axisymmetric equilibria.
  Every array result must be a device array and match NumPy to 1e-12.
- test_equil_fake_cupy runs without a GPU, in a subprocess on the fake CuPy. The fake CuPy cannot launch CUDA kernels, so host_geometry_kernels() runs the four geometry entry kernels with their pyccel version on the host buffers of the fake arrays. The CUDA geometry kernels have their own parity and emulation tests. The emulated_launches() of CUDA version of the linear_vlasov_ampere accumulation #705 would also work here, but it compiles every launch, and these cases make hundreds of domain evaluations.
- test_equil_cupy runs the same check on a GPU.

Model-specific changes:

None.

Documentation changes:

CUDA_STRATEGY.md:

New section "MHD equilibria on CuPy (MHD equilibria on the CuPy backend #696)".
The open question "Polar splines and MHD equilibria" now covers polar splines only.
New open question: equilibrium splines on the device.

Testing (macOS, no GPU; pyccel kernels compiled with GNU/Fortran; GPU tests not run):

	before (`devel`)	after
`fields_background/tests geometry/tests`	92 passed, 12 skipped	113 passed, 31 skipped
`kinetic_background/tests post_processing/tests/test_orbits_tools.py feec/tests/test_field_init.py`	50 passed	50 passed

The 21 new passes are the 19 fake-CuPy cases, test_all_equilibria_have_cases and test_host_call_numpy_is_plain_call. The 19 new skips are the GPU tests (test_equil_cupy). The DESC case takes about 35 s, because it loads W7-X from the desc package.

🤖 Generated with Claude Code

Equilibria with host-only setup (AdhocTorus q_kind 1/2, AdhocTorusQPsi, EQDSKequilibrium) run their __init__ on the NumPy backend (setup_on_host) and hold only host data. SciPy spline evaluations go through host_call, GVEC/DESC evaluations through @evaluate_on_host: device arguments are copied to the host and the result back, once per call; NumPy arguments are evaluated as before. Tokamak no longer builds its default EQDSKequilibrium on the NumPy backend itself. New test_equils_cupy.py compares every equilibrium of equils on CuPy (GPU, or cunumpy's fake CuPy with host geometry kernels) with NumPy. Solves #696. Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>

Stack the PR on #713. Conflict in CUDA_STRATEGY.md (porting-order gates): kept #714's removal of the MHD equilibria gate and #705's note on 6D views; marker sorting is on the device since #712. Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>

Stack the PR on #714. Conflicts in CUDA_STRATEGY.md: kept the linear_vlasov_ampere checklist entry and #690's decision in the Next entry (without the steps done below), the marker exchange note of #712, and the polar-splines open question of #714. Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>

… child tests - GVECequilibrium evaluates markers point by point: the coordinates go to gvec as xarray DataArrays with a shared dimension, so bv/jv at N markers no longer return N x N x N arrays and absB0 at markers works (also on NumPy). With use_boozer=True, marker evaluation raises, since gvec computes the Boozer transform per flux surface. New test_gvec_equil.py. - Use xp.host_call, xp.evaluate_on_host and xp.setup_on_host from cunumpy 0.6.2 instead of the copies in fields_background/base.py; pin cunumpy 0.6.2. - run_fake_cupy_child: under mpirun only rank 0 starts the fake-CuPy child (every rank started the same child at once, and on CI one of two concurrent GVEC children died with SIGILL and no output); the child runs with faulthandler and one OpenMP thread, and failures report the signal and the end of stdout/stderr. Used by the four fake-CuPy subprocess tests. Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>

…GVEC tests GVEC dies with SIGILL on the CI runners (in gvec's library init), which took down the serial and the 2-rank unit test jobs. Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>

max-models · 2026-10-09T15:16:31Z

Remaining host work in the equilibria: what's needed to run everything on the GPU

On the CuPy backend, every analytic equilibrium is already built and evaluated fully on the device (covered by test_equils_cupy.py). Five equilibria still drop to the host, through xp.setup_on_host, xp.host_call and xp.evaluate_on_host:

Equilibrium	Host setup	Host evaluation, on every call
`AdhocTorus` (only `q_kind` 1 or 2)	`quad` in a Python loop, `UnivariateSpline`	`_psi_i(r, nu)` and `_p_i(r)`
`AdhocTorusQPsi`	`odeint`, `fsolve`, `UnivariateSpline`	`_psi_i(r, nu)`
`EQDSKequilibrium`	file reader, `RectBivariateSpline`, 3× `UnivariateSpline`, `minimize`	`psi(R,Z,dR,dZ)`, `q/g/p_psi`
`GVECequilibrium`	`gvec.State` (Fortran)	`bv`, `jv`, `p0`, `n0` all go through gvec
`DESCequilibrium`	DESC (JAX)	5 methods through DESC

The geometry setup also runs on the host: interp_mapping uses SciPy splu/spsolve, and field_line_tracing (used by the Tokamak domain) does a scalar Newton solve for each (s, ξ) point in a Python loop.

This matters at run time, not just at setup. pic/base.py calls equil.b_cart(markers) in boundary handling, and pic/particles.py calls equil.psi_r(r) on markers. For EQDSK, GVEC, DESC and the spline-based AdhocTorus variants, each of those calls is a device→host→device round trip.

1. Evaluation on the GPU (the important part)

Spline-based analytic equilibria (AdhocTorus, AdhocTorusQPsi, EQDSK). Keep the fit on the host, then convert the SciPy splines to device B-splines and evaluate those.

UnivariateSpline(s=0, ext=3) is a plain B-spline. Take t, c, k from spl._eval_args, rebuild it as a B-spline, clip the input to [t[k], t[-k-1]] to reproduce ext=3, and evaluate with nu=der.
RectBivariateSpline is a tensor-product B-spline (spl.tck, kx, ky) evaluated at scattered points with derivatives dx, dy.
After this, the only host work in these classes is setup.

GVEC and DESC. Evaluate B, J, p (and n) once at setup, interpolate them onto 3D tensor-product splines, and evaluate those on the device. The mapping is already handled this way, since GVECunit and DESCunit interpolate it with interp_mapping, so this extends the same approach to the fields.

The cost is interpolation error instead of the exact gvec/desc values. This could be an option, e.g. evaluate="spline" | "native".
DESC could also stay native through JAX on the GPU with DLPack zero-copy, but that is more fragile.

2. Setup on the GPU

field_line_tracing: worth doing. Replace the double loop of scalar newton calls with one batched Newton iteration over all (i, j) points at once. Once psi evaluates on the device (item 1), this runs on the GPU, and it will be much faster on the CPU too. The kron_lusolve_2d solve is small and can stay on the host or use cupyx.scipy.sparse.linalg.splu.
interp_mapping: could use cupyx.scipy.sparse.linalg (spsolve/splu). It's small and runs once, so the gain is low.
quad loops, odeint, fsolve, minimize (1D/scalar, about 50 points): leave these on the host. They are latency-bound and would be slower on a GPU, and cupyx has no integrate or optimize. The quad loop in AdhocTorus can become a spline antiderivative if we want it vectorised.
EQDSK file reading, gvec.State: necessarily host-only.

So "all setup on the GPU" isn't sensible, but everything that scales with grid or marker count, plus all evaluation, can move to the GPU.

Suggested cunumpy additions

Backend-agnostic B-spline evaluation: xp.interpolate.BSpline(t, c, k, extrapolate="clip") and NdBSpline / tensor-product with per-axis derivatives, evaluated at scattered points. It would wrap cupyx.scipy.interpolate where that covers it, with a pure-xp de Boor fallback so it also works under fake CuPy and Metal. To check: whether the minimum supported CuPy version has NdBSpline and 2D derivative evaluation.
xp.to_device_spline(scipy_spline): converts UnivariateSpline, RectBivariateSpline, BSpline or make_interp_spline output into (1), mapping ext to clip or extrapolate. This turns each host_call(self._psi_i, …) into a one-line change at setup.
Batched scalar root-finding: xp.optimize.newton(f, x0_array) with a vectorised secant/Newton step and a convergence mask, for field_line_tracing.
Fake CuPy support for cupyx.scipy.interpolate and cupyx.scipy.sparse.linalg (if not there yet), so test_equils_cupy.py can check items 1–3 without a GPU.

Suggested order

cunumpy (1) and (2), then switch AdhocTorus, AdhocTorusQPsi and EQDSK over. Their evaluation then runs fully on the device, and the existing parity tests check it.
Batched field_line_tracing with cunumpy (3).
Spline-projected evaluation for GVEC and DESC.

🤖 Generated with Claude Code

max-models · 2026-10-09T18:16:21Z

SciPy functions used in struphy that `cupyx.scipy` lacks

Checked struphy's SciPy imports (outside the tests) against CuPy's current reference pages.

cupyx.scipy has fft, fftpack, interpolate, linalg, ndimage, signal (+ windows), sparse (+ linalg, csgraph), spatial (+ distance), special and stats. These SciPy subpackages have no cupyx equivalent at all: integrate, optimize, cluster, io, odr, differentiate, datasets, constants, misc. The subpackages that do exist cover only part of SciPy.

What struphy uses that has no cupyx equivalent

Missing in cupyx	Where in struphy	Matters?
`scipy.linalg.solve_circulant`	`feec/preconditioner.py`, `linear_algebra/linalg_kron.py`	Yes. This is in the preconditioner, so it probably runs on every solver iteration. It's the only gap on a hot path.
`scipy.special.jv` (only `j0`, `j1` exist)	`initial/perturbations.py` (Bessel perturbations)	Runs once at initialisation; could use `host_call`.
`scipy.special.erfi`	`dispersion_relations/utilities.py`	Post-processing and analytic dispersion relations; stays on the host.
`RectBivariateSpline`	EQDSK	Already handled: `make_interp_spline` per axis + `NdBSpline` (matches `RectBivariateSpline(s=0)` to round-off; clip points first, since `NdBSpline` extrapolates).
`integrate.quad`, `odeint`	AdhocTorus, AdhocTorusQPsi setup	Scalar, run once, stay on the host.
`optimize.fsolve`, `minimize`	AdhocTorusQPsi, EQDSK setup; `dispersion_relations/analytic.py`	Scalar, stay on the host.
`optimize.newton`, `root`	`geometry/utilities.py` (`field_line_tracing`)	Covered by a new batched `xp.optimize.newton` in cunumpy.

What struphy uses that cupyx has

UnivariateSpline, sparse (csr_matrix, csc_matrix, kron, identity), sparse.linalg.splu and spsolve, fft, signal.argrelextrema, and special.erf, ndtri and yn.

Suggestion: `solve_circulant` in cunumpy

The only gap that matters for GPU runs is solve_circulant. It is easy to add to cunumpy: a circulant system diagonalises under the FFT, so on either backend it is

x = xp.fft.ifft(xp.fft.fft(b, axis=0) / xp.fft.fft(c)[:, None], axis=0).real

or the rfft version for real input. Batching over the right-hand-side columns, as linalg_kron.py needs, comes for free. It could live in cunumpy, e.g. as xp.scipy.linalg.solve_circulant filled in on the CuPy backend, so struphy doesn't need its own FFT helper in two places.

jv and erfi could get the same treatment, but they're off the hot path and xp.host_call is fine for them.

🤖 Generated with Claude Code

max-models · 2026-10-09T18:20:30Z

Plan: struphy changes to evaluate all equilibria on the GPU

This builds on max-models/cunumpy#93, which adds:

a fake cupyx.scipy, so xp.scipy code runs its CuPy path without a GPU, with only the interpolate names CuPy really has;
xp.optimize.newton, a batched Newton/secant solver.

It corrects my first comment above: we don't need our own B-spline evaluator. cupyx.scipy.interpolate has UnivariateSpline, make_interp_spline and NdBSpline, and xp.scipy.interpolate already picks SciPy or cupyx by backend.

Rule for all of the following: fit splines on the active backend with xp.scipy.interpolate.<name>, accessed at the point of use. Don't import them by name at module level, which would pin one backend. Only the scalar setup steps (quad, odeint, fsolve, minimize, file reading) stay on the host, under with xp.use_backend("numpy"):, and their results are copied to the device once with xp.to_cunumpy.

1. `AdhocTorus` (`q_kind` 1, 2)

Remove @xp.setup_on_host.
Run the two quad loops on the host.
Fit _psi_i and _p_i with xp.scipy.interpolate.UnivariateSpline(r_i, psis, k=psi_k, s=0.0, ext=3) on the active backend.
In psi_r and p_r, call the spline directly instead of xp.host_call(self._psi_i, r, nu=der).

2. `AdhocTorusQPsi`

Same as 1, with odeint and fsolve on the host and _psi_i fitted on the device.

3. `EQDSKequilibrium`

File reading and subsampling ([::smooth_step]) stay on the host.
Fit _g_i, _p_i and _q_i with xp.scipy.interpolate.UnivariateSpline(..., s=0.0, ext=3).

RectBivariateSpline is not in cupyx. With s=0 it is the same spline as

bx = xp.scipy.interpolate.make_interp_spline(R, psi, k=kx)
by = xp.scipy.interpolate.make_interp_spline(Z, bx.c.T, k=ky)
psi_i = xp.scipy.interpolate.NdBSpline((bx.t, by.t), by.c.T, (kx, ky))

It has the same knots, with values and derivatives (including (2,0), (0,2), (1,1)) equal to round-off.

RectBivariateSpline clamps points outside the grid, while NdBSpline extrapolates. So psi() must clip R and Z to range_R / range_Z before evaluating at xp.stack([R, Z], -1) with nu=(dR, dZ).
Run minimize (magnetic axis) against a host copy of the spline built the same way under use_backend("numpy"). Otherwise every optimizer step is a device evaluation plus a sync.

4. `field_line_tracing` (`Tokamak` domain)

The Newton guess for surface i comes from surface i-1 on the same ray (for i >= ns//2). So loop over the surfaces and solve all rays at once with xp.optimize.newton: ns batched solves instead of ns * nx scalar ones.
After that, drop the use_backend("numpy") around it in geometry/domains/tokamak/tokamak.py.
The small kron_lusolve_2d can stay on the host.

5. `GVECequilibrium` / `DESCequilibrium` (separate PR)

At setup, evaluate B, J, p and n natively on a grid, interpolate onto 3D tensor-product splines (make_interp_spline per axis + NdBSpline, or the existing spline_interpolation_nd), and evaluate those on the device.
Keep native evaluation as an option, e.g. evaluate="spline" | "native", since this trades exactness for speed.

6. Tests

test_equils_cupy.py already compares every equilibrium against NumPy on the fake CuPy. With cunumpy#93, the fake now also checks the xp.scipy.interpolate path and rejects names CuPy lacks.
Add xp.profiling.assert_no_transfers() around evaluation for 1–3, to prove no host round trips remain.
Needs a real GPU run: that CuPy's UnivariateSpline and NdBSpline match SciPy to the test's 1e-12. The fake runs SciPy itself, so it can't show this.
Possible catch on real CuPy: whether cupyx splines accept Python floats (self.psi(R0 + a, 0.0)). The fake accepts them. If CuPy doesn't, wrap scalars with xp.asarray.

Not part of the equilibria, but found along the way

scipy.linalg.solve_circulant (feec/preconditioner.py, linear_algebra/linalg_kron.py) has no cupyx equivalent and is on a hot path. Proposed as an FFT-based version in cunumpy (separate PR).
scipy.special.jv (initial/perturbations.py) and erfi (dispersion_relations/) are missing in cupyx, but run at initialisation or in post-processing only, so xp.host_call is fine.

Order: require a cunumpy release that includes #93, then 1 → 2 → 3 (with the test changes), then 4, then 5 separately.

🤖 Generated with Claude Code

spossann · 2026-10-09T19:24:02Z

Plan: struphy changes to evaluate all equilibria on the GPU

This builds on max-models/cunumpy#93, which adds:

a fake cupyx.scipy, so xp.scipy code runs its CuPy path without a GPU, with only the interpolate names CuPy really has;

xp.optimize.newton, a batched Newton/secant solver.

It corrects my first comment above: we don't need our own B-spline evaluator. cupyx.scipy.interpolate has UnivariateSpline, make_interp_spline and NdBSpline, and xp.scipy.interpolate already picks SciPy or cupyx by backend.

Rule for all of the following: fit splines on the active backend with xp.scipy.interpolate.<name>, accessed at the point of use. Don't import them by name at module level, which would pin one backend. Only the scalar setup steps (quad, odeint, fsolve, minimize, file reading) stay on the host, under with xp.use_backend("numpy"):, and their results are copied to the device once with xp.to_cunumpy.

1. AdhocTorus (q_kind 1, 2)

Remove @xp.setup_on_host.

Run the two quad loops on the host.

Fit _psi_i and _p_i with xp.scipy.interpolate.UnivariateSpline(r_i, psis, k=psi_k, s=0.0, ext=3) on the active backend.

In psi_r and p_r, call the spline directly instead of xp.host_call(self._psi_i, r, nu=der).

2. AdhocTorusQPsi

Same as 1, with odeint and fsolve on the host and _psi_i fitted on the device.

3. EQDSKequilibrium
File reading and subsampling ([::smooth_step]) stay on the host.

Fit _g_i, _p_i and _q_i with xp.scipy.interpolate.UnivariateSpline(..., s=0.0, ext=3).
RectBivariateSpline is not in cupyx. With s=0 it is the same spline as
bx = xp.scipy.interpolate.make_interp_spline(R, psi, k=kx)
by = xp.scipy.interpolate.make_interp_spline(Z, bx.c.T, k=ky)
psi_i = xp.scipy.interpolate.NdBSpline((bx.t, by.t), by.c.T, (kx, ky))
It has the same knots, with values and derivatives (including (2,0), (0,2), (1,1)) equal to round-off.
RectBivariateSpline clamps points outside the grid, while NdBSpline extrapolates. So psi() must clip R and Z to range_R / range_Z before evaluating at xp.stack([R, Z], -1) with nu=(dR, dZ).

Run minimize (magnetic axis) against a host copy of the spline built the same way under use_backend("numpy"). Otherwise every optimizer step is a device evaluation plus a sync.
4. field_line_tracing (Tokamak domain)

The Newton guess for surface i comes from surface i-1 on the same ray (for i >= ns//2). So loop over the surfaces and solve all rays at once with xp.optimize.newton: ns batched solves instead of ns * nx scalar ones.

After that, drop the use_backend("numpy") around it in geometry/domains/tokamak/tokamak.py.

The small kron_lusolve_2d can stay on the host.

5. GVECequilibrium / DESCequilibrium (separate PR)

At setup, evaluate B, J, p and n natively on a grid, interpolate onto 3D tensor-product splines (make_interp_spline per axis + NdBSpline, or the existing spline_interpolation_nd), and evaluate those on the device.

Keep native evaluation as an option, e.g. evaluate="spline" | "native", since this trades exactness for speed.

6. Tests

test_equils_cupy.py already compares every equilibrium against NumPy on the fake CuPy. With cunumpy#93, the fake now also checks the xp.scipy.interpolate path and rejects names CuPy lacks.

Add xp.profiling.assert_no_transfers() around evaluation for 1–3, to prove no host round trips remain.

Needs a real GPU run: that CuPy's UnivariateSpline and NdBSpline match SciPy to the test's 1e-12. The fake runs SciPy itself, so it can't show this.

Possible catch on real CuPy: whether cupyx splines accept Python floats (self.psi(R0 + a, 0.0)). The fake accepts them. If CuPy doesn't, wrap scalars with xp.asarray.

Not part of the equilibria, but found along the way

scipy.linalg.solve_circulant (feec/preconditioner.py, linear_algebra/linalg_kron.py) has no cupyx equivalent and is on a hot path. Proposed as an FFT-based version in cunumpy (separate PR).

scipy.special.jv (initial/perturbations.py) and erfi (dispersion_relations/) are missing in cupyx, but run at initialisation or in post-processing only, so xp.host_call is fine.

Order: require a cunumpy release that includes #93, then 1 → 2 → 3 (with the test changes), then 4, then 5 separately.

🤖 Generated with Claude Code

I agree with this strategy. The device-agnostic spline interpolation serms like a good idea.

Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>

…e tracing Implements items 1-4 and 6 of the plan in #714: - AdhocTorus (q_kind 1, 2), AdhocTorusQPsi, EQDSKequilibrium: the splines are fitted on the active backend with xp.scipy.interpolate and evaluated there, without host round trips. Only the scalar setup steps (quad, odeint, fsolve, the EQDSK file, minimize on a host copy of the flux spline) run on the host. RectBivariateSpline, which cupyx.scipy lacks, is replaced by make_interp_spline along R and Z as an NdBSpline, with the points clamped to the domain of the spline. - field_line_tracing: one batched xp.optimize.newton over all rays per flux surface instead of one scalar Newton per point; Tokamak runs its setup on the active backend. - test_equils_cupy: the evaluation of all equilibria except GVEC/DESC runs under xp.profiling.assert_no_transfers(). - cunumpy pinned to fb0f7c0a (#93: fake cupyx.scipy, xp.optimize.newton). Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>

Co-Authored-By: Claude Sonnet 5.5 <noreply@anthropic.com>

…ged) Co-Authored-By: Claude Sonnet 5.5 <noreply@anthropic.com>

…pe checking) Co-Authored-By: Claude Sonnet 5.5 <noreply@anthropic.com>

**Stack:** part 1 of 14, based on `devel`, next: #708. Full order: #709 → #708 → #711 → #712 → #705 → #713 → #714 → #718 → #720 → #721 → #724 → #725 → #723 → #722. --- Two small bugs that stop particle runs early on `CUNUMPY_BACKEND=cupy`: - **`SplineFunction.__init__`** (`feec/psydac_derham.py`): the knots of `derham.V0fem` are NumPy arrays, and `xp.ascontiguousarray` (= `cupy.ascontiguousarray`) does not accept them: `TypeError: Argument 'a' has incorrect type (expected cupy._core.core._ndarray_base, got numpy.ndarray)`. They are now made contiguous on the host and copied to the device with `xp.asarray`. - **`SortingBoxes._set_boxes`** (`pic/sorting.py`): `xp.sqrt` of a Python int returned a 0-d CuPy array, which `round()` rejects (`TypeError: type ndarray doesn't define __round__ method`). It is a host scalar, so this now uses `math.sqrt`. The NumPy backend is unaffected. ### Note With these fixes, `examples/VlasovAmpereOneSpecies/strong_Landau_damping` gets further on CuPy but still stops in `initialize_neighbours`: box sorting and the `vlasov_maxwell` accumulation are not on the GPU yet (see `CUDA_STRATEGY.md`, "Next"). This PR does not try to change that.

**Stack:** part 2 of 14, based on #709 (merge that first), next: #711. Full order: #709 → #708 → #711 → #712 → #705 → #713 → #714 → #718 → #720 → #721 → #724 → #725 → #723 → #722. --- ## Fuse the pusher preamble into one CUDA kernel On the CuPy backend, the three column-slice assignments at the start of `Pusher._push` ```python markers[:, first_pusher_idx:first_shift_idx] = markers[:, : 3 + vdim] markers[:, first_shift_idx:residual_idx] = 0.0 markers[:, residual_idx:-2] = 0.0 ``` each sweep the whole row-major marker buffer, because the slices are strided. With 4M markers (an 8M-row buffer of 1.66 GB), this took 4.6 ms per pusher call, more than the push kernels themselves. ## Changes - New src/struphy/pic/pushing/prepare_push_cuda.cu: one thread per entry of the contiguous column range [first_pusher_idx, n_cols - 2), so all three assignments are done in a single coalesced pass. - Pusher._push calls it on the CuPy backend; the kernel is compiled once, on first use. The NumPy path is unchanged. This roughly speeds up the preamble by 2x. The next big improvement for the markers would be #707 --------- Co-authored-by: Claude Opus 5.5 <noreply@anthropic.com>

…kernels) (#711) **Stack:** part 3 of 14, based on #708 (merge that first), next: #712. Full order: #709 → #708 → #711 → #712 → #705 → #713 → #714 → #718 → #720 → #721 → #724 → #725 → #723 → #722. --- **Solves the following issue(s):** Solves #704, tracked in #650. Until now a CUDA kernel compiled on its first call, which put the compile time into the first time step and its profiling regions. Now `Simulation.run()` calls the new `Simulation.compile_cuda_kernels()` after `allocate()`, so the kernels compile before the time stepping. **Core changes:** - **`Simulation.compile_cuda_kernels()`** (`simulation/sim.py`): - It runs after `allocate()` in `run()`, inside a new profiling region, `setup: compile cuda kernels`. - On the CuPy backend it calls `Kernel.compile()` once for each kernel of `Simulation.kernels()` and returns the compiled kernels. - On the NumPy backend it does nothing and returns `()`. NumPy runs are unaffected. - **Fail fast:** on CuPy, if any kernel the time stepping calls has no CUDA version, the method raises `NotImplementedError` before compiling anything. - The error lists every such kernel and the path where its `.cu` file is expected. For example, `VlasovAmpereOneSpecies` names `vlasov_maxwell`. - Without this check, the run fails at the first call of the first missing kernel, in the time loop. - Kernels with `missing_cuda="fallback"` are not counted as missing. - **`Simulation.kernels()`:** collects the kernels from the objects that keep and call them, through a new `kernels()` method on each owner. Each kernel is listed once (`utils/kernel_compilation.collect_kernels`). The owners: - `Domain`: the four geometry kernels. - `Derham`: the three spline evaluation kernels of `SplineFunction`, and feectools' `stencil_dot_3d`, `stencil_transpose_3d`, `stencil_inner_3d` and `stencil_axpy_3d` (called by the FEEC operators and solvers). - `Particles`: `reflect`, if an axis has reflecting boundary conditions. - `Pusher`: the pusher kernel and the kernels of its init/eval `KernelSetup`s. `KernelSetup`, `Accumulator` and `AccumulatorVector` list their own kernel. - `Propagator.kernels()`: by default it collects from the propagator's attributes. That covers `Kernel`s kept as attributes and objects with a `kernels()` method, also inside lists, tuples and dicts. It does not go into other objects. A propagator that calls a kernel it does not keep can override the method; none needs to today. - **Not listed:** kernels that only diagnostics or the initial solves call. They run during setup ("setup: initial diagnostics", "model post allocate") anyway. **Model-specific changes:** None. I checked the kernel lists against what the time loop actually calls. For each of 19 models (all particle models, plus Maxwell and LinearMHD), I ran the default parameter file for one step on NumPy while recording every `Kernel` call during `model.integrate`. Every kernel called in the time loop is in `Simulation.kernels()` for all 19 models. (This was an ad-hoc script; the committed test does the same check for `VlasovAmpereOneSpecies`.) **Documentation changes:** `CUDA_STRATEGY.md` has a short new section, "Compiling the CUDA kernels at setup (struphy#704)". The PR 16 notes now point to it from the sentence "when to compile them is decided later". **Tests** (new `simulation/tests/test_compile_cuda_kernels.py`): - `collect_kernels` walks owners and containers, lists each kernel once, and `missing_cuda` finds the kernels without a CUDA version. - For a small `VlasovAmpereOneSpecies`, `Simulation.kernels()` is exactly the domain + Derham kernels plus `push_eta_stage`, `vlasov_maxwell` and `push_v_with_efield`, and the only kernel missing CUDA is `vlasov_maxwell`. Reflecting particles add `reflect`. - Every `Kernel` called in the time loop of a one-step run is listed. - On NumPy, `compile_cuda_kernels()` is a no-op (`Kernel.compile` patched to fail). - With the backend switched to CuPy by monkeypatching: every kernel is compiled, and a missing CUDA version raises before anything is compiled. - The same CuPy path runs in a subprocess on cunumpy's fake CuPy (`CUNUMPY_FAKE_CUPY=1`, clean MPI environment), with `CudaKernel.compile` recorded. - GPU test (skipped here): `Vlasov` compiles all its kernels at setup, and `VlasovAmpereOneSpecies` raises naming `vlasov_maxwell`. **Testing** (macOS, no GPU; pyccel kernels compiled with GNU/Fortran; **GPU tests not run**): The suites were `simulation/tests`, `pic/tests/{test_kernel_setup,test_reflect,test_pushers,test_kernel_backends,test_cuda_emulation}.py`, `feec/tests/test_derham_gpu.py` and `geometry/tests/test_domain.py`. | | before (`devel`) | after | | --- | --- | --- | | those suites | 336 passed, 55 skipped | 343 passed, 56 skipped (7 new passed, the new GPU test skipped) | The "before" numbers are the same suites without the new test file; no existing test file was changed. Not verified without a GPU: the real NVRTC compile through `Kernel.compile()` and the GPU test. On the CPU, only the call path up to `CudaKernel.compile` is tested. 🤖 Generated with [Claude Code](https://claude.com/claude-code) --------- Co-authored-by: Claude Opus 5.5 <noreply@anthropic.com>

**Stack:** part 4 of 14, based on #711 (merge that first), next: #705. Full order: #709 → #708 → #711 → #712 → #705 → #713 → #714 → #718 → #720 → #721 → #724 → #725 → #723 → #722. --- **Solves the following issue(s):** Solves #698, part of #650. Markers move between ranks without going through the host on the CuPy backend: `Particles.mpi_sort_markers` exchanges the marker rows as device buffers with CUDA-aware MPI. **Core changes:** - **Device exchange of marker rows** (`Particles._sendrecv_markers_device` in `pic/base.py`, called by `_sendrecv_markers` on CuPy): - The rows to send (already gathered on the device per destination by `_sendrecv_get_destinations`) and one device receive buffer go to `Isend`/`Irecv` through `cunumpy.mpi.mpi_buffer`. With CUDA-aware MPI this passes the device arrays directly, after `cunumpy.mpi.synchronize_for_mpi`, as feectools' data exchangers have done since struphy-hub/feectools#86. Without CUDA-aware MPI the same buffers are staged through (pinned) host memory, with a one-time warning, instead of handing device pointers to a host-only MPI. - Rows from rank `i` arrive in a contiguous slice of the receive buffer. The holes are filled in the same order as on NumPy, but by one device scatter instead of one per source rank. - Messages with a count of zero are skipped. Both sides know every count after the `Alltoall`, so they skip the same messages. - Whether MPI is CUDA-aware: the answer already recorded by `xp.mpi.mpi_is_cuda_aware()`/`set_mpi_cuda_aware()`, otherwise a collective probe of `mpi_comm` at the first exchange. This is safe because `mpi_sort_markers` is collective. The probe copies 32 bytes to the host once. - **Counts on the host without copies:** `send_info`/`recv_info` are small NumPy arrays on every backend. The send counts are sizes of `nonzero` results, which are host integers already (CuPy synchronizes to size them). So the `Alltoall` of the counts passes host buffers, and the receive buffer, hole slices and load-imbalance check need no device reads. Before this PR, CuPy used device count arrays: one host write per rank, `Alltoall` on device buffers, and `list(recv_info)`, `recv_info[i]` and `first_hole[i]` read back element by element. - **`alpha`:** `_sendrecv_determine_mtbs` checks a tuple or list `alpha` on the host. Before, it asserted on device comparisons, which made two implicit reads. A uniform `alpha` (the usual case: 1.0 or 0.5 everywhere) becomes a Python scalar, so the 3-vector is no longer copied to the device. The values are the same on NumPy. - **NumPy path unchanged:** the same messages in the same order, and the same hole filling. The only differences are the NumPy count arrays, which `xp` already gave on NumPy, and the scalar `alpha`, which gives the same values. - **Tests:** new `pic/tests/test_sorting_device.py`: - `test_mpi_sort_markers_fake_cupy[bc, cuda_aware]` runs under `mpirun -n N env CUNUMPY_FAKE_CUPY=1 pytest --with-mpi ...` and is skipped unless the fake CuPy is active. It sorts the same markers on NumPy and on the fake CuPy across all ranks, two shifts in a row, so most markers change rank, with periodic and remove boundaries. It compares the whole marker arrays and `valid_mks` exactly, then sorts again with `do_test=True`. - With CUDA-aware MPI (the probe, or forced), it checks that `count_transfers()` is 0 during `mpi_sort_markers` and that no fake CuPy array is read by the host (`int`/`bool`/`float`/`__index__`/iteration/`.get`/`.item`/`.tolist` are patched to count). Before this PR there were 4 such reads per call. - With staging (`set_mpi_cuda_aware(False)`), it checks for exactly one `to_host` per non-empty send buffer, carrying `n_send * n_cols * 8` bytes, and one `to_device` when something is received. - Fake CuPy arrays are host memory, so the non-CUDA-aware Open MPI here takes them through `__cuda_array_interface__`. - `test_mpi_sort_markers_gpu[bc]` (`requires_cupy`, `mpirun -n 2` with one GPU per rank): the same check with the probed MPI, without the host-read patching. - CI: a step in `test-PR-unit-mpi.yml` runs the fake CuPy test with the same `n-procs` matrix (1–4). **Model-specific changes:** None. Not covered by this PR: - The SPH ghost-box exchange (`_communicate_boxes` → `_sendrecv_markers_boxes`) still passes device arrays to MPI without `synchronize_for_mpi`, and it reads `_recv_info_box` element by element. It can use the same pattern in a separate PR. - `apply_kinetic_bc`, which `mpi_sort_markers` calls first, was not changed. With the boundaries tested here, it makes no host/device copy and no host read. - A two-rank GPU run (the last item of #698) has not been done: no GPU here. **Documentation changes:** - `CUDA_STRATEGY.md`: the "MPI + GPUs" open question now refers to the new section "Marker exchange implementation notes (#698)". - Docstrings of `mpi_sort_markers`, `_sendrecv_all_to_all` and the new `_sendrecv_markers_device`/`_mpi_cuda_aware`. **Testing** (macOS, Open MPI 5.0.9 (not CUDA-aware), pyccel kernels compiled with GNU/Fortran; **GPU tests not run**): | | before (`devel`) | after | |---|---|---| | `mpirun -n 2 pytest --with-mpi` `test_sorting.py test_pushers.py test_draw_parallel.py test_neighbor_ranks.py test_accum_vec_H1.py` | 205 passed | 205 passed | | same, `mpirun -n 3` | 205 passed | 205 passed | | `mpirun -n 2 pytest --with-mpi test_sph.py test_kernel_setup.py test_tesselation.py` | – | 743 passed | | `mpirun -n 2 env CUNUMPY_FAKE_CUPY=1 pytest --with-mpi test_sorting_device.py` | 6 failed, 2 skipped | 6 passed, 2 skipped (GPU) | | same with `-n 1`, `-n 3`, `-n 4` | – | 6 passed, 2 skipped | On `devel`, the fake CuPy test fails because of the implicit host reads (`__bool__` ×3, `__iter__`) with CUDA-aware MPI, and because the device buffers are never staged when MPI is not CUDA-aware. 🤖 Generated with [Claude Code](https://claude.com/claude-code) --------- Co-authored-by: Claude Opus 5.5 <noreply@anthropic.com>

**Stack:** part 5 of 14, based on #712 (merge that first), next: #713. Full order: #709 → #708 → #711 → #712 → #705 → #713 → #714 → #718 → #720 → #721 → #724 → #725 → #723 → #722. --- **Solves the following issue(s):** Part of #688 (blocked matrix accumulations), tracked in #650. First matrix accumulation on the GPU: `linear_vlasov_ampere`, the accumulation of `EfieldWeightsCoupling` in `LinearVlasovAmpereOneSpecies` (step 2 of the porting order in `CUDA_STRATEGY.md`). It uses the `Array6D<double>` views of cunumpy 0.6.1. **Core changes:** - **CUDA kernel** `pic/accumulation/kernels/linear_vlasov_ampere/linear_vlasov_ampere_cuda.cu`: same arguments in the same order as the pyccel kernel. It launches one thread per marker row and skips holes and boundary particles. Per marker it calls `df`, `matrix_inv`, `matrix_vector`, `outer` and then `m_v_fill_b_v1_symm`. The six blocks `mat11 … mat33` are `Array6D<double>` (strided), the vectors are `Array3D<double>`, and `f0_values` is `const double*`, as in `push_weights_with_efield_lin_va`. - **Device helpers, 1:1 with pyccel** (same names, arguments and loop order): - `fill_mat` and `fill_mat_vec` in `filler_kernels.cuh` - `m_v_fill_b_v1_symm` in `particle_to_mat_kernels.cuh` - `outer` in `linalg_kernels.cuh` Every addition into matrix or vector data is a `cunumpy_atomic_add`, as in `charge_density_0form`. The spline values live in the thread-local `SplineScratch`. - **`Accumulator` matrix path on CuPy:** no code change was needed, and the new test checks it. On CuPy, `StencilMatrix._data`/`StencilVector._data` are device arrays. The accumulator zeroes them in place and passes the same arrays to the kernel. The ghost-region exchange, the transposed blocks of the symmetric matrix (feectools' device `transpose`), the scaled operator built by `EfieldWeightsCoupling` for `SchurSolver`, its `dot` and the Schur solve all run on them. - **Tests:** - `PARITY_CASES["linear_vlasov_ampere"]` covers 129 markers (with a hole and a boundary particle) in Cuboid, Colella, HollowTorus and a 3d spline mapping. The data are shaped like the real accumulation matrices (`v1_symm_accumulation_data()` in `kernel_test_args.py`). Tolerances are `rtol=1e-12`, `atol=1e-10`: atomics change the summation order and entries reach 1e4. The CPU emulation catches a transposed `DF`. - `m_v_fill_b_v1_symm` is compared with pyccel through a wrapper kernel: `test_v1_symm_filler` on a GPU, `test_emulated_v1_symm_filler` by CPU emulation. - New `test_accum_matrix_cupy.py` accumulates as `EfieldWeightsCoupling` does (Colella, clamped in eta1) and compares with NumPy: - It checks the nine kernel arrays, the dense blocks, the accumulated vector, `BC.dot(x)` and three CG iterations of the Schur solve. - It also checks that the kernel writes into the blocks' own `_data` (no copies), and that nothing is copied to the host from the accumulation up to `BC.dot(x)`. - On a GPU it runs with `assert_no_transfers`. Without a GPU it runs in a subprocess on cunumpy's fake CuPy (`CUNUMPY_FAKE_CUPY=1`, clean MPI environment). There, every `CudaKernel` launch is emulated on the CPU in place on the fake device arrays: the new `emulated_launches()` in `pic/tests/cuda_emulation.py` covers struphy's and feectools' kernels, struct arguments included. This test also fails on a transposed `DF`. **Model-specific changes:** None. `LinearVlasovAmpereOneSpecies` cannot run end to end on the GPU yet. Still needed: - **CG inner products:** feectools' `StencilVector.inner` returns its result on the host (`xp.to_numpy`), so each CG inner product copies a scalar. The Schur solve is outside the transfer guard because of this. - **Default preconditioner:** `MassMatrixPreconditioner` (the default of `EfieldWeightsCoupling`) cannot be created on CuPy, because `is_circulant` asserts a host array. The device Kronecker solve is a separate PR. - **Parallel feectools PRs:** the stencil 6D views and the device Kronecker solve. - **H100 run:** the first run on an H100 (#687), then the model end to end (#689). **Documentation changes:** `CUDA_STRATEGY.md`: - `linear_vlasov_ampere` is marked as ported in the checklist and the porting order. - 6D views are no longer a gate. - New section: implementation notes for `linear_vlasov_ampere`. **Fix added after review (commit 1979ea9): the test hole was not a hole for accumulation kernels.** `kernel_test_args.marker_arguments()` marked the hole (row 0) only in column 8 (`first_init_idx`), which the pushers test. Accumulation kernels, including `linear_vlasov_ampere`, test column 0, so this PR's parity case had no real hole. Row 0 is now `-1` in every column, as struphy marks holes. After the fix: `test_cuda_emulation.py` + `test_accum_matrix_cupy.py` 72 passed, 1 skipped; `test_cuda_parity.py -k linear_vlasov_ampere` passed on NumPy. **Testing** (macOS, no GPU; pyccel kernels compiled with GNU/Fortran; **GPU tests not run**): | | before (`devel`) | after | |---|---|---| | `test_kernel_backends.py test_cuda_parity.py test_device_helpers.py test_cuda_emulation.py` | 83 passed, 317 skipped | 85 passed, 322 skipped | | `test_accum_matrix_cupy.py` (new) | – | 1 passed (fake CuPy), 1 skipped (GPU) | | `mpirun -n 2`: `test_accum_vec_H1.py test_pushers.py test_verif_LinearVlasovAmpereOneSpecies.py` | 110 passed | 110 passed | The new skips are GPU tests: 4 parity cases, `test_v1_symm_filler` and `test_accumulator_matrix_on_cupy`. 🤖 Generated with [Claude Code](https://claude.com/claude-code) --------- Co-authored-by: Claude Opus 5.5 <noreply@anthropic.com>

**Stack:** part 6 of 14, based on #705 (merge that first), next: #714. Full order: #709 → #708 → #711 → #712 → #705 → #713 → #714 → #718 → #720 → #721 → #724 → #725 → #723 → #722. **Stack fix:** `simulation/tests/test_compile_cuda_kernels.py` (from #711) now expects `VlasovAmpereOneSpecies` to have no kernel without a CUDA version (`vlasov_maxwell` has one now). The fail-fast check uses `push_bxu_Hdiv` instead. The `CUDA_STRATEGY.md` example was updated to match. The `linear_vlasov_ampere` (#705) and `vlasov_maxwell` parity cases and imports are both kept. --- **Solves the following issue(s):** Part of #688 (blocked matrix accumulations), tracked in #650. Ports `vlasov_maxwell` to CUDA. It is the accumulation of `VlasovAmpereCoupling` in `VlasovAmpereOneSpecies` and `VlasovMaxwellOneSpecies` (#3 in the porting order of `CUDA_STRATEGY.md`). It uses the `Array6D<double>` views of cunumpy 0.6.1. This PR is based on `devel` and does not depend on #705 (`linear_vlasov_ampere`). The two PRs need the same device helpers. Here they are copied **byte-identically** from #705 (`origin/cuda-linear-vlasov-ampere` at 1ad4cd6), so `filler_kernels.cuh`, `particle_to_mat_kernels.cuh`, `kernel_test_args.py`, `test_device_helpers.py` and `test_cuda_emulation.py` merge without conflicts. Whichever PR merges second gets small conflicts in three files: `cuda_parity_cases.py`, `test_cuda_parity.py` and `CUDA_STRATEGY.md`. In all three, both sides add lines in the same place, so keep both (checked with `git merge-tree`). **Core changes:** - **CUDA kernel** `pic/accumulation/kernels/vlasov_maxwell/vlasov_maxwell_cuda.cu`: - Same arguments in the same order as the pyccel kernel, one thread per marker row. - Per marker it calls `df`, `matrix_inv`, `transpose`, `matrix_matrix` and `matrix_vector`. Then `m_v_fill_b_v1_symm` adds `A_p = w_p DF^{-1} DF^{-T}` and `B_p = w_p DF^{-1} v_p`. - The six blocks `mat11 … mat33` are `Array6D<double>` (any strides) and the vectors are `Array3D<double>`. - As in pyccel, only holes are skipped. Boundary particles (`markers[ip, -1] == -2`) are accumulated, unlike in `linear_vlasov_ampere`. - **Device helpers, 1:1 with pyccel, byte-identical to #705:** - `fill_mat` and `fill_mat_vec` in `filler_kernels.cuh` - `m_v_fill_b_v1_symm` in `particle_to_mat_kernels.cuh` Every addition into matrix or vector data is a `cunumpy_atomic_add`. `outer` from #705 is not needed here and is not copied. - **Tests:** - `PARITY_CASES["vlasov_maxwell"]` covers 129 markers in five mappings: Cuboid, Colella, HollowTorus, ShafranovDshapedCylinder and a 3d spline mapping (`geometry_domain` 0, 2, 5, 9, 12). - Row 0 is a real hole (`markers[0, 0] = -1`) and row 1 is a boundary particle. - The data are shaped like the real accumulation matrices (`v1_symm_accumulation_data()` from #705). - The CPU emulation agrees with pyccel to 7e-12 or better. Entries reach 6e5 in the spline mapping, because of `G^{-1}`. - Tolerances are `rtol=1e-12`, `atol=1e-8`, for the atomic summation order on the GPU. `atol` is scaled up from `linear_vlasov_ampere`'s 1e-10 because the entries are about 60 times larger. - Mutation checks: the emulated parity test fails for a transposed `DF`, for skipped boundary particles and for holes that are not skipped. - `vlasov_maxwell` is added to `test_vlasov_kernel_coverage`. - From #705 (identical): `test_v1_symm_filler` (GPU) and `test_emulated_v1_symm_filler` compare `m_v_fill_b_v1_symm` with pyccel. - **Not in this PR:** #705's `test_accum_matrix_cupy.py` and `emulated_launches()`. They test the matrix path of `Accumulator` on CuPy and are not copied, to keep this PR small. I ran an ad-hoc copy of that test with `vlasov_maxwell` in place of `linear_vlasov_ampere`, not committed. It runs on the fake CuPy with emulated launches and matches NumPy: the nine kernel arrays, the dense blocks, the vector, `BC.dot(x)` and three Schur CG iterations, with no host copies from the accumulation to `BC.dot(x)`. After both PRs are merged, that test could get a `vlasov_maxwell` case. **Model-specific changes:** None. Status of `VlasovAmpereOneSpecies` and `VlasovMaxwellOneSpecies` on the GPU: | Part | Kernels / operations | GPU status | |---|---|---| | `PushEta` | `push_eta_stage`, `reflect` | ported | | `PushVxB` (analytic and implicit) | `push_vxb_analytic`, `push_vxb_implicit` | ported | | `VlasovAmpereCoupling` | `vlasov_maxwell` | **this PR** | | | `push_v_with_efield` | ported | | | Schur solve | needs the FEEC parts below | | Initial `PoissonSolve` | `charge_density_0form` | ported | | | weights through the control variate (`update_weights`, `Maxwellian` on `xp`, `jacobian_det` through the geometry kernels) | ported | | | stiffness solve | needs the FEEC parts below | | `MaxwellWeakAmpere` (VlasovMaxwell only) | FEEC only (`curl`, `M1`, `M2`, Schur solve) | needs the FEEC parts below | | Scalars | `KineticEnergyPIC` (`xp.sum`) | runs on the device | | | `BilinearEnergyFEEC` | goes through `StencilVector.inner`, so a host scalar | | | `gauss_error` diagnostic | uses `toarray()`, a host copy | Still missing for an end-to-end GPU run: - **Preconditioner:** `MassMatrixPreconditioner`, the default of `VlasovAmpereCoupling`, `MaxwellWeakAmpere` and `PoissonSolve`, cannot be created on CuPy (`is_circulant` asserts a host array). The device Kronecker solve is a separate PR. - **CG inner products:** feectools' `StencilVector.inner` returns a host scalar, so each CG iteration and the FEEC energy scalars copy a scalar from the device. - **FEEC operators:** the derham operators and the mass matrices on the device need the feectools CUDA stack (stencil 6D views, device Kronecker solve, PR 17). - **Marker sorting:** multi-rank marker sorting still goes through the host. - **H100 run:** the first run on an H100 (#687), then the models end to end (#689). **Documentation changes:** `CUDA_STRATEGY.md`: - Porting order: #3 `vlasov_maxwell` is marked ✓. - New short section: `vlasov_maxwell` implementation notes. The checklist and the "6D views are a gate" paragraph are left to #705, to keep the conflict small. **Testing** (macOS, no GPU; pyccel kernels compiled with GNU/Fortran; **GPU tests not run**): | | before (`devel`) | after | |---|---|---| | `test_kernel_backends.py test_cuda_parity.py test_device_helpers.py test_cuda_emulation.py` | 83 passed, 317 skipped | 85 passed, 323 skipped | | `mpirun -n 2`: `test_verif_VlasovMaxwellOneSpecies.py test_verif_VlasovAmpereOneSpecies.py test_accum_vec_H1.py` | – | 7 passed | - The 2 new passes are `test_emulated_parity[vlasov_maxwell]` and `test_emulated_v1_symm_filler`. - The 6 new skips are GPU tests: 5 `vlasov_maxwell` parity cases and `test_v1_symm_filler`. - The MPI tests only use the pyccel kernels, which this PR does not change. I ran them only on this branch. 🤖 Generated with [Claude Code](https://claude.com/claude-code) --------- Co-authored-by: Claude Opus 5.5 <noreply@anthropic.com>

… generation) (#718) **Stack:** part 8 of 14, based on #714 (merge that first), next: #720. Full order: #709 → #708 → #711 → #712 → #705 → #713 → #714 → #718 → #720 → #721 → #724 → #725 → #723 → #722. --- **Solves the following issue(s):** Solves #690 (the decision part), tracked in #650. **Decision: every CUDA kernel is written by hand. Kernels are not generated from the pyccel source** (no numba, `cupyx.jit` or other code generation). To decide, one guiding-center kernel, `bstar_parallel_3form`, was ported by hand and through a numba.cuda prototype generator. The prototype has been removed from this PR now that the decision is made. The short version of the reasons is in `CUDA_STRATEGY.md`, section *Hand-written CUDA (decision for #690)*: - **The shared code is already ported.** Most of what a generator would save is the helper chain: mappings, metric, splines, linear algebra. PRs 10–19 have already ported and tested it by hand. What is left per kernel is its body: here 36 lines of CUDA plus one 5-line helper. - **A generator would be a second compiler to maintain, and its mistakes are silent.** It changes array lengths when a size is known only at run time, and it cannot tell which writes need atomics. The prototype translated only 31 of the 60 pusher and accumulation entry kernels. - **It does not fit the existing infrastructure.** Structs, `CudaKernel` signature checks, the CPU emulation of the real `.cu` and the parity tests are all built for CUDA C++. numba also adds about 175 MB of dependencies. - **Drift between pyccel and CUDA is caught by the parity cases.** **Core changes:** - **Hand-written port of `bstar_parallel_3form`:** - `pic/pushing/kernels/bstar_parallel_3form/bstar_parallel_3form_cuda.cu`, 1:1 with pyccel: same names and arguments, one thread per marker row, holes skipped as in pyccel (`markers[ip, 0] == -1`). `numpy.mod(eta, 1.0)` becomes `eta - floor(eta)` (Fortran `MODULO`, which pyccel generates). - It reuses `df` (every mapping), `det`, `get_spans`/`SplineScratch` and `eval_spline_mpi_kernel`. - New device helper `eval_0form_spline_mpi` in `bsplines/evaluation_kernels_3d.cuh`. - The pyccel kernel and its callers (`KernelSetup` in the guiding-center propagators) are unchanged. - **Parity cases:** 14 cases in `pic/tests/cuda_parity_cases.py`. - One per mapping: the 10 analytic mappings plus IGAPolarCylinder, IGAPolarTorus, Tokamak and a 3d spline. - α = 1, 0 and mixed. Shifts make `eta + shift` wrap around in `mod`, row 0 has a negative `eta_n` and row 2 is a hole. - Two output columns; launch size `args_markers.n_markers`. - The CPU emulation agrees with pyccel to at most 5.7e-14. Removing the `mod` from the CUDA kernel makes the emulation test fail. **Model-specific changes:** None. `bstar_parallel_3form` is one of the init kernels of `PushGuidingCenterBxEstar` and `PushGuidingCenterParallel`. Neither propagator runs on the GPU yet: the rest of their chain is not ported (`driftkinetic_hamiltonian`, `unit_b_1form`, `bstar_2form`, the discrete-gradient pushers, `gc_density_0form`). **Documentation changes:** `CUDA_STRATEGY.md`: - New section: *Hand-written CUDA (decision for #690)*. - The "Next" checklist entry, the step 3 note of the porting order and the open question on single-source alternatives point to the new section. **Testing** (macOS, no GPU; pyccel kernels compiled with GNU/Fortran; **GPU tests not run**, so the new kernel has not run on a GPU): | | before (`devel`) | after | |---|---|---| | `test_cuda_emulation.py`, `test_cuda_parity.py`, `test_kernel_backends.py`, `test_device_helpers.py` | 400 collected (83 passed, 317 skipped) | 415 collected: 84 passed, 331 skipped (+1 emulation test with 14 cases, +14 GPU parity tests skipped) | | `mpirun -n 1 pytest test_kernel_setup.py test_pushers.py` | 181 passed | 181 passed | The "before" pass/skip split is derived from the collected tests: this PR adds exactly one emulation test and 14 GPU parity tests. 🤖 Generated with [Claude Code](https://claude.com/claude-code) --------- Co-authored-by: Claude Opus 5.5 <noreply@anthropic.com>

**Stack:** part 10 of 14, based on #720 (merge that first), next: #724. Full order: #709 → #708 → #711 → #712 → #705 → #713 → #714 → #718 → #720 → #721 → #724 → #725 → #723 → #722. **Stack merge:** in `SplineFunction.__init__` (`feec/psydac_derham.py`), the `SplineArguments`/`CudaSplineArguments` are now built after #709's degree check and its host-to-device copy of the knots. --- **Solves the following issue(s):** Solves #678 (the remaining point: #667 (comment)). Part of #650. **Core changes:** `SplineFunction` no longer builds `(kind, pn, tn1, tn2, tn3, starts)` tuples for the spline evaluation kernels. It now uses an argument class pair, like `DerhamArguments`/`CudaDerhamArguments`: - New pyccel class `kernel_arguments/spline_args_kernels.SplineArguments(kind, pn, tn1, tn2, tn3, starts)` and its CUDA twin `kernel_arguments/spline_args_cuda.CudaSplineArguments` (a `CudaStructArguments` with struct `SplineArgs`). The header `kernel_arguments/spline_args.cuh` is generated by the new `write_spline_header()`. `SPLINE_STRUCTS` is added to `CUDA_STRUCTS`. - `SplineFunction.__init__` creates one argument object per component: `SplineArguments` on NumPy, `CudaSplineArguments` on CuPy. They are stored in `_args_spline`, which replaces `_args_eval`. The CUDA degree check (degrees 1 to 8) moves from `SplineFunction` into `CudaSplineArguments`, the same as in `CudaDerhamArguments`. - `eval_spline_mpi_markers`, `eval_spline_mpi_matrix` and `eval_spline_mpi_sparse_meshgrid` (pyccel and CUDA, still 1:1) take `args_spline` instead of the six loose arguments. The shared helper `eval_spline_mpi` (pyccel and `__device__`) keeps its flat signature, because other kernels call it too. - `CudaSplineArguments` accepts only C-contiguous knot vectors. `SplineFunction` already passes contiguous ones, so `test_evaluation_parity` now uses contiguous knots instead of strided views. No CUDA kernel other than these three takes the spline arguments, and pushers and accumulation already use `DerhamArguments`. One possible follow-up: `eval_spline_mpi_tensor_product_fixed` is a pyccel-only kernel in `evaluation_kernels_3d.py` that `SplineFunction` still calls with loose `kind, pn, starts` (no knots). Tests: - `test_kernel_backends.py`: the new pair is added to `ARGUMENT_PAIRS`, so the mirror check (and the struct layout check on a GPU) covers it. New test `test_generated_spline_header`. The owner test checks that `SplineFunction` creates `SplineArguments` on NumPy. - `kernel_test_args.spline_evaluation_arguments` returns `(_data, args_spline)`, so the parity and emulation cases now pass the argument object. `test_cuda_emulation.CUDA_CLASSES` includes `CudaSplineArguments`. Local tests (macOS, no GPU), after recompiling the changed kernels: - `bsplines/tests`, `feec/tests/test_eval_field.py`, `feec/tests/test_derham_gpu.py`, `pic/tests/test_cuda_emulation.py`, `pic/tests/test_kernel_backends.py`, `pic/tests/test_cuda_parity.py`, `pic/tests/test_pushers.py`, `pic/tests/test_accum_vec_H1.py`: 258 passed, 315 skipped. On devel, the same files gave 256 passed, 314 skipped; the differences are the 2 new tests and 1 new GPU-only skip. - `mpirun -n 2` on `bsplines/tests/test_eval_spline_mpi.py`, `feec/tests/test_eval_field.py`, `pic/tests/test_pushers.py`, `pic/tests/test_accum_vec_H1.py`: 139 passed (139 on devel). - CPU emulation of the three CUDA kernels against pyccel passes (`test_emulated_parity[eval_spline_mpi_*]`). - With cunumpy's fake CuPy, an ad-hoc check built `CudaSplineArguments` and packed the struct, and the degree check and the rejection of host arrays both work. The fake CuPy cannot construct a full `Derham`, so `SplineFunction` on CuPy was not run end-to-end. **GPU tests were not run** (no GPU available). This affects the struct layout, the parity tests and `test_spline_function_backends[cupy]`. **Model-specific changes:** None. **Documentation changes:** `CUDA_STRATEGY.md`: the new pair is added to the argument class table and the current-state table, plus a short note in the PR 16 implementation notes. 🤖 Generated with [Claude Code](https://claude.com/claude-code) --------- Co-authored-by: Claude Opus 5.5 <noreply@anthropic.com>

**Stack:** part 11 of 14, based on #721 (merge that first), next: #725. Full order: #709 → #708 → #711 → #712 → #705 → #713 → #714 → #718 → #720 → #721 → #724 → #725 → #723 → #722. **Stack fixes (for #721 below):** `SplineFunction.eval_tp_fixed_loc` now takes `kind`, `pn` and `starts` from `self._args_spline[i]`, because the tuple `_args_eval` no longer exists. The parity case of `eval_spline_mpi_tensor_product_fixed` unpacks `spline_evaluation_arguments` the same way. `test_catalog_signatures` now counts 4 spline evaluation kernels. Follow-up, not done here: switch `eval_spline_mpi_tensor_product_fixed` to the `SplineArguments` class of #721. --- **Solves the following issue(s):** Part of #697 (field evaluation and diagnostics on the device), tracked in #650. On the CuPy backend, the diagnostics of each output step now run on the device. Only the reduced results reach the host: the scalars in one copy per update, and every other saved dataset in one counted copy per output step. The marker diagnostics kernels in `pic/diagnostics/kernels` (used by the 5D and hybrid models) are not ported yet. They are listed as follow-up work, so this PR does not close #697. **Core changes:** I found the problems below without a GPU, on cunumpy's fake CuPy (`CUNUMPY_FAKE_CUPY=1`, clean environment) with `cunumpy.profiling.count_transfers`. CUDA launches were emulated on the CPU with `emulated_launches()` from #705's worktree, in ad-hoc probes that are not committed. In those probes, implicit syncs (`float()`, `bool()`, `__index__`, `.get()` of a device array) were counted by patching the fake `ndarray`. | Path | Before (fake CuPy) | After | |---|---|---| | `Scalars.update` (`FunctionScalarPIC` etc.) | one `float()` per function scalar (implicit sync) | no syncs; **one** `xp.to_numpy` of all values (8 bytes per scalar) | | `DataContainer.save_data` | one raw `.get()` per device dataset, not seen by `count_transfers`; scalars and time copied one by one | one counted `xp.to_numpy` per device dataset; scalars and time are already on the host, so they need no copy | | `Particles.binning` | **fails**: `xp.count_nonzero(list)` | runs `xp.histogramdd` on the device, no transfers | | `update_markers_to_be_saved` | `xp.count_nonzero` used as a slice bound (`__index__` sync) | the row count is taken from the gather; no extra sync | | `SplineFunction.__call__` (meshgrid/point) | six `bool()` syncs per call (process bounds compared as device scalars) | bounds cached on the host once; flagging uses `xp.copyto(..., where=...)` | | `SplineFunction.eval_tp_fixed_loc` | **fails**: pyccel kernel called with device arrays; `xp.array(degree)` uploaded per call | CUDA kernel, no transfers | | `Derham.prepare_eval_tp_fixed` (setup) | **fails** on CuPy (pyccel helpers given device buffers) | computed on the host, copied to the device once | - **Scalars** (`models/scalars.py`): `Scalars` keeps one device buffer. Each contained scalar's `value` is a view into it, so `update()` ends with `to_host()`, a single `xp.to_numpy`. `host_value(name)` returns a NumPy view of the host copy. `FunctionScalarFEEC/PIC/SPH` store the function's result through `_scalar_value` instead of `float()`, so a device reduction stays on the device. `StruphyModel.print_scalar_quantities` and the HDF5 setup read the host copy. - **Output** (`io/output_handling.py`): `DataContainer._as_numpy_array` is now `xp.to_numpy`. This is the one place where output copies device data. - **Time state** (`simulation/sim.py`): `time_state` is NumPy on both backends. The loop reads it with `float()` every step, which was a device sync on CuPy. `TimeDependentSource` still gets a scalar it can pass to `xp.sin`/`xp.cos`. - **Binning** (`pic/base.py`): the valid markers are gathered once, and the Jacobian is evaluated once instead of twice. On NumPy the results are the same values as before. - **Kernel** `eval_spline_mpi_tensor_product_fixed`: moved from `bsplines/evaluation_kernels_3d.py` into `bsplines/kernels/eval_spline_mpi_tensor_product_fixed/` (folder rule). It has a new `_cuda.cu` (one thread per grid point, reusing `eval_spline_mpi_kernel`) and parity cases in `pic/tests/cuda_parity_cases.py` (additive; #705/#713 also add entries there). `eval_tp_fixed_loc` takes `kind/pn/starts` from `SplineFunction._args_eval` instead of building them for every call. - **Not changed, still implicit syncs:** boolean-mask gathers (`markers[valid_mks]`, also inside `KineticEnergyPIC`), one per gather. `BilinearEnergyFEEC`/`VolumeFormEnergyFEEC` still copy one scalar per inner-product block inside feectools until #716 (struphy-hub/feectools#98) is in. With #716 these become device scalars, and this PR's container already handles them. **Model-specific changes:** None. `LinearVlasovAmpereOneSpecies._compute_en_w` now stays on the device because `FunctionScalarPIC` no longer calls `float()`. **Documentation changes:** `CUDA_STRATEGY.md`: one checklist entry and a short section "Diagnostics on the device (#697)" with the follow-ups: CUDA versions of the marker diagnostics kernels in `pic/diagnostics/kernels/`, and SPH kernel-density plots (step 6). **Tests:** - New `src/struphy/models/tests/test_diagnostics_on_device.py` has five checks (scalars, binning, saved markers, output, point flagging) with mock particles. They need no CUDA kernel. Each check runs: - on NumPy (`test_on_numpy`); - on fake CuPy in a subprocess (`test_on_fake_cupy`), compared with NumPy, counting transfers, with implicit syncs turned into errors; - on a real GPU (`test_on_cupy`, skipped here). `test_eval_tp_fixed_loc_on_cupy` (GPU, skipped here) compares the CUDA `eval_tp_fixed_loc` with NumPy. - New parity cases `eval_spline_mpi_tensor_product_fixed` (2 cases), run by CPU emulation in `test_cuda_emulation.py`. Local tests run (macOS, no GPU, pyccel kernels compiled with GNU/Fortran; **GPU tests not run**): - `models/tests/test_diagnostics_on_device.py`, `io/tests/test_output_handling.py`, `pic/tests/test_binning.py`, `models/tests/test_scalars.py`: 24 passed, 6 skipped (GPU), 2 failed. The two failures are `test_en_fB_ignores_ghost_markers` and `test_guiding_center_en_fB_does_not_divide_by_Np` in `test_scalars.py`. Their mocks have no `save_magnetic_background_energy`, which `GuidingCenter._compute_en_fB` calls, and the failures come from code this PR does not touch, so I think they also fail on `devel`. I did not confirm this with a separate `devel` run. - `pic/tests/test_cuda_emulation.py -k tensor_product` and the parity-case/signature checks of `test_cuda_parity.py`: passed (GPU parity skipped). - `models/tests/verification/test_verif_LinearVlasovAmpereOneSpecies.py`: 2 passed, serial and with `mpirun -n 2 ... --with-mpi`. - Not run locally: the full suite, the other model tests, and `test_variational_dissipation.py` (the models that use `eval_tp_fixed_loc`). CI covers these. 🤖 Generated with [Claude Code](https://claude.com/claude-code) --------- Co-authored-by: Claude Opus 5.5 <noreply@anthropic.com>

**Stack:** part 12 of 14, based on #724 (merge that first), next: #723. Full order: #709 → #708 → #711 → #712 → #705 → #713 → #714 → #718 → #720 → #721 → #724 → #725 → #723 → #722. --- **Solves the following issue(s):** Solves #411, tracked in #650. The FEEC assembly kernels are the last entry kernels that were still called as `PyccelKernel(getattr(module, name))` from shared modules. They now follow the folder rule of `CUDA_STRATEGY.md`, like the pusher, accumulation and geometry kernels. This is a pure restructuring: no numerical change. **Core changes:** - **One folder per kernel:** the 16 entry kernels of `feec/mass_kernels.py` and `feec/basis_projection_kernels.py` each move into `feec/kernels/<name>/`, which holds `<name>_kernels.py` and an `__init__.py` declaring `<name> = Kernel.from_folder(__name__, structs=CUDA_STRUCTS, missing_cuda="fallback")`. There is no `.cu` file, because none of these kernels has a CUDA port yet. Both old modules are removed. - From `mass_kernels.py`: `kernel_{1,2,3}d_mat`, `kernel_{1,2,3}d_vec`, `kernel_{1,2,3}d_eval`, `kernel_3d_matrixfree`, `kernel_3d_diag`, `surface_kernel_3d_mat`, `surface_kernel_3d_vec` - From `basis_projection_kernels.py`: `assemble_dofs_for_weighted_basisfuns_{1,2,3}d` - **The function bodies are unchanged.** I checked this with an AST comparison of all 16 functions against `devel`. Each file imports only what it uses (`numpy as np` and/or `shape`). - **Call sites** in `feec/mass.py`, `feec/boundary_mass.py` and `feec/basis_projection_ops.py` import the kernels. For kernels that depend on the dimension, they pick the right one from a small module-level dict (`_MAT_KERNELS[ldim]`, `_EVAL_KERNELS`, `_MATRIXFREE_KERNELS`, `_DIAG_KERNELS`, `_ASSEMBLY_KERNELS`) instead of `getattr` on a module. The direct pyccel calls of `mass_kernels.kernel_3d_vec` in `L2Projector` also go through the `Kernel` now. - **`missing_cuda="fallback"`:** on the CuPy backend these kernels keep the behaviour of the former `PyccelKernel` calls: they run on host copies of the arrays, with a warning the first time. With the default (`"raise"`), FEEC setups that currently assemble mass matrices on CuPy would start to fail. When a kernel is ported, its port removes the fallback for that kernel. - **`collect_kernel_files()`** picks up the new files with no change, since their names contain `kernels`. It now lists 142 files, 17 of them under `feec/kernels`. - **Tests:** `test_catalog_signatures` expects 17 kernels in `feec.kernels` (was 1). `feec.kernels` was already in `test_cuda_parity.PACKAGES`, so the new folders get the signature check and the folder-declaration check. **Model-specific changes:** None. **Documentation changes:** - `feec/kernels/__init__.py`: the docstring lists the assembly kernels and the fallback. - `CUDA_STRATEGY.md`: - a checklist entry for this PR - the kernel count in "Current state" - a short section "FEEC assembly kernels (#411)" saying that these kernels are now ready for CUDA ports (add `<name>_cuda.cu` and parity cases, then drop the fallback) **Testing** (macOS, no GPU; pyccel kernels compiled in the worktree with GNU/Fortran; **GPU tests not run**; following the agreed scope, I ran only a few directly relevant tests, and CI covers the rest): - `feec/tests/test_mass_reassembly.py`: 2 passed - `feec/tests/test_basis_ops.py::test_some_basis_ops` and `::test_transposed_update_weights_drops_zero_blocks`: passed - `feec/tests/test_boundary_integrals.py -k "unit_cube_constant or hcurl_per_face"`: 12 passed - `feec/tests/test_utilities_kernels.py`, `pic/tests/test_kernel_backends.py::test_catalog_signatures`, `pic/tests/test_cuda_parity.py::test_folders_declare_their_kernels`: passed - An ad-hoc script (not committed) on Colella with 6×5×4 elements and degree (2,3,2): - Matrix-free `M0`–`M3` `dot` agree with the assembled operators to about 1e-15 relative (`kernel_3d_matrixfree` against `kernel_3d_mat`). - Matrix-free and assembled `diagonal()` agree exactly (`kernel_3d_diag`). - `L2Projector.get_dofs` (`kernel_3d_vec`) and `eval_quad` (`kernel_3d_eval`) run. Not run locally: the full test suites, `test_mass_matrices.py` (32³ and 48³ grids, slow) and MPI runs. No MPI code changed. 🤖 Generated with [Claude Code](https://claude.com/claude-code) Co-authored-by: Claude Opus 5.5 <noreply@anthropic.com>

**Stack:** part 13 of 14, based on #725 (merge that first), next: #722. Full order: #709 → #708 → #711 → #712 → #705 → #713 → #714 → #718 → #720 → #721 → #724 → #725 → #723 → #722. **Stack fixes:** added `OUTPUTS` to the 16 FEEC assembly folders of #725, using the same arguments as the former `PyccelKernel(..., outputs=...)` wraps in `mass.py`/`boundary_mass.py`/`basis_projection_ops.py`, plus `kernel_*d_vec` and `surface_kernel_3d_vec`. Also added `OUTPUTS` for `eval_spline_mpi_tensor_product_fixed` (#724). The indices of `eval_spline_mpi_{markers,matrix,sparse_meshgrid}` were moved to the `SplineArguments` signatures of #721 (3, 5, 5). `argument_arrays` now knows `CudaSplineArguments`. The `test_cuda_emulation.py` merge keeps the tests added below (`test_emulated_v1_symm_filler`). Local re-run on the stacked branch: `test_cuda_parity.py`, `test_kernel_backends.py`, `test_evaluation_cuda.py` (59 passed, 336 skipped), `test_cuda_emulation.py` (74 passed). --- **Solves the following issue(s):** Solves #342, part of the CUDA work tracked in #650. Struphy's `PyccelKernel`s now declare which arguments they write to (`outputs=`), which cunumpy 0.6.1 already supports (`PyccelKernel(..., outputs=...)`, `Kernel.from_folder(..., host_options={"outputs": ...})`, and `assert_kernels_agree` uses `kernel.host_kernel.outputs` by default). Before this PR, no struphy kernel declared outputs. Without a declaration, the CuPy fallback copies every converted array back to the device after each call, and the parity tests compare every array. **Core changes:** - **Kernel folders (all 90):** each `<package>/kernels/<name>/__init__.py` declares ```python # the arguments the kernel writes to, by position: args_markers OUTPUTS = (2,) push_eta_stage = Kernel.from_folder(__name__, structs=CUDA_STRUCTS, host_options={"outputs": OUTPUTS}) ``` The declarations come from a static analysis of the pyccel sources, which follows writes through aliases (`markers = args_markers.markers`) and through called kernel helpers. I reviewed the results by hand. The breakdown: - 41 pushers write only `args_markers` (only its `markers` field). - The accumulation kernels write their `mat*`/`vec*` arrays. - The diagnostics write `markers` or `res`. - Evaluation, geometry, SPH and local-projector kernels write their `out`/`values`/`mat_f`/`f_eval*` argument. The scratch arrays `bn1 … bd3` of `DerhamArguments` are written by some kernels but are not outputs. They are not in the CUDA struct either. - **Other `PyccelKernel(...)` wraps** pass `outputs=` too: - `WeightedMassOperator` assembly (`kernel_{1,2,3}d_mat`, `-1`), `eval_quad` (`kernel_*d_eval`, `-1`), the matrix-free `dot` (`kernel_3d_matrixfree`, `-2`, which writes `out._data`; the last argument is `v._data`) and `diagonal` (`kernel_3d_diag`, `-1`), in `feec/mass.py` - `BoundaryMassOperator` (`surface_kernel_3d_mat`, `-1`), in `feec/boundary_mass.py` - `BasisProjectionOperator` (`assemble_dofs_for_weighted_basisfuns_*d`, `0`), in `feec/basis_projection_ops.py` Negative indices keep one declaration valid for every dimension. The test-only wraps in `test_mat_vec_filler.py` and `test_kernel_setup.py` are unchanged. - **Parity tests use the declarations:** - `test_cuda_parity.py`: `KernelCatalog.from_package` builds its own `Kernel`s and does not see the folders' `OUTPUTS`. The test therefore takes the kernels from the folder modules (`DECLARED`), and the GPU `test_parity` (`assert_kernels_agree`) now compares only the declared outputs. A side effect, found by reading the code but not run on a GPU: before, `assert_kernels_agree` also collected the host `DerhamArguments` scratch arrays (`argument i.bn1`, …). The CUDA struct has no such fields, so the host and CUDA key sets could not match for kernels taking `args_derham`. - New `test_kernels_declare_outputs` (no GPU needed): every kernel folder declares non-empty, in-range `OUTPUTS`. - New `test_declared_outputs` (NumPy backend, every parity case): the pyccel kernel changes at least one declared output and no other argument array. Argument objects are compared through their CUDA struct fields, so scratch is ignored. The helper `argument_arrays` lives in `cuda_parity_cases.py`. - `test_cuda_emulation.test_emulated_parity` now compares only the outputs with the pyccel kernel, and also asserts that the emulated CUDA kernel leaves every non-output array unchanged. I checked both new checks by mutation: declaring `push_eta_stage` with `OUTPUTS = (4,)` makes both tests fail. - `CUDA_STRATEGY.md`: short note on `OUTPUTS`. **Model-specific changes:** None. Behaviour on the NumPy backend is unchanged, because `PyccelKernel` uses `outputs` only when it converts CuPy arrays. **Documentation changes:** A short note in `CUDA_STRATEGY.md`. **Local tests run** (macOS, no GPU; GitHub CI covers the rest): - `pic/tests/test_cuda_parity.py`: 18 passed, 213 skipped (the GPU tests) - `pic/tests/test_cuda_emulation.py -k emulated_parity`: all 15 kernels passed. This ran as two invocations: the first hit my 5-minute limit after 13 kernels, and the remaining two (`kernel_pullpush_pic`, `kernel_pullpush`) passed separately in 80 s. - `pic/tests/test_kernel_backends.py` and `test_kernel_setup.py`: 89 passed, 22 skipped - `feec/tests/test_boundary_integrals.py` and `test_mass_reassembly.py`: 25 passed GPU tests were not run. On the GPU, `test_parity` now compares only the declared outputs. The `outputs=` of the FEEC mass/projection wraps only take effect on the CuPy fallback path. No local test exercises them, neither on a GPU nor with the fake CuPy. 🤖 Generated with [Claude Code](https://claude.com/claude-code) Co-authored-by: Claude Opus 5.5 <noreply@anthropic.com>

**Stack:** part 14 of 14, based on #723 (merge that first), top of the stack. Full order: #709 → #708 → #711 → #712 → #705 → #713 → #714 → #718 → #720 → #721 → #724 → #725 → #723 → #722. **Stack fixes:** applied the `struphy_cuda::<pyccel module>` rule to the CUDA code of the PRs below: `linear_vlasov_ampere` (#705), `vlasov_maxwell` (#713), `bstar_parallel_3form` (#718), `eval_spline_mpi_tensor_product_fixed` (#724) and the `SplineArgs` versions of `eval_spline_mpi_*` (#721). Each of them has `using namespace struphy_cuda;` and qualified helper calls. The helpers `outer`, `fill_mat`, `fill_mat_vec` and `m_v_fill_b_v1_symm` (#705/#713) are inside their module namespaces, and `m_v_fill_b_v1_symm` calls `filler_kernels::`/`evaluation_kernels_3d::`. The `fill_v1_symm` test wrapper is qualified. The "open PRs that need the same treatment" below are therefore covered by this stack. Local re-run on the stacked branch: `test_cuda_parity.py`, `test_kernel_backends.py`, `test_evaluation_cuda.py`, `test_device_helpers.py` (59 passed, 419 skipped), `test_cuda_emulation.py` (74 passed; it compiles every `.cu`/`.cuh` with the host C++ compiler). --- **Solves the following issue(s):** Closes #672, tracked in #650. The CUDA device helpers now live in C++ namespaces that mirror the pyccel modules they port. A pyccel call `linalg_kernels.matrix_inv(...)` becomes `linalg_kernels::matrix_inv(...)` in CUDA, so the 1:1 correspondence is visible at every call site and two helpers with the same name in different modules (`det`/`det_df`, `f`/`df`, `push`/`pull`) can no longer clash. **Core changes:** - **One namespace per pyccel module**, nested in `struphy_cuda` as proposed in the issue (`struphy_cuda::linalg_kernels::matrix_inv`). The namespace is named after the module the header ports: - `bsplines_kernels`, `evaluation_kernels_2d`, `evaluation_kernels_3d` (`bsplines/`) - `evaluation_kernels`, `spline_mappings_kernels`, `transform_kernels` (`geometry/`) - one per mapping, named after its pyccel module: `cuboid_kernels`, `orthogonal_kernels`, `colella_kernels`, `hollow_cylinder_kernels`, `powered_elliptic_cylinder_kernels`, `hollow_torus_kernels`, `shafranov_shift_cylinder_kernels`, `shafranov_sqrt_cylinder_kernels`, `shafranov_dshaped_cylinder_kernels` (`geometry/domains/<name>/<name>_cuda.cuh`) - `linalg_kernels`, `filler_kernels`, `particle_to_mat_kernels`, `pusher_utilities_kernels` - **CUDA-only helpers and constants** go into the namespace of the module whose calculation they extract: `bsplines_kernels::MAX_SPLINE_DEGREE`, `evaluation_kernels_3d::SplineScratch`, `spline_mappings_kernels::spline_index_row`/`first_plane`, `transform_kernels::column_norm`, `pusher_utilities_kernels::reflect_velocity`. - **Call sites:** - Inside a header, calls into another module are qualified (`linalg_kernels::det(...)` in `evaluation_kernels::det_df`). Calls within the same module are not, as in pyccel. - Every `.cu` file writes `using namespace struphy_cuda;` after its includes and qualifies every helper call (`evaluation_kernels_3d::get_spans(...)`, `evaluation_kernels::df(...)`). - The CUDA sources embedded in `test_device_helpers.py` and `geometry/tests/spline_mapping_cases.py` use fully qualified names. - **Left global, on purpose:** - `struphy_cuda::pi` (`geometry/domains/constants_cuda.cuh`). It stands in for `numpy.pi`, not for a struphy module. - The generated argument structs `MarkerArgs`, `DerhamArgs`, `DomainArgs` and `LocalProjectorsArgs` (`kernel_arguments/*.cuh`, written by cunumpy's `write_cuda_header`). cunumpy matches the `extern "C" __global__` signatures against the struct names, and the CPU emulation (`pic/tests/cuda_emulation.py`) rebuilds the structs by name, so the generators and the emulation are unchanged. - **`evaluation_kernels_3d.cuh`:** it had three `namespace struphy_cuda` blocks with an `#include` in between. It now has one block. `eval_spline_mpi` was outside the namespace and is now `evaluation_kernels_3d::eval_spline_mpi`, like its pyccel counterpart. - Pure refactor: no function body, signature or entry kernel changes. Lines changed by the renaming were rewrapped at 120 columns. **Open PRs that need the same treatment after merge:** #705, #713 and #718 add `.cu`/`.cuh` code that still calls `struphy_cuda::<name>`. Whichever of them and this PR merges second needs the rule applied: `using namespace struphy_cuda;` in new `.cu` files, and new helpers in their module's namespace. Their branches were not touched. **Model-specific changes:** None. **Documentation changes:** `CUDA_STRATEGY.md`: a new convention bullet under "Target layout" documents the rule: the namespace is `struphy_cuda::<pyccel module>`, CUDA-only helpers go into their module, cross-module calls are qualified, `.cu` files use `using namespace struphy_cuda;`, and the exceptions are `pi` and the generated structs. **Local tests** (macOS, no GPU; pyccel kernels compiled with GNU/Fortran; **GPU tests not run**, so nothing was compiled with NVRTC). The CPU emulation compiles every `.cu` file and every helper header as C++17 with the host compiler. That is the syntax check for the namespaces; it caught one misplaced brace while I was working on this. - `src/struphy/pic/tests/test_cuda_parity.py` - `src/struphy/pic/tests/test_cuda_emulation.py` - `src/struphy/pic/tests/test_device_helpers.py` - `src/struphy/bsplines/tests/test_evaluation_cuda.py` - `src/struphy/pic/tests/test_kernel_backends.py` Result: 107 passed, 392 skipped (the GPU-only tests). `devel` gives the same counts. One open point for the GPU run: the nested namespace definition (`namespace struphy_cuda::bsplines_kernels {`) needs C++17. NVRTC uses C++17 by default since CUDA 11, and the emulation compiles with `-std=c++17`. 🤖 Generated with [Claude Code](https://claude.com/claude-code) Co-authored-by: Claude Opus 5.5 <noreply@anthropic.com>

max-models and others added 2 commits October 7, 2026 21:11

max-models changed the base branch from devel to cuda-vlasov-maxwell October 7, 2026 22:06

max-models added this pull request to stack #728 October 8, 2026 05:45

max-models removed this pull request from stack #728 October 8, 2026 06:00

max-models added this pull request to stack #729 October 8, 2026 06:00

Merge cuda-vlasov-maxwell (#713) into mhd-equilibria-cupy

2e92d25

max-models marked this pull request as draft October 8, 2026 09:14

max-models and others added 5 commits October 8, 2026 11:46

Merge branch 'cuda-vlasov-maxwell' into mhd-equilibria-cupy

2e30c40

Merge parent into mhd-equilibria-cupy

d269195

Skip the new GVEC tests by default (with_gvec=False), like the other …

8223abf

…GVEC tests GVEC dies with SIGILL on the CI runners (in gvec's library init), which took down the serial and the 2-rank unit test jobs. Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>

Merge branch 'cuda-vlasov-maxwell' into mhd-equilibria-cupy

39e0e8d

max-models and others added 2 commits October 10, 2026 07:59

Merge cuda-vlasov-maxwell into mhd-equilibria-cupy

dcedf86

Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>

max-models mentioned this pull request Oct 10, 2026

Fake cupyx.scipy: copy and deepcopy of proxy objects max-models/cunumpy#97

Merged

max-models and others added 3 commits October 10, 2026 16:11

Pin cunumpy to the commit with copyable fake cupyx splines (cunumpy#97)

70557fb

Co-Authored-By: Claude Sonnet 5.5 <noreply@anthropic.com>

Pin cunumpy to devel with copyable fake cupyx splines (cunumpy#97 mer…

822af64

…ged) Co-Authored-By: Claude Sonnet 5.5 <noreply@anthropic.com>

Pin cunumpy to devel (f996393: copyable fake cupyx splines, faster ty…

dce63a9

…pe checking) Co-Authored-By: Claude Sonnet 5.5 <noreply@anthropic.com>

max-models marked this pull request as ready for review October 11, 2026 22:00

Merge branch 'cuda-vlasov-maxwell' into mhd-equilibria-cupy

6ca659b

spossann reviewed Oct 12, 2026

View reviewed changes

Comment thread src/struphy/fields_background/equils.py

Comment thread src/struphy/fields_background/equils.py

spossann approved these changes Oct 12, 2026

View reviewed changes

max-models merged commit 952222a into devel Oct 12, 2026
28 checks passed

max-models deleted the mhd-equilibria-cupy branch October 12, 2026 10:03

max-models mentioned this pull request Oct 14, 2026

MHD equilibria on the CuPy backend #696

Closed

3 tasks

Navigation Menu

Repository navigation

MHD equilibria on the CuPy backend - #714

MHD equilibria on the CuPy backend#714
max-models merged 14 commits into
cuda-vlasov-maxwellfrom
mhd-equilibria-cupy

max-models commented Oct 7, 2026 •

edited

Loading

Uh oh!

max-models commented Oct 9, 2026

Uh oh!

max-models commented Oct 9, 2026

Uh oh!

max-models commented Oct 9, 2026

Uh oh!

spossann commented Oct 9, 2026

Plan: struphy changes to evaluate all equilibria on the GPU

1. `AdhocTorus` (`q_kind` 1, 2)

2. `AdhocTorusQPsi`

3. `EQDSKequilibrium`

4. `field_line_tracing` (`Tokamak` domain)

5. `GVECequilibrium` / `DESCequilibrium` (separate PR)

6. Tests

Not part of the equilibria, but found along the way

Uh oh!

Uh oh!

Uh oh!

Uh oh!

Reviewers

Assignees

Labels

Projects

Milestone

Development

Uh oh!

2 participants

Conversation

max-models commented Oct 7, 2026 • edited Loading Uh oh! There was an error while loading. Please reload this page.

Uh oh!

Uh oh!

max-models commented Oct 9, 2026

Remaining host work in the equilibria: what's needed to run everything on the GPU

1. Evaluation on the GPU (the important part)

2. Setup on the GPU

Suggested cunumpy additions

Suggested order

Uh oh!

max-models commented Oct 9, 2026

SciPy functions used in struphy that cupyx.scipy lacks

What struphy uses that has no cupyx equivalent

What struphy uses that cupyx has

Suggestion: solve_circulant in cunumpy

Uh oh!

max-models commented Oct 9, 2026

Plan: struphy changes to evaluate all equilibria on the GPU

1. AdhocTorus (q_kind 1, 2)

2. AdhocTorusQPsi

3. EQDSKequilibrium

4. field_line_tracing (Tokamak domain)

5. GVECequilibrium / DESCequilibrium (separate PR)

6. Tests

Not part of the equilibria, but found along the way

Uh oh!

spossann commented Oct 9, 2026

Plan: struphy changes to evaluate all equilibria on the GPU

1. AdhocTorus (q_kind 1, 2)

2. AdhocTorusQPsi

3. EQDSKequilibrium

4. field_line_tracing (Tokamak domain)

5. GVECequilibrium / DESCequilibrium (separate PR)

6. Tests

Not part of the equilibria, but found along the way

Uh oh!

Uh oh!

Uh oh!

Uh oh!

Reviewers

Assignees

Labels

Projects

Milestone

Development

Uh oh!

2 participants

max-models commented Oct 7, 2026 •

edited

Loading

SciPy functions used in struphy that `cupyx.scipy` lacks

Suggestion: `solve_circulant` in cunumpy

1. `AdhocTorus` (`q_kind` 1, 2)

2. `AdhocTorusQPsi`

3. `EQDSKequilibrium`

4. `field_line_tracing` (`Tokamak` domain)

5. `GVECequilibrium` / `DESCequilibrium` (separate PR)

1. `AdhocTorus` (`q_kind` 1, 2)

2. `AdhocTorusQPsi`

3. `EQDSKequilibrium`

4. `field_line_tracing` (`Tokamak` domain)

5. `GVECequilibrium` / `DESCequilibrium` (separate PR)