[pancetta]

Part of moving the reduced-precision study (#784) into pySDC proper rather than keeping it in project-level problem twins. The goal it serves: showing that an expensive implicit solver may run, and return its result, at significantly reduced precision while pySDC still converges to full precision. This PR gives the library a solver that genuinely does so, independently of how the state is stored.

What changes

GenericNDimFinDiff, and through it heatNd_unforced, heatNd_forced and advectionNd, gains two options:

• solve_precision: the precision the implicit solve runs at, independently of the state's dtype. None, the default, runs the existing solve unchanged. When set, the solve genuinely works in that precision: 'direct', 'GMRES' and 'CG' on a copy of the operator held at it, and the result is widened back to dtype.
• solver_type='FFT', for periodic grids. The operator is circulant, so the FFT diagonalises it and a solve is two transforms and a division. That is exact, works for non-symmetric operators such as advection, and is the cheapest solve on a GPU.

Half precision is only reachable through 'FFT', since SciPy and CuPy sparse matrices stop at single precision:

• On a GPU the transforms are cuFFT's complex32 ones (power-of-two grids), with every value stored and every product rounded in half precision.
• On a CPU there is no half-precision FFT, so float16 is emulated: values are rounded through it around a single-precision transform. The docstring says so; it is optimistic.
• Range: the right-hand side is scaled to a maximum of one and the transforms are unitary, which keeps everything inside float16's range. An unscaled 1024² transform reaches 1e6; float16 stops at 65504.

Complex states solve at complex precision. That is ParaDiag's case: complex128 state, complex64 solve, complex factor.

Two fixes along the way:

• NaN from reduced-precision Krylov solves: such a solve is not asked for less than 100·ε of its own precision. Asked for more, fp32 CG on a T4 divided by a dot product that had rounded to zero and returned NaN.
• advectionNd now passes dtype and useGPU through. It never did, so the GPU project's statement that every FD problem became GPU-capable did not hold for advection.

Tests

New pySDC/tests/test_problems/test_FD_solve_precision.py:

• Accuracy at every precision, for every solver type: heat (symmetric) and advection (non-symmetric), in 1D and 2D, against a direct double-precision solve. Each check has a lower bound as well as an upper one: a solve that silently stayed in double would pass the upper bound, so only the lower bound shows the reduced precision is real.
• The complex ParaDiag case, and the two misconfigurations: float16 with a sparse solver, and 'FFT' on a non-periodic grid.
• What the option is for: delta-form SDC (delta_implicit, linear_implicit=True) with a half-precision FFT solve reaches the double-precision residual floor. Stock generic_implicit with the same solve stalls above 1e-6. That second case is the control.
• GPU (cupy): fp32 CG and genuine fp16 FFT on the device.

Local: black and ruff over pySDC, and 1594 base tests across problems, sweepers, transfers, tutorials and the problem helper. On 4×T4 through etc/modal_gpu_tests.py: the 17 cupy tests of test_problems pass.

🤖 Generated with Claude Code

[github-actions]

Coverage report

Click to see where and how coverage changed

File	Statements	Missing	Coverage	Coverage (new stmts)	Lines missing
pySDC/implementations/problem_classes
AdvectionEquation_ND_FD.py
HeatEquation_ND_FD.py
generic_ND_FD.py					456, 499
Project Total

This report was generated by python-coverage-comment-action