ForceFreeStates - BUGFIX! - Trigger the Gaussian reduction on the conditioning of U₁ and factor U₁ in the plasma response matrix

[jhalpern30]

Closes #490 and enables #481. Stacked on #489 (base bugfix/ffs-crossing-resonant-elimination); retarget to develop once #489 merges.

Accurate ideal energies close to rational surfaces, and trustworthy perturbed equilibria at high n.

Regime	Before → after	Case
Smaller singfac_min	8500× → 1e-6	Solovev et₁ at singfac_min 1e-8. Before, anything below ~1e-5 was unusable on this equilibrium; now accurate down to 1e-8. Prerequisite for #481, whose crossing defect scales with singfac_min.
High-α surfaces	4.1 → 7e-8	Solovev q₀ ≥ 1.85, singfac_min 1e-6
High-n perturbed equilibrium	1e95 → 2e-4	DIIID n = 6 and 8 overflowed the old 100-slot reduction storage and silently corrupted ξ, island widths and coupling outputs (n = 10 overflows too). Eigenvalues were unaffected.
Tolerance convergence (#490)	1e-3 → ≤ 1e-6	Solovev tolerance scan; ucrit now insensitive from 10 to 1e8

Standard DIIID-like decks at default settings move ≤ 1e-8, and the Riccati path is untouched. Errors are relative to the same build reducing at every step (ucrit = 1).

Review package with the figures, the #490 rerun, the ucrit ladder and the decision table: https://claude.ai/artifact/JqVHgVrx3xHbZLRJmqNAPe

Release note

• Audience: users
• Numerical impact: forward-integrator free-boundary energies and perturbed-equilibrium outputs. Near rational surfaces (small singfac_min) and at high-α surfaces, energies that were up to 8500× wrong now agree to ≤ 2e-6; at high n, perturbed-equilibrium outputs no longer overflow the reduction storage; regular decks move ≤ 1e-5 relative. (harness @ 1dc6138, 4-case gist; full run at merge time, see below)
• Migration: none. ucrit keeps its default of 1e3, but in the forward path it now means the U₁ condition-number limit (Riccati keeps the entry-size meaning).

You can now move the ideal jump much closer to rational surfaces (singfac_min down to 1e-8) and get accurate energies, and high-n perturbed-equilibrium runs no longer silently corrupt their output. The forward solve now applies a Gaussian reduction when the U₁ solution columns lose independence, not when their norms spread apart; the plasma response matrix W = U₂U₁⁻¹ is also formed stably.

Regression report

Gist only: 4 cases, d50a5e5b6 (#489 head) vs 1dc6138aa, trimmed to the headline rows. The full harness is deferred to merge time and must be run on the final head before merging, because this PR waits on #489 and a full report taken now would go stale. The harness stamp will not match later commits until then.

All regular-deck moves are at the ODE-node level (et₁ ≤ 1e-5 relative). Solovev Im(et₁) is a 2e-5 numerical zero on an under-resolved deck (etol 1e-7). Riccati is untouched.

Notes for reviewers

• Trigger: compute_solution_norms! reduces when column_condition(U₁) (1-norm estimate of the column-normalized U₁, from QR via cond(UpperTriangular(R), 1)) exceeds ucrit. Singular surfaces still always reduce.
• Storage: the reduction arrays grow by doubling (grow_reduction_storage!). The old 100-slot cap overflowed silently at n ≥ 6 and corrupted PE output (island widths off by 1e73).
• W = U₂ / lu(U₁) in free_run, free_compute_total and build_result: U₂ / U₁ factors U₁ᵀ, so partial pivoting runs across U₁'s column norms, which span ~1e110 just after a reduction. When the last event before the edge is a reduction, W was garbage (DIIID n = 4, δm = 24: et₁ −0.026 instead of 0.157). Older than this PR; the new trigger makes it more likely.
• Decisions to weigh in on: default ucrit 1e3 rather than 1e6 (1e6 is equally accurate in the forward path; 1e3 keeps Riccati, which shares the key, on its tested value); redefining ucrit rather than adding a key; cond(U₁) rather than cond[U₁;U₂]; keeping elimination rather than an (S, I) or QR reset; no dedicated test for the lu fix (no small matrix reproduces it).

[jhalpern30]

cc-ing @logan-nc @matt-pharr since this would basically fundamentally change the decades-old approach of how fixups are done in ForceFreeStates. After a lot of investigation/tests, it seems that the growth oriented reduction does not properly capture when a fixup is needed. Instead, the condition number of the displacement columns is more accurate, producing good numerics. Specifically, this fixes the Solov'ev case, which has a higher alpha (Frobenius coefficient in the expansion) at the singular surfaces, and was giving garbage as singfac_min was decreased

And on a related note, I've verified and tested that the 3D code also has this same issue (but worse) and is fixed by this method

[logan-nc]

I'm happy to implement tested improvemlike this! this! The fixup is the magic sauce, so having our own recipe is a good result of the rewrite effort

[matt-pharr]

@jhalpern30 be sure to post some plots highlighting the results, a table is nice but a plot is a much better record!

[logan-nc]

@matt-pharr figures are in the artifact linked in in the description

[d-burg]

The mechanism is convincing to me: columns that grow alike and turn parallel are invisible to a spread test, and the new unit tests separate the two triggers nicely. The storage overflow and the U₂ / lu(U₁) fix are clear bug fixes. I'm less sure the evidence is there yet for replacing the reduction trigger itself.

My and claude's brains merged into a disconcerting synthetic hybrid network and yielded the following items:

Blocking

1. Migration. The title carries ! and Migration starts with "none", which is what fails the title check, and substantively a forward-path deck with a hand-set ucrit is reinterpreted without a warning. Suggested text: "Migration: forward-integrator decks that set ucrit by hand should return it to the default 1e3. The value is now compared with a condition estimate of the column-normalized U₁, not with the spread of column-norm growth; small values such as 10 now reduce very often, and values above about 1e10 lose digits. Riccati decks need no change."
2. Regression report. The code block is empty and the stamp (1dc6138aa) predates the last two commits. This changes when every forward run reduces and how W is formed for every integrator, so I'd like the full harness now, d50a5e5b6 against the current head as two git refs, especially the Riccati, perturbed-equilibrium, kinetic Solovev and multi-n cases. Rerunning at merge time is cheap; reviewing blind is not.

Evidence I'd want for "changes the decades-old approach"

• Every error in the table is measured against the same build with ucrit = 1, which shows self-consistency, not correctness. One independent reference per headline row would help: the Riccati integrator on the same Solovev decks at singfac_min 1e-8, or a tolerance-converged or Fortran value where available. The figures sit behind a link "shared on request"; the numbers the argument rests on should be in the PR.
• Per case: accepted steps, number of reductions, the largest estimate at the moment of reduction, wall time and peak memory against #489's head. Without reduction counts nobody can judge 1e3 against 1e6, and with unbounded storage a run that reduces every other step is a memory problem the old 100-slot cap used to hide.
• Cases not yet covered: an ideal-unstable equilibrium where det U₁ passes through zero inside the plasma (U₁ is singular there for a physical reason, so the trigger fires whatever its threshold), a kinetic run, and the largest multi-n deck.

On your five decisions

1. 1e3 vs 1e6: undecided until the reduction counts are shown; leaning 1e6 for the forward path. The stated reason for 1e3 is about Riccati sharing the key, not about the forward path.
2. One key, two meanings: I'd add a forward-path key. Decks are reinterpreted silently, switching integrator changes what the same number means, and the shared key is already forcing a default neither path would pick on its own. The ! in the title makes a new key cheap now.
3. cond(U₁) rather than cond([U₁; U₂]): agree; the elimination pivots on U₁ and W inverts U₁.
4. Elimination rather than an (S, I) or QR reset: agree for this PR. Worth a follow-up issue: the elimination never rescales columns, which is why the norms span 1e110 in the first place.
5. No test for the lu fix: disagree; a 4×4 matrix reproduces it (inline).

Suggestion: the storage fix and the lu fix stand on their own and are easy to verify; the trigger change is a design decision that needs the evidence above. Splitting them would let the two bug fixes go in quickly on #489.

Small, outside the diff

• The OdeState field docstrings in EulerLagrange.jl still give the reduction arrays fixed shapes ((numpert_total, numunorms_init), "length numunorms_init"); they are initial capacities now. docs/src/stability.md still calls numunorms_init a "renorm checkpoint budget".
• The comment at Riccati/Driver.jl:77 still quotes wp = u[:,:,2] / u[:,:,1].
• The description's notes refer to column_condition, which the last two commits inlined.

[d-burg]

Two claude things about this trigger.

A 1-norm condition estimate of a column-normalised matrix has a floor that grows with the matrix size. Applying this elimination to random complex matrices leaves the estimate around 3e2 at N = 34 and above 1e3 from about N = 64, so for large numpert_total (multi-n, high n with a wide m range) the default 1e3 could fire on every second step. Real solution matrices are structured and may sit well below that, so this is a question: could you post odet.ifix and odet.step for each harness case, in particular the multi-n deck and DIIID n = 6/8/10? If ifix is close to step/2 anywhere, I'd compare against the value right after the last reduction, or scale the threshold with N, rather than raise the default.

It also does a QR plus about five N×N allocations on every accepted step, where the old test reused norms already computed. Standalone that is 0.05 ms at N = 34 and 7.7 ms at N = 400; probably small next to the right-hand side, but could you post benchmarks/benchmark_git_branches.jl numbers for the DIIID ideal and multi-n decks against #489's head? I'd also bring back the named helper the description mentions: it is the natural place for a reused work matrix and a direct unit test of the estimate.

[d-burg]

The claim checks out, and it is reproducible at 4×4: take L[i,j] = 1 on the diagonal and 0.8cis(i + 2j) below it, U1 = (L .* exp10.(range(110, 0; length=4))')[4:-1:1, :], A = [cis(0.3i*j)/(i+j) for i in 1:4, j in 1:4], W = A + A', and U2 = W * U1. Then U2 / U1 is off by 0.15 relative and U2 / lu(U1) by 2e-16 (other Hermitian W give 0.1 to 1). With random dense U₁ of the same column scales both forms are fine, which is probably why the search came up empty.

Since the expression now appears three times (here, line 279 and Result.jl:243), I'd put it in one small helper used by free_run, free_compute_total and build_result, and unit-test it with that matrix at rtol = 1e-12.

This path is shared with Riccati, so "the Riccati path is untouched" should be backed by the Riccati rows of the harness table. I expect rounding-level moves only.

[d-burg]

The docstring says the rounding error scales as eps·ucrit, but the check runs once per accepted step and is skipped on the first step after a reduction, so the conditioning when the elimination actually runs is at least ucrit and can be well above it where growth per step is large, which is the small-singfac_min regime this PR targets. Do you have the distribution of the estimate at the moment each reduction fires for the Solovev singfac_min = 1e-8 case? If the overshoot is routinely a few decades, I'd say "at least", and also that this is a 1-norm estimate whose healthy baseline grows with the number of modes.

[d-burg]

Doubling zero is still zero, so numunorms_init = 0 throws a BoundsError on the first reduction now that the docstring says the storage grows as needed.

Suggested change

	nfix = size(odet.fixfac, 3)
	nfix = max(size(odet.fixfac, 3), 1)

[d-burg]

Why were 2.0 and 5.0 dropped from this ladder? The rtol is unchanged, but removing the smallest values narrows what the test accepts, and 5.0 was the value that fails on develop for the #489 bug. The reference for your whole results table is the ucrit = 1 run, and the description claims insensitivity from 10 to 1e8. If 2.0 still passes at 1e-6, please keep it and add 1.0 and 1e8 so the test covers the claimed range; if it doesn't, we should understand that before trusting the reference.

[d-burg]

Nothing pins the results that motivate this PR. This testset uses the stock Solovev deck (singfac_min = 1e-4), where the old and new triggers agree. Since it already rewrites the deck per run, a second loop at singfac_min = 1e-8 checking et₁ across the same ucrit ladder is a few lines, and a harness case for one high-n perturbed-equilibrium quantity (an island width at n = 6) would pin the storage fix end to end.