Skip to content

Write down the correction-solve contract, and keep the sweeper's side of it - #786

Merged
pancetta merged 2 commits into
masterfrom
delta-solve-contract
Oct 8, 2026
Merged

pancetta merged 2 commits into
masterfrom
delta-solve-contract

Conversation

@pancetta

@pancetta pancetta commented Oct 8, 2026 •

Copy link
Copy Markdown
Member

The delta-form sweepers exist so that the implicit solve -- typically the expensive, often external part -- can work in, and return, a much lower precision while pySDC still converges to full precision. This PR writes down what the sweeper and the solver guarantee each other for that to hold, and makes the sweeper keep its side. It is the second step in moving the reduced-precision study (#784) into the library, after #785.

The contract

It is documented on Problem, replacing the note on the two optional delta-form methods, since a problem author looks there.

Routes, both existing:

  • the stock solve_system with linear_implicit=True, for linear or affine operators;
  • solve_system_delta(rhs, factor, base, f_base, t) for nonlinear ones. This is the contract for an external nonlinear solver.

The sweeper guarantees the solver:

  • a correction equation, with a right-hand side that shrinks with the iteration and a zero initial guess;
  • on the linear route, a right-hand side of unit size (new);
  • base and f_base at the level's precision, read only, and factor a plain scalar;
  • everything outside the solve kept at the level's precision, so each sweep is a step of iterative refinement around the solve.

The solver guarantees the sweeper:

  • the correction, in any precision (new: the sweeper widens it);
  • a relative accuracy η below what the method needs;
  • relative tolerances only, and none below what its own precision can deliver;
  • no boundary data or affine terms re-imposed on the correction.

The η values measured so far are given: SDC about 1e-4, an MLSDC fine level about 1e-6, an MLSDC coarse level about 1e-2.

Code

DeltaFormMixin._solve_correction's linear route goes through a new _solve_scaled, which:

  • divides the right-hand side by its maximum and multiplies the result back. That is exact for a linear solve, and keeps a reduced-precision format clear of underflow: the correction it hands over shrinks below float16's smallest subnormal. It also turns an absolute solver tolerance into a relative one;
  • removes the affine part after scaling. The stock solve adds it back unscaled, and the two only cancel at the same scale. The existing test_linear_implicit_handles_affine_operator caught the first version, which got this backwards;
  • widens a result that comes back in another precision before scaling it, so the multiplication never happens at the solver's precision. This applies only to datatypes whose init names a dtype; PETSc and FEniCS types keep today's copy;
  • returns a zero correction for a zero right-hand side without calling the solver.

At full precision the results are unchanged up to round-off: the scaling costs two roundings per solve. The node-parallel sweeper reuses _solve_correction, so it gets the same behaviour.

Decisions taken

  • No scale argument on solve_system_delta. On the nonlinear route the scaled unknown enters f, so the sweeper can't scale for the solver. Instead, the contract says the solver scales internally if its format needs it, so existing implementations are unaffected.
  • The required η is documented, not passed to the solver. The values are measured on the heat equation only; passing them as parameters can wait until they are known across problems and methods.
  • solve_jacobian is not a third route, for now.

Tests

  • test_half_precision_solver_reaches_full_precision: a heat problem whose solver rounds the right-hand side to float16, solves in single precision and returns float16. In delta form it reaches the double-precision floor; stock SDC with it stalls above 1e-6. With the scaling removed the test fails, at a floor of 1.2e-5.
  • test_half_precision_result_is_widened: the stored corrections are float64.
  • On a GPU, test_half_precision_solve_reaches_double_precision_in_delta_form (from Let the finite-difference problems solve at a precision of their own, and by FFT #785, CPU-only until now, where fp16 is emulated) also runs under the cupy marker. There the solve is cuFFT's genuine half-precision arithmetic, and delta-form SDC must still reach the double-precision floor while stock SDC with the same solve stalls. It passes on 4×T4 through etc/modal_gpu_tests.py. By hand at 256²: 4.9e-14 against 6.2e-14 in double, and the stock control stalls at 1.7e-1.

Local: black and ruff over pySDC; 1543 base tests across sweepers, transfers, problems, tutorials, test_core, test_reduced_precision and test_problem; and the delta-form MPI tests (test_delta_form_MPI, test_BaseTransferDelta_MPI) at every rank count they declare.

The docstrings of delta_form, delta_form_MPI and BaseTransferDelta now write the collocation residual as r, keeping ε for a format's precision, as the DeltaSDC README does.

🤖 Generated with Claude Code

… of it

The delta form exists so that the implicit solve -- typically the expensive,
often external part -- can work in, and return, a much lower precision while
the iteration still converges to full precision. What the sweeper and the
solver guarantee each other for that to hold is now documented on Problem,
replacing the note on the two optional methods.

The sweeper's side, on the linear route: it hands the solve its right-hand
side at unit size and multiplies the result back, which is exact and keeps any
format clear of underflow -- the correction it hands over shrinks with the
iteration, below float16's smallest subnormal by 1e-10 -- and it widens a
result that comes back in any precision before scaling it, so scaling never
happens at the solver's precision. The affine part f(0, t) is removed after the
scaling: the stock solve adds it back unscaled, and the two only cancel at the
same scale, which an existing test caught. On the nonlinear route the scaled
unknown would enter f, so there the solver scales internally, as documented.

Tested with a solver that rounds to float16, solves in single precision and
returns float16: in delta form it reaches the double-precision floor, and
fails at 1.2e-5 with the scaling removed; stock SDC with it stalls. On a T4,
delta-form SDC with #785's genuine half-precision FFT solve reaches 4.9e-14
against 6.2e-14 in double.

The docstrings write the collocation residual as r, keeping epsilon for a
format's precision, as the DeltaSDC README now does.

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

github-actions Bot commented Oct 8, 2026

Copy link
Copy Markdown

Coverage report

Click to see where and how coverage changed

FileStatementsMissingCoverageCoverage
(new stmts)
Lines missing
  pySDC/core
  problem.py
  pySDC/implementations/sweeper_classes
  delta_form.py
  delta_form_MPI.py
Project Total  

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

…metic

test_half_precision_solve_reaches_double_precision_in_delta_form ran on the
CPU only, where a float16 FFT solve is an emulation around a single-precision
transform. It now also runs under the cupy marker, where the solve is cuFFT's
half-precision arithmetic -- the case the correction-solve contract is written
for: delta-form SDC must still reach the double-precision floor, and stock SDC
with the same solve must still stall. One parametrisation, the same assertions.

Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
@pancetta
pancetta merged commit 1d1f143 into master Oct 8, 2026
73 checks passed
@pancetta
pancetta deleted the delta-solve-contract branch October 8, 2026 22:05
pancetta added a commit that referenced this pull request Oct 8, 2026
Brings in #785 (solve_precision and the FFT solve on the finite-difference
problems) and #786 (the correction-solve contract), which the project moves
onto next.

Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
pancetta added a commit that referenced this pull request Oct 8, 2026
…lve contract

#785 gave the finite-difference problems a solve precision of their own and an
FFT solve, and #786 made the delta-form sweeper scale the linear route's
right-hand side for any solver. The project's own versions of both go:

- heat_delta keeps only what no library solver has, half precision on a sparse
  solver, which it emulates; every other precision is the library's genuine
  solve. Its CG and FFT routes and its `normalize` knob are gone.
- heat_paradiag goes entirely: ParaDiag runs on the stock heat problem with
  dtype='complex128' and the library's solve_precision. paradiag.py keeps the
  transform at reduced precision, and circulant_eigenvalues goes with it.
- run_gpu.py uses heatNd_unforced and heatNd_forced directly.
- The project's FFT tests, which the library's test_FD_solve_precision now
  covers, and the GPU fingerprint of the FFT solve go.

The heat ladder's "unnormalised fp16 solve" control no longer fails, because
the sweeper now scales for that solver -- the contract working. It is removed
from the table, and lives on in the library's test with the scaling switched
off. The heat figure's control becomes stock SDC with the same half-precision
solve (it stalls at 0.25); Allen-Cahn keeps its unscaled solver, since on the
nonlinear route the solver scales itself. The README says so.

Every iteration count in the heat table is unchanged; the fp32 rows' floors move
within round-off (7.2e-14 to 8.8e-14), their solves now genuine rather than
emulated. ParaDiag's table is unchanged.

Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment

Labels

None yet

Projects

None yet

Development

Successfully merging this pull request may close these issues.

1 participant