ForceFreeStates - BUGFIX - Decouple the ideal-crossing jump width from singfac_min

Aperta
  • Qualcosa non va?
#481 3 commenti 0 reazioni 1 assegnatario Vedi su GitHub

I maintainer di solito rispondono entro 3 giorni

@jhalpern30 ci sta già lavorando.

Dal 6/10/2026.

Valutazione

Difficoltà
5/5
Tempo stimato
Più di una settimana
Idoneità per principianti
30/100
Tipo di issue
Bug
Chiarezza
Abbastanza chiara
Stato di attività
Attiva
Stack tecnologico
julia
Ambito
backend

Direzione di ricerca

Start by reading cross_ideal_singular_surf! and the sing_get_ca / sing_get_ua routines mentioned in the issue, along with the discussion on #93. Compare the candidate approaches for preserving U₁†U₂ − U₂†U₁ = 0 across each crossing. Done means the crossing no longer introduces a persistent non-Hermiticity defect; the issue does not specify a chosen approach or a test file.

Scritto dal modello di indicizzazione a partire dal testo della issue.

Descrizione

Rewritten 2026-10-13. The original body blamed the trapezoid step and proposed re-projection, ODE sub-steps or extended precision; none of those addresses the actual source. The original text is in the edit history.

Problem

At each ideal singular surface, both integrators (forward and Riccati) stop at ψ_s − d and jump to ψ_s + d, where d = singfac_min/(n·|q′|). The jump zeroes the resonant component of the regular solutions, but that component only vanishes like x at the surface. Zeroing it at x = d instead of 0 leaves an error in W, δW and the energies that is first order in singfac_min:

Case singfac_min 1e-4 (default) 1e-6 1e-8
DIIID-like n=1 forward et₁ 0.801550 0.798864 0.798834
  • 100× smaller jump → 100× smaller change in et₁ (DIIID 2.69e-3 → 3.0e-5; Solovev 3.40e-3 → 3.4e-5), so first order, as expected.
  • Riccati et₁ drifts by 0.7% (n=1) to 10% (n=6) on DIIID-like at the default.
  • The non-Hermiticity of W after each crossing (the original symptom here) is the same truncation error.

Why lowering singfac_min isn't the fix

singfac_min sets four things at once:

  1. the jump half-width (wants to be small);
  2. the Riccati Δ′ matching distance (wants to be large: precision degrades as tol·x^(−2α); DIIID Riccati breaks at 1e-7);
  3. the chunk ends and kinetic chunking;
  4. the knot-free equilibrium grid bracket around each rational (GridRefinement).

After #520 the forward path stays accurate at small singfac_min (Solovev to 1e-8), so forward-only runs can lower it as a workaround, at the cost of moving the grid, the per-surface ca_left/ca_right stubs and kinetic chunking. Riccati cannot.

Proposed fix

Give the jump its own distance and leave singfac_min for everything else. Two candidate designs; prototype both and pick one:

A. ODE legs (prototyped). Add δ_jump (default 1e-6/(min(n)·|q′|)). Inside both crossings (cross_ideal_singular_surf!, riccati_cross_ideal_singular_surf!), integrate the same ODE from ψ_s − d to ψ_s − δ_jump, jump across ±δ_jump, and integrate back out to ψ_s + d. Every handoff point stays at d.

B. Crossing in asymptotic-coefficient space (not prototyped). Write u = V c_A (Glasser 2016 Eq. 50), set the large resonant coefficient to exactly zero, carry the non-resonant coefficients across (Glasser–Kolemen 2018 Eq. 36) and rebuild with ua on the right. This removes the O(d) error at today's 1e-4 with no legs and no new key, and also the d^(2α) residual that resonant-row elimination leaves for α < ½. Open risk: extracting coefficients through an ill-conditioned ua.

Prototype evidence (design A, δ_jump = 1e-6)

Case Today (et₁ / flagged steps) With fix
Forward DIIID n=1 0.801550 / 14 0.798864 / 0
Forward multi-n DIIID −0.575053 / 64 −0.575094 / 0
Riccati DIIID n=1 0.804133 0.798892
  • DIIID-like n=1–6: ≤ 1.1e-3 error at δ_jump = 1e-6; forward and Riccati extrapolations in δ_jump agree within 2e-5.
  • Riccati Δ′ matrix is bit-identical in every case (the STRIDE BVP reads only ua at the matching distance and the chunk propagators). Kinetic is identical.
  • Cost: +5–10% runtime at 1e-6.
  • Developer page with figures and the draft diff: https://claude.ai/artifact/MaMkwYpbRxALD99cT6GAWB

Caveat: the prototype predates #520's final trigger (cond(U₁) on the ucrit key, U₂ / lu(U₁)); re-run the table on #520's head. The earlier "DIIID n=3 forward 2.7× wrong" was mostly forward basis collapse, which #520 fixes, not this jump error (Riccati at the default is within ~1%).

Required corrections before a PR (design A)

  • Capture ca_left/ca_right and the Δca stub at d, outside the legs. They go to HDF5 and SLAYER's fallback, and ua is far worse conditioned at δ_jump.
  • Riccati: the column zeroing assumes (S, I) form, and the left leg ends without a renormalization, so the large solution is effectively removed at d. Renormalize at ψ_s − δ_jump before zeroing, or reuse #489's resonant-row elimination. This is the likely cause of the remaining 7e-5 forward/Riccati gap at DIIID n=3.
  • Scale δ_jump with min(n) in multi-n.
  • Attribute non-Hermiticity flags with a fixed axis start (frobenius_psi_max = 0); the default Frobenius start adds its own m=0 floor (#464).
  • The kinetic crossing (cross_kinetic_singular_surf!) has the same structure at singfac_min; its convergence is untested. Leave it unchanged and say so, or open a follow-up issue.
  • Check reverse-shear double rationals (#492) where two surfaces sit within d.

Impact

  • Every free-boundary energy moves toward convergence by roughly singfac_min relative (DIIID et₁ −0.34%), so the PR takes !. Shipped Fortran inputs use 1e-4, so results move away from Fortran benchmarks by about that amount.
  • Re-pin the harness and add a small-jump Solovev case.

Dependencies

  • #489 (crossing pivots on the resonant row) → #520 (cond(U₁) reduction trigger, resizable reduction storage, stable W). Both open; XB stacks on #520. Without #520, a small jump collapses the forward basis after high-α surfaces.
  • Independent of the Δ′-input work (#506, #517, field-line and grid PRs), since Δ′ is unchanged by construction. Those PRs do move the equilibrium, so take the final before/after numbers on the base this lands on. If design B is chosen, land #517 (trace-aware exponents) first, since B relies on ua.
  • The grid-bracket PR must keep its bracket ≥ d so the matching points and legs stay in one cubic cell.

Related

  • #93: the original report, closed by #480.
  • #490: forward basis collapse undetected by ucrit; fixed by #520 (open).
  • #464: Frobenius axis-start floor; separate from this.
Lingua principale
Julia
Stelle
11
Fork
7
Merge medio
3g 14h
PR unite (30g)
69

Preparare l'ambiente

  • Nessun Dockerfile né file Docker Compose
  • Ha un modello di pull request
  • Nessuna guida per i contributori

Come iniziare

  1. Leggi tutta la issue e poi la guida ai contributi del progetto.
  2. Commenta sulla issue per dire che te ne occupi tu — evita che due persone facciano lo stesso lavoro.
  3. Fai un fork del repository e lavora su un branch.
  4. Apri una pull request che faccia riferimento al numero della issue.

Altre issue di OpenFUSIONToolkit/GPEC

Tutte le issue di OpenFUSIONToolkit/GPEC

Issue simili

Altre issue su Julia

Ricevi le nuove issue nella tua casella

Un breve riepilogo di issue GitHub adatte ai principianti.