Skip to content

Add MIN-SR-S and MIN-SR-FLEX diagonal sweepers to SDC - #4347

Draft
singhharsh1708 wants to merge 3 commits into
SciML:masterfrom
singhharsh1708:sdc-min-sr
Draft

Add MIN-SR-S and MIN-SR-FLEX diagonal sweepers to SDC#4347
singhharsh1708 wants to merge 3 commits into
SciML:masterfrom
singhharsh1708:sdc-min-sr

Conversation

@singhharsh1708

Copy link
Copy Markdown
Contributor

Stacked on #4208 and the adaptivity PR. Until those merge the diff here also
shows their commits; the change in this one is the single commit "Add MIN-SR-S
and MIN-SR-FLEX diagonal sweepers".

Adds the two remaining diagonal preconditioners from Čaklović, Lunet, Götschel
and Ruprecht, SIAM Journal on Scientific Computing 47 (2025) A430-A453. A
diagonal removes the only node-to-node coupling inside a sweep, so these are
what make parallel across the nodes worth doing; they run serially here.

SDCSweeper.MIN_SR_S makes the stiff-limit iteration matrix I - QΔ⁻¹Q
nilpotent. SDCSweeper.MIN_SR_FLEX uses diag(τ)/k on sweep k and falls back
to MIN_SR_S past sweep M.

Implementation notes

MIN_SR_S has no closed form. Its entries solve a nilpotency condition
numerically, seeded incrementally in M, so rather than port that solve into a
tableau constructor they are tabulated for 2 to 10 nodes across both node
distributions and all four quadrature types, which is how the library ships
tableau data. Asking for more nodes raises an error naming the limit. The
residual spectral radius grows with M, from about 1e-8 at M = 2 to 1e-1 at
M = 10, because the paper offers no guarantee that an exact solution exists.

MIN_SR_FLEX needed the preconditioner to vary per sweep, so the tableau now
holds Vector{Matrix} accessed through sdc_qdelta_for(tab, k). Every other
sweeper stores a single entry, so nothing else pays for it.

Testing

Coefficients match qmat to 1e-13, including the fallback past sweep M.

The tests assert the defining properties rather than only the numbers:
(Q - QΔ)^M vanishes for MIN_SR_NS (Theorem 2.8), the stiff-limit spectral
radius is near zero for MIN_SR_S (Definition 2.11), and the product of the
stiff-limit iteration matrices over sweeps 1 to M vanishes for MIN_SR_FLEX
(Theorem 2.13).

Behaviour on Prothero-Robinson at lambda = -1e4, M = 3, dt = 0.1:

sweeper K=3 K=4 K=6 K=8
LU (serial) 3.6e-3 8.2e-4 2.9e-7 1.2e-9
MIN_SR_S 2.3e+0 3.1e-3 2.0e-6 1.3e-9
MIN_SR_FLEX 7.9e-2 5.7e-2 1.8e-4 2.6e-7
BEpar 1.4e+1 9.9e+0 1.5e+0 2.7e-1

MIN_SR_S reaches serial LU accuracy once the iteration has converged while
staying decoupled across the nodes, which is the claim of the paper. BEpar, the
naive diagonal choice, never gets there. MIN_SR_NS diverges on this problem, as
expected for a preconditioner optimised for the non-stiff limit, and there is a
test asserting that so the distinction does not later look like a bug.

AI Disclosure

Claude assisted with this change.

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