MassMatrixPreconditioner on the CuPy backend

[max-models]

Solves the following issue(s):

Part of #689 (LinearVlasovAmpereOneSpecies on the GPU), tracked in #650. EfieldWeightsCoupling uses MassMatrixPreconditioner by default, and the preconditioner could not even be created on the CuPy backend.

The failure (reproduced with cunumpy's fake CuPy, CUDA launches emulated on the CPU):

File "src/struphy/feec/preconditioner.py", line 234, in __init__
    if is_circulant(M_arr):
File "src/struphy/feec/preconditioner.py", line 969, in is_circulant
    assert isinstance(mat, xp.ndarray)
AssertionError

StencilMatrix.toarray() returns a host (NumPy) array on every backend, because feectools builds it through a scipy COO matrix. is_circulant and FFTSolver.__init__ asserted xp.ndarray, which is cupy.ndarray on CuPy. The next steps would have failed as well: xp.nonzero(M_arr) on a host array, and per-entry writes of host scalars into the device _data of the local stencil matrix.

Core changes:

• The 1d setup data stay on the host, on every backend. These are the dense 1d mass matrices and the 1d solvers (FFTSolver for circulant matrices, SparseSolver otherwise). KroneckerLinearSolver expects host 1d solvers. feectools#96 builds its dense device inverses from these host solvers (M = solver.solve(I) on host arrays), so this matches the device Kronecker solve as well.
• Only the process-local 1d stencil matrices go to the device. They are the factors of the KroneckerStencilMatrix. Each one is filled in a host array and copied to M_local._data in one step.
• New helper _solver_and_local_matrix_1d. The 1d setup (dense matrix, solver, local stencil matrix, consistency check) was duplicated in MassMatrixPreconditioner and MassMatrixDiagonalPreconditioner. Both now call this helper, which uses NumPy explicitly. Both classes had the same problem, and both are fixed.
• is_circulant accepts host or device matrices and checks on the host. It returns a Python bool.
• FFTSolver keeps its circulant column on the host and accepts a host or device matrix. solve with device right-hand sides (or out) solves on the host and copies the result back. This is what BandedSolver/SparseSolver do in the current feectools, and it stays until the device Kronecker solve of feectools#96 is in. With Only run the unit tests which touch changed code in the ci #96, FFTSolver is only called on host arrays (to build the dense inverse), so it needs no change there. except xp.linalg.LinAlgError is now np.linalg.LinAlgError, which is what scipy raises.
• NumPy path: same arithmetic as before. The one structural difference is that M_local is no longer copied a second time (M_local.copy()), because it is a fresh object.
• The other classes in the module (MassMatrixDiagonalPreconditioner, FFTSolver) are covered above. The module has no other preconditioners.

Tests: new feec/tests/test_preconditioner_cupy.py.

• check_preconditioners_on_cupy builds MassMatrixPreconditioner and MassMatrixDiagonalPreconditioner for M0, M1 and M2 (Colella, 6x5x4 elements, degrees (2, 3, 1)) on both backends. It covers periodic bcs and clamped bcs (Dirichlet in eta1/eta3, so the boundary-operator path and both FFTSolver and SparseSolver are used).
• It applies them to the same random vector, dot and dot(out=...), checks that the results live on the device, and compares with NumPy at rtol=1e-12.
• Without a GPU, the check runs in a subprocess on cunumpy's fake CuPy (clean serial environment). There, every CudaKernel launch (feectools' stencil_transpose_1d, struphy's kernel_evaluate) is emulated on the CPU through a small private _emulated_launches(). It is a minimal version of emulated_launches() from CUDA version of the linear_vlasov_ampere accumulation #705 and can be replaced by it once CUDA version of the linear_vlasov_ampere accumulation #705 is merged.
• On a GPU, the same check runs directly (skipped here).

Model-specific changes:

None. Still missing for #689: the device Kronecker solve (feectools#96 / #703, so that applying the preconditioner makes no host copies), CG inner products returning host scalars, and an H100 run.

Documentation changes:

CUDA_STRATEGY.md: a short note in the feectools section on where the preconditioner data live and how they relate to feectools#96.

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

	before (devel preconditioner.py)	after
test_preconditioner_cupy.py (fake CuPy)	2 failed (AssertionError in is_circulant), 2 skipped (GPU)	2 passed, 2 skipped (GPU)
test_preconditioner_transpose.py	3 passed	3 passed
test_mass_matrices.py::test_mass_preconditioner, ::test_mass_preconditioner_array_weights_mpi, test_poisson.py::test_poisson_1d, test_curl_curl.py::test_convergence_1d, test_gyrokinetic_poisson.py::test_poisson_M1perp_1d, test_saddlepoint_massmatrices.py::test_saddlepointsolver_uzawa_small	not rerun	all passed
total of the run above		80 passed, 2 skipped (24.5 min)

Not run locally: test_mass_preconditioner_polar (more than 30 min per case here, with the venv's uncompiled feectools kernels) and the MPI runs. CI covers them.

🤖 Generated with Claude Code