Polar splines on the CuPy backend

[max-models]

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

---

Solves the following issue(s):

Solves #695, part of #650. Derham(polar_splines=True) raised NotImplementedError on CuPy (since PR 19, #683). The polar extraction operators applied SciPy sparse matrices to stencil data that lives on the device. The first failure was csr_matrix of a device array in PolarExtractionBlocksC1.

Core changes:

• Blocks on the host (polar/extraction_operators.py): PolarExtractionBlocksC1 builds every block on every backend with NumPy/SciPy, from a host copy of the control points (xp.to_numpy(domain.cx)). This covers the basis and DOF extraction blocks and the polar grad/curl/div blocks. The blocks are setup data, so .T, .toarray() and the restart left inverse keep working on them. The legacy PolarSplines_C0_2D/C1_2D classes are unchanged.

• Products on the device (polar/linear_operators.py):

  • PolarExtractionOperator and PolarLinearOperator keep their blocks_* as host SciPy matrices.
  • On CuPy, dot uses device copies: DeviceSparseMatrix, a cupyx.scipy.sparse.csr_matrix obtained through cunumpy's xp.scipy.sparse.
  • The copies are made once per block list, on first use, and the block setters drop them.
  • The blocks are small: they act only on the first 2–3 radial rings (at most 3·n2 × 3·n2 in eta1–eta2, n3 × n3 in eta3). Blocks with no non-zeros are not copied; their products are zeros.
  • Inputs are made contiguous before the cuSPARSE product, because kron_matvec_2d passes transposed views.
  • The NumPy path is unchanged: on NumPy, _backend_blocks returns the SciPy blocks themselves.
  • set_device_sparse_module(module) swaps the sparse backend. This is only for tests without a GPU.

• MPI: the ring reductions in dot_inner_tp_rings, dot_parts_of_polar and PolarVector.toarray call cunumpy.mpi.synchronize_for_mpi before Allreduce. On CuPy these are device buffers (CUDA-aware MPI).

• Derham: the NotImplementedError guard for polar splines on CuPy is removed. SplineFunction._restart_extraction_op builds its pseudo-inverse with NumPy (host blocks).

• Tests: new polar/tests/test_polar_cupy.py. It builds the same polar Derham on NumPy and CuPy for IGAPolarCylinder and IGAPolarTorus and compares:

  • E, P (DOF extraction) and their transposes for all five spaces
  • grad/curl/div and their transposes
  • PolarVector arithmetic (+ - * neg += -= *= dot copy toarray)
  • the polar projectors P0..P3
  • polar SplineFunction coefficients (extract_coeffs, restart inverse)
  • the polar mass matrices M0..Mv

  It also checks that every result lives on the expected backend.

  • Without a GPU (test_polar_fake_cupy, subprocess with CUNUMPY_FAKE_CUPY=1): cupyx.scipy.sparse is replaced by a dense device stand-in, and kernels run in their host version, because the fake cannot launch CUDA kernels. The mass matrices are left out here: their weights need struphy's CUDA geometry kernels, which take CUDA argument objects. They could be added once emulated_launches() from CUDA version of the linear_vlasov_ampere accumulation #705 is on devel. A mutation check confirmed that the test catches host blocks applied to device data.
  • On a GPU (test_polar_on_cupy, skipped here): the full comparison, mass matrices included, with real cupyx.scipy.sparse.
  • test_device_blocks_numpy_backend checks that the NumPy path uses the SciPy blocks unchanged.
  • The obsolete test_polar_splines_not_supported_on_cupy is removed.

Not covered: SplineFunction.__call__ (field evaluation is not ported to CuPy for any space yet), and the polar mass-matrix preconditioners on CuPy.

Model-specific changes:

None.

Documentation changes:

CUDA_STRATEGY.md:

• Polar splines are no longer listed as host-only (current state, gates, open questions, PR 19 notes).
• New short section "Polar splines on CuPy (Polar splines on the CuPy backend #695)".

Local testing (macOS, no GPU, pyccel kernels compiled with GNU/Fortran; GPU tests not run; GitHub CI covers the rest):

• polar/tests/test_polar_cupy.py, feec/tests/test_derham_gpu.py, test_l2_projectors.py::test_l2_projectors_polar: 6 passed, 5 skipped (GPU)
• polar/tests/test_polar.py (NumPy): 13 passed. test_mass_matrices.py::test_mass_polar (first two cases): 2 passed.
• mpirun -n 2 on test_polar.py::test_polar_adjoints_small_nel2, test_extraction_ops_and_derivatives and test_restart_polar: 9 passed.
• I stopped test_mass_matrices.py::test_mass_preconditioner_polar locally: one case takes about 22 min on this machine (it passed on unchanged devel). I did not run test_basis_ops_polar after the change.

🤖 Generated with Claude Code