ForceFreeStates - Repair or retire the Frobenius axis start

Open
  • Is something wrong?
#464 2 comments 0 reactions 0 assignees View on GitHub

Maintainers usually reply within 3 days

Nobody has claimed this yet.

Assessment

Difficulty
5/5
Estimated time
Over a week
Newbie friendliness
30/100
Issue type
Bug
Clarity
Mostly clear
Activity status
Active
Tech stack
julia
Domain
backend

Research direction

Start by reading compute_axis_init and the investigations in #440 and #460; the issue identifies the diagonal 2×2 truncation and regular-branch selection as suspects. Verify the start against its documented asymptotic limit on an analytic equilibrium, or make it opt-in and remove that promise. After #460 merges, carry forward the fixture-based tests and promote the @test_broken lines when they pass.

Written by the indexing model from the issue text.

Description

Summary

compute_axis_init implements the Frobenius axis start [Glasser 2016 Eq. 51], whose docstring promises a solution behaving as ~ψ_low^(|m|/2). It does not deliver that limit at any starting surface anyone actually uses, and the resulting initial condition is wrong by enough to move — and on one configuration, to destroy — the free-boundary energies.

Two independent investigations have now measured this from different directions. #460 bounds where the start is used; neither it nor the superseded #440 repairs the start itself. This issue tracks the repair.

Evidence

Away from the axis (#460, DIII-D ideal example, against Fortran DCON):

ψ_low innermost coupling row, Julia/Fortran norm lowest W_t, Julia vs Fortran
1e-4 1.000 1e-4
0.01 1.010 0.03 %
0.03 1.042 0.1 %
0.10 1.243 2.3 %
0.30 2.140 7 %

At ψ_low = 0.01 (#440, synthetic diverted DIII-D-like ramp-up, q0 = 3.41, n = 1, no wall, ldp grid mpsi 128, riccati). The returned U₁ diagonal is 0.11…0.23 for m = −14…−3 and +6.1/−9.8/−13.7 for m = 3/4/5 — nowhere near ψ_low^(|m|/2). Consequences on that equilibrium:

  • one eigenvalue of W_p passes through zero, and et[1] reports a spurious −1e5 "instability" where Fortran DCON is stable;
  • on a clean neighbouring slice the same stiff eigenvalue is +2.76e5 against DCON's +1.61e4, a factor 17.

The second point is the important one: the eigenvalue is wrong everywhere in that scan, and the apparent "pole" is merely where it crosses zero. A threshold on the starting ψ does not address this.

Measured on #460's head (e93df8f0) with #440's fixture, which starts at exactly ψ_N = 0.01 and so stays on the Frobenius path under the frobenius_psi_max = 0.01 default:

setting et[1] max eig W_p negative W_p eigs
default (frobenius_psi_max = 0.01) −1.29e5 4.70e3 1
frobenius_psi_max = 0 (fixed start) +1.504 1.6090e4 0

Fortran DCON on the same deck, same mlow/mhigh/psilim/dmlim: stiff W_p +1.6092e4, no negative eigenvalue.

Suspects

Both named in #440, neither investigated:

  • the diagonal 2×2 truncation of A₀ — compute_axis_init builds one 2×2 Frobenius block per mode from the diagonal entries of F̄⁻¹, K and G, discarding all mode coupling;
  • the regular-branch selection rule — the eigenvector of the 2×2 block is picked by dominant |U₁| component, with a fallback to U₁ = 1, U₂ = 0 when the second component underflows.

Resolution

Either the start delivers its documented asymptotic limit, verified against that limit on an analytic equilibrium, or it is demoted to an opt-in with the docstring's promise removed. The present situation — a start that is documented as asymptotically correct, used by default near the axis, and measurably not — is the part to end.

Once it is repaired, the @test_broken lines in the fixture-based test (to be carried forward from #440 after #460 merges) report an unexpected pass and should be promoted.

Caveats on the numbers above

  • #440's fixture uses grid_type = "ldp" with mpsi = 128, not the current auto default; the ψ_low = 0.01 numbers should be re-confirmed on auto before anyone leans on them.
  • #440's regression table dates from 2026-09-02 and is stale against current develop. The two-row table in this issue was measured 2026-09-17 on #460's head.

References

  • #460 — bounds the Frobenius start by starting ψ (frobenius_psi_max); the mitigation, not the repair.
  • #440 — superseded and closed; carries the full investigation, the DCON crit(ψ) trace comparison and the scan figures.
Dominant language
Julia
Stars
11
Forks
7
Avg merge
3d 22h
Merged PRs (30d)
70

Getting set up

  • No Dockerfile or Docker Compose file
  • Has a pull request template
  • No contributing guide

First steps

  1. Read the whole issue, then the project's contributing guide.
  2. Comment on the issue to say you are picking it up — it saves two people doing the same work.
  3. Fork the repository and make your change on a branch.
  4. Open a pull request that references the issue number.

More from OpenFUSIONToolkit/GPEC

All issues in OpenFUSIONToolkit/GPEC

Similar issues

More Julia issues

Get new issues in your inbox

A short digest of beginner-friendly GitHub issues.