ForceFreeStates - Gaussian reduction trigger misses near-parallel solutions

[d-burg]

Summary

In the standard (forward) Euler–Lagrange integration of ForceFreeStates (FFS), the Gaussian reduction fires when the spread of the solution columns' growth since the previous reduction exceeds ucrit. That trigger cannot see columns that grow together and turn nearly parallel. By the time it fires, the reduction loses digits separating them, and the free-boundary energies carry an error that depends on where the ODE steps happened to land.

At the default and deck values of ucrit (1e3 on the Solovev decks, 1e4 elsewhere), that error is up to ~1e-5 relative, with no trend in the Euler–Lagrange tolerance. At ucrit = 1e5 it reaches 20 %.

Evidence

Case: examples/Solovev_ideal_example, n = 1.

• Reference: et[1] = 0.6773399618, from Vern9 at reltol 1e-12, abstol 1e-14.
• Agreement: every configuration below that converges agrees with the reference to ≤ 3e-10.
• Code: develop.

Relative error of et[1]:

eulerlagrange_tolerance	ucrit = 1e2	ucrit = 1e3 (deck)
1e-7	3.3e-8	−2.6e-7
1e-8	−5.2e-9	4.6e-10
1e-9	3.5e-10	−1.1e-9
5e-11	−1.2e-10	−2.0e-5
8e-11	−1.5e-10	9.8e-7
1e-10	−8.9e-11	5.3e-10
1.2e-10	−1.2e-10	−5.4e-6
1.4e-10	−2.6e-10	−8.0e-6
1e-11	−1.4e-11	−1.9e-6

At ucrit = 1e2 the error falls steadily as the tolerance tightens. At 1e3, tightening the tolerance does not reliably help. ucrit = 1e4 behaves like 1e3, and 1e5 gives a 20 % error.

Why the trigger misses it. Instrumenting the reduction shows the mechanism. At the two rational-surface crossings the growth spread was only 35–250, while every column had grown by 1e8–1e31 since the previous reduction. The columns were dominated by the same fastest-growing solution, and the norm-spread trigger cannot detect that.

The same holds on #480. #480 makes ucrit = 1e3 the default everywhere and refreshes a per-column abstol every step. On its head (3d6acd813), same deck and reference:

eulerlagrange_tolerance	ucrit = 1e2	ucrit = 1e3 (#480 default)
1e-9	3.1e-8	−1.8e-6
5e-11	7.5e-9	9.5e-6
8e-11	1.4e-10	−5.8e-6
1e-10	−1.4e-10	−8.3e-10
1.4e-10	4.1e-9	−7.6e-6
1e-11	−6.6e-10	2.4e-6

So the per-column tolerance does not remove this error. It comes from the reduction trigger, not from the ODE tolerance.

To reproduce, give a harness case an [overrides] table with "ForceFreeStates.ucrit" and "ForceFreeStates.eulerlagrange_tolerance" on the Solovev deck, and compare et[1] across tolerances.

Temporary fix

1. Merge the crossing fix, ForceFreeStates - BUGFIX - Remove the resonant solution at ideal crossings regardless of reduction timing #489, first. Lowering ucrit without it is unsafe. Develop at ucrit = 2–5 removes the wrong solution at a crossing whenever a reduction lands a step or two before it, and gives et[1] between −30.8 and +3.4 on this deck.
2. Then lower ucrit to 1e2.
   • On the Solovev decks this is the table's left column: clean convergence at no measurable runtime cost, since the step counts are the same.
   • The DIII-D decks and the default need the same tolerance scan before they change. They are 1e4 on develop and 1e3 under ForceFreeStates - BUGFIX! - Scale the ODE absolute tolerance to each solution column #480. That scan has not been run yet.

Longer term

• A better trigger: make the reduction track loss of linear independence rather than norm spread. For example, trigger on the growth since the last reduction (maximum(unorm)), or on a condition estimate of the first-component columns. Periodic orthonormalization is another option.
• Docstring: on develop the ucrit docstring says lower values give "excess renorms without accuracy gain". That does not hold for the standard path.
• Not measured: the Riccati path uses ucrit through renormalize_riccati_inplace!, a different mechanism, and was not measured here.