Equilibrium - BUGFIX! - ⚠️ Recover the magnetic axis when psi is not smooth near it, without relaxing the convergence test by d-burg · Pull Request #444 · OpenFUSIONToolkit/GPEC

d-burg · 2026-09-05T18:13:03Z

Release note

Audience: users
Numerical impact: none (harness @ bf4e298)
Migration: smooth equilibria are unchanged. A file whose magnetic-axis search lands off the ψ grid or on an X-point (det ≤ 0) now raises an error instead of returning a wrong axis silently.

Two TJ-written 257×257 geqdsks (two points of an 81-point circular β scan) used to abort reading with Jacobian matrix is singular near (…) or Failed to find magnetic axis after 200 iterations. The cause is the files, not the spline: their ψ has a curvature spike in the grid cells at the axis (∂²ψ/∂R² ≈ 13–19 against ≈ 3–10 on the neighbouring cells, in the raw data), so Newton from the midplane guess meets a singular or indefinite Hessian. The axis search now recovers such files without relaxing what "found" means, and says loudly that the file looks wrong.

Regression report

Compared develop @ 6d655c74c against this branch @ bf4e298e1: diiid_n1 53 unchanged, diiid_n1_riccati 17 unchanged, efit_fixedbdy_separatrix 5 unchanged, diiid_slayer_n1 17 unchanged, and solovev_n1 22 unchanged. All numeric differences are zero; checksums are identical, including the SLAYER results in this run. Focused equilibrium, axis-Newton, and grid-refinement tests passed.

Run on ITER SDCC through Slurm, Julia 1.12.6, CI-pinned manifest, 2 Julia threads / 1 BLAS thread per case (16 threads for diiid_error_field).

Full regression report

================================================================
Case: diiid_n1 — DIII-D-like equilibrium, n=1, ideal + perturbed equilibrium
================================================================

Regression Report: diiid_n1
===========================================================================================================================
Ref 1: develop  @ 6d655c74c (2026-10-27)
       env: julia 1.12.6, x86_64-linux-gnu, manifest 4c9d54fa (pinned), 2 threads/1 BLAS, jit build
Ref 2: bugfix/equilibrium-axis-search-robustness  @ bf4e298e1 (2026-10-27)
       env: julia 1.12.6, x86_64-linux-gnu, manifest 4c9d54fa (pinned), 2 threads/1 BLAS, jit build
---------------------------------------------------------------------------------------------------------------------------
Quantity                                      develop          bugfix/equilibrium-axis-search-robustness  Diff       Status
---------------------------------------------------------------------------------------------------------------------------
total energy Re(et[1])                        8.006122e-01     8.006122e-01                               0.0e+00    OK    
total energy Im(et[1])                        7.924888e-05     7.924888e-05                               0.0e+00    OK    
plasma energy Re(ep[1])                       -1.372261e+00    -1.372261e+00                              0.0e+00    OK    
vacuum energy Re(ev[1])                       2.172874e+00     2.172874e+00                               0.0e+00    OK    
vacuum matrix min eigenvalue                  1.873975e-01     1.873975e-01                               0.0e+00    OK    
plasma energy (all)                           [35 elem]        [35 elem]                                  0.0e+00    OK    
vacuum energy (all)                           [35 elem]        [35 elem]                                  0.0e+00    OK    
total energy (all)                            [35 elem]        [35 elem]                                  0.0e+00    OK    
ODE steps (saved)                             792              792                                        0.0e+00    OK    
ODE steps (total)                             787              787                                        0.0e+00    OK    
q0                                            1.204212e+00     1.204212e+00                               0.0e+00    OK    
q95                                           4.781723e+00     4.781723e+00                               0.0e+00    OK    
beta_t                                        1.327024e-02     1.327024e-02                               0.0e+00    OK    
beta_n                                        1.372511e+00     1.372511e+00                               0.0e+00    OK    
internal inductance li1                       8.842268e-01     8.842268e-01                               0.0e+00    OK    
internal inductance li2                       7.080757e-01     7.080757e-01                               0.0e+00    OK    
internal inductance li3                       7.304340e-01     7.304340e-01                               0.0e+00    OK    
poloidal beta betap1                          6.680739e-01     6.680739e-01                               0.0e+00    OK    
poloidal beta betap2                          5.349836e-01     5.349836e-01                               0.0e+00    OK    
poloidal beta betap3                          5.518764e-01     5.518764e-01                               0.0e+00    OK    
# singular surfaces                           5                5                                          0.0e+00    OK    
singular psi locations                        [5 elem]         [5 elem]                                   0.0e+00    OK    
singular q values                             [5 elem]         [5 elem]                                   0.0e+00    OK    
current beta betaj                            4.236480e-01     4.236480e-01                               0.0e+00    OK    
plasma volume                                 1.829472e+01     1.829472e+01                               0.0e+00    OK    
plasma current                                1.152130e+00     1.152130e+00                               0.0e+00    OK    
mpert                                         35               35                                         0.0e+00    OK    
npert                                         1                1                                          0.0e+00    OK    
toroidal field bt0                            2.006573e+00     2.006573e+00                               0.0e+00    OK    
wall field bwall                              3.880145e-01     3.880145e-01                               0.0e+00    OK    
aspect ratio                                  2.845746e+00     2.845746e+00                               0.0e+00    OK    
elongation kappa                              1.708322e+00     1.708322e+00                               0.0e+00    OK    
q profile (checksum)                          6da947d0754e...  6da947d0754e...                            identical  OK    
pressure profile (checksum)                   c9c6ca4072ca...  c9c6ca4072ca...                            identical  OK    
Mercier D_I profile (checksum)                26babbe22505...  26babbe22505...                            identical  OK    
resistive interchange D_R profile (checksum)  c46830151a0b...  c46830151a0b...                            identical  OK    
ballooning Delta' profile (checksum)          c27c6779804a...  c27c6779804a...                            identical  OK    
island half-widths                            [5 elem]         [5 elem]                                   0.0e+00    OK    
Chirikov parameter                            [5 elem]         [5 elem]                                   0.0e+00    OK    
||resonant area-weighted field||              5.288461e-04     5.288461e-04                               0.0e+00    OK    
dominant-coupling singular values             [3 elem]         [3 elem]                                   0.0e+00    OK    
|forcing overlap with dominant mode|          1.411548e-04     1.411548e-04                               0.0e+00    OK    
|delta_nominal| of coil set 1                 7.034622e-05     7.034622e-05                               0.0e+00    OK    
||ddelta/d(shift)|| over coil sets            1.245446e-04     1.245446e-04                               0.0e+00    OK    
||ddelta/d(tilt)|| over coil sets             4.177337e-07     4.177337e-07                               0.0e+00    OK    
PE plasma energy                              3.320466e+00     3.320466e+00                               0.0e+00    OK    
PE vacuum energy                              3.174319e+00     3.174319e+00                               0.0e+00    OK    
PE surface energy                             5.508585e+00     5.508585e+00                               0.0e+00    OK    
PE toroidal torque                            1.119355e-04     1.119355e-04                               0.0e+00    OK    
NTV torque FGAR [N·m]                         5.480580e-01     5.480580e-01                               0.0e+00    OK    
NTV kinetic energy dW FGAR [J]                8.163220e-02     8.163220e-02                               0.0e+00    OK    
Runtime (s)                                   295.6s           295.5s                                                --    
||forcing b~|| (root-area-weighted)           4.683171e-04     4.683171e-04                               0.0e+00    OK    
resonant area-weighted field b^r              [5 elem]         [5 elem]                                   0.0e+00    OK    
===========================================================================================================================
Summary: 53 unchanged


================================================================
Case: diiid_n1_riccati — DIII-D-like equilibrium, n=1, Riccati integrator Δ' matrix
================================================================

Regression Report: diiid_n1_riccati
========================================================================================================
Ref 1: develop  @ 6d655c74c (2026-10-27)
       env: julia 1.12.6, x86_64-linux-gnu, manifest 4c9d54fa (pinned), 2 threads/1 BLAS, jit build
Ref 2: bugfix/equilibrium-axis-search-robustness  @ bf4e298e1 (2026-10-27)
       env: julia 1.12.6, x86_64-linux-gnu, manifest 4c9d54fa (pinned), 2 threads/1 BLAS, jit build
--------------------------------------------------------------------------------------------------------
Quantity                       develop        bugfix/equilibrium-axis-search-robustness  Diff     Status
--------------------------------------------------------------------------------------------------------
delta prime (BVP diagonal)     [5 elem]       [5 elem]                                   0.0e+00  OK    
delta prime (raw side-major)   [10 elem]      [10 elem]                                  0.0e+00  OK    
edge coil response delta_coil  [10 elem]      [10 elem]                                  0.0e+00  OK    
total energy Re(et[1])         8.031971e-01   8.031971e-01                               0.0e+00  OK    
plasma energy Re(ep[1])        -1.368534e+00  -1.368534e+00                              0.0e+00  OK    
vacuum energy Re(ev[1])        2.171732e+00   2.171732e+00                               0.0e+00  OK    
total energy (all)             [35 elem]      [35 elem]                                  0.0e+00  OK    
# singular surfaces            5              5                                          0.0e+00  OK    
singular psi locations         [5 elem]       [5 elem]                                   0.0e+00  OK    
singular q values              [5 elem]       [5 elem]                                   0.0e+00  OK    
ODE steps (saved)              95             95                                         0.0e+00  OK    
ODE steps (total)              2633           2633                                       0.0e+00  OK    
mpert                          35             35                                         0.0e+00  OK    
npert                          1              1                                          0.0e+00  OK    
q0                             1.204212e+00   1.204212e+00                               0.0e+00  OK    
q95                            4.781723e+00   4.781723e+00                               0.0e+00  OK    
beta_n                         1.372511e+00   1.372511e+00                               0.0e+00  OK    
Runtime (s)                    244.0s         250.2s                                              --    
========================================================================================================
Summary: 17 unchanged


================================================================
Case: efit_fixedbdy_separatrix — Fixed-boundary EFIT (box hugs LCFS), Newton separatrix-find regression
================================================================

Regression Report: efit_fixedbdy_separatrix
====================================================================================================
Ref 1: develop  @ 6d655c74c (2026-10-27)
       env: julia 1.12.6, x86_64-linux-gnu, manifest 4c9d54fa (pinned), 2 threads/1 BLAS, jit build
Ref 2: bugfix/equilibrium-axis-search-robustness  @ bf4e298e1 (2026-10-27)
       env: julia 1.12.6, x86_64-linux-gnu, manifest 4c9d54fa (pinned), 2 threads/1 BLAS, jit build
----------------------------------------------------------------------------------------------------
Quantity                    develop       bugfix/equilibrium-axis-search-robustness  Diff     Status
----------------------------------------------------------------------------------------------------
psio (total poloidal flux)  3.344390e-01  3.344390e-01                               0.0e+00  OK    
q0                          1.155202e+00  1.155202e+00                               0.0e+00  OK    
q95                         5.741022e+00  5.741022e+00                               0.0e+00  OK    
beta_t                      3.945172e-05  3.945172e-05                               0.0e+00  OK    
beta_n                      1.009505e-01  1.009505e-01                               0.0e+00  OK    
Runtime (s)                 62.3s         61.6s                                               --    
====================================================================================================
Summary: 5 unchanged


================================================================
Case: diiid_slayer_n1 — DIII-D-like H-mode equilibrium, n=1, SLAYER tearing-mode analysis (uncoupled per-surface, AMR, validity-gated)
================================================================

Regression Report: diiid_slayer_n1
========================================================================================================
Ref 1: develop  @ 6d655c74c (2026-10-27)
       env: julia 1.12.6, x86_64-linux-gnu, manifest 4c9d54fa (pinned), 2 threads/1 BLAS, jit build
Ref 2: bugfix/equilibrium-axis-search-robustness  @ bf4e298e1 (2026-10-27)
       env: julia 1.12.6, x86_64-linux-gnu, manifest 4c9d54fa (pinned), 2 threads/1 BLAS, jit build
--------------------------------------------------------------------------------------------------------
Quantity                            develop   bugfix/equilibrium-axis-search-robustness  Diff     Status
--------------------------------------------------------------------------------------------------------
SLAYER surface indices              [6 elem]  [6 elem]                                   0.0e+00  OK    
SLAYER poloidal m                   [6 elem]  [6 elem]                                   0.0e+00  OK    
SLAYER toroidal n                   [6 elem]  [6 elem]                                   0.0e+00  OK    
SLAYER minor radius rs              [6 elem]  [6 elem]                                   0.0e+00  OK    
SLAYER r-based shear                [6 elem]  [6 elem]                                   0.0e+00  OK    
SLAYER Lundquist S                  [6 elem]  [6 elem]                                   0.0e+00  OK    
SLAYER D_norm                       [6 elem]  [6 elem]                                   0.0e+00  OK    
SLAYER P_perp                       [6 elem]  [6 elem]                                   0.0e+00  OK    
SLAYER tauk                         [6 elem]  [6 elem]                                   0.0e+00  OK    
SLAYER iota_e                       [6 elem]  [6 elem]                                   0.0e+00  OK    
SLAYER Q_e                          [6 elem]  [6 elem]                                   0.0e+00  OK    
SLAYER Q_i                          [6 elem]  [6 elem]                                   0.0e+00  OK    
SLAYER Q_root [2/1,3/1,4/1]         [3 elem]  [3 elem]                                   0.0e+00  OK    
SLAYER ω_Hz [2/1,3/1,4/1]           [3 elem]  [3 elem]                                   0.0e+00  OK    
SLAYER γ_Hz [2/1,3/1,4/1]           [3 elem]  [3 elem]                                   0.0e+00  OK    
SLAYER no_root flags [2/1,3/1,4/1]  [3 elem]  [3 elem]                                   0.0e+00  OK    
SLAYER enabled flag                 1         1                                          0.0e+00  OK    
Runtime (s)                         277.0s    280.0s                                              --    
========================================================================================================
Summary: 17 unchanged


================================================================
Case: solovev_n1 — Solovev analytical equilibrium, n=1, ideal stability
================================================================

Regression Report: solovev_n1
===========================================================================================================
Ref 1: develop  @ 6d655c74c (2026-10-27)
       env: julia 1.12.6, x86_64-linux-gnu, manifest 4c9d54fa (pinned), 2 threads/1 BLAS, jit build
Ref 2: bugfix/equilibrium-axis-search-robustness  @ bf4e298e1 (2026-10-27)
       env: julia 1.12.6, x86_64-linux-gnu, manifest 4c9d54fa (pinned), 2 threads/1 BLAS, jit build
-----------------------------------------------------------------------------------------------------------
Quantity                      develop          bugfix/equilibrium-axis-search-robustness  Diff       Status
-----------------------------------------------------------------------------------------------------------
total energy Re(et[1])        1.090560e-02     1.090560e-02                               0.0e+00    OK    
total energy Im(et[1])        1.309258e-10     1.309258e-10                               0.0e+00    OK    
plasma energy Re(ep[1])       -1.041640e+01    -1.041640e+01                              0.0e+00    OK    
vacuum energy Re(ev[1])       1.042730e+01     1.042730e+01                               0.0e+00    OK    
vacuum matrix min eigenvalue  2.170703e+00     2.170703e+00                               0.0e+00    OK    
plasma energy (all)           [32 elem]        [32 elem]                                  0.0e+00    OK    
vacuum energy (all)           [32 elem]        [32 elem]                                  0.0e+00    OK    
total energy (all)            [32 elem]        [32 elem]                                  0.0e+00    OK    
ODE steps (saved)             287              287                                        0.0e+00    OK    
ODE steps (total)             285              285                                        0.0e+00    OK    
q0                            1.900006e+00     1.900006e+00                               0.0e+00    OK    
q95                           3.147422e+00     3.147422e+00                               0.0e+00    OK    
beta_t                        4.621027e-02     4.621027e-02                               0.0e+00    OK    
beta_n                        3.217775e+00     3.217775e+00                               0.0e+00    OK    
# singular surfaces           2                2                                          0.0e+00    OK    
singular psi locations        [2 elem]         [2 elem]                                   0.0e+00    OK    
singular q values             [2 elem]         [2 elem]                                   0.0e+00    OK    
mpert                         32               32                                         0.0e+00    OK    
npert                         1                1                                          0.0e+00    OK    
q profile (checksum)          2340b888d410...  2340b888d410...                            identical  OK    
pressure profile (checksum)   108b947087e6...  108b947087e6...                            identical  OK    
ca_left (checksum)            21c2d085d11b...  21c2d085d11b...                            identical  OK    
Runtime (s)                   242.8s           241.4s                                                --    
===========================================================================================================
Summary: 22 unchanged

What changed

Smooth files take exactly the old path. Plain, uncapped Newton from the midplane guess runs first, with the same step-size test (|δR|, |δZ| ≤ 1e-12·R). Nothing on that path changed, so healthy files are bit-identical by construction.
On failure, one fallback. Newton restarts, with each step capped at one march step, from a point located by alternating 1-D bisections on B_z (along R) and B_r (along Z). This uses first derivatives only. Each bisection checks that its bracket straddles a sign change, doubles the bracket only within the grid, and errors if it finds none.
Same acceptance test on both paths. The restart must pass the same 1e-12 step-size test. Otherwise Failed to find magnetic axis is a hard error again. No unconverged point is ever accepted. A successful fallback issues a @warn that ψ is likely not smooth near the axis.

Notes for reviewers

Premise corrected from the first version of this PR. That version described the files as "perfectly good" and the ψ-spline Hessian as unreliable. It also accepted the bisection point when Newton failed, with only an @info, which quietly changed the convergence criterion. An adversarial review found both problems, and both are fixed here. The spline is faithfully interpolating defective data, so the files themselves should be regenerated. That is for the TJ geqdsk writer to fix, not something this PR can do.
The axis the fallback finds sits ~3e-4 m from the header's (2.0, 0.0). It is at (2.0001, 2.4e-4) and (1.99974, 3.1e-4), with Z ≠ 0 despite up-down symmetry. It is a genuine converged critical point of the spline, and well inside one grid cell (3.4 mm). The offset is another symptom of the defect, which is why the test's 1e-3 bound comes with that reason.
Fixtures: 1.1 MB each, synthetic TJ output (not experimental data), provenance now recorded in test/test_data/README_test_data.md. They are the only files known to reproduce each failure mode.
Interaction with the local bugfix/build-reproducible-arithmetic branch: that branch edits the same Newton loop (an extra polish step after the stop test), so the two will conflict textually. The polish belongs inside _axis_newton so both attempts get it.

Audit round (2026-09-25)

Every returned axis is now checked, on both the plain and the fallback path. It must lie inside the ψ grid, and it must be an O-point (det ∂(B_R,B_Z)/∂(R,Z) > 0).
- The in-grid check catches a real failure. On a cusped synthetic grid, plain Newton converged silently to a false critical point 0.44 m outside the grid.
No midplane-march R-bracket test. An earlier head required the axis to lie inside the R-bracket of the midplane B_z sign change. The axis of an up-down asymmetric equilibrium need not share its R with that sign change, and the bracket rejected correct axes on both paths: a smooth Solovev ψ on a Z-shifted grid (plain Newton), and the same grid with a cusp at the axis (the restart found the axis to 0.05 cell). It is removed.
- The plain path accepts exactly what develop accepts, minus off-grid and X-point results. The restart is certified by the in-grid, O-point and residual checks.
Tests do not depend on which path Newton takes. Newton on the cusp fixture converges, cycles or leaves the grid depending on the machine's rounding: the same commit behaved differently on macOS arm64, the CI runners and an x86 cluster. The cusp test is now a sweep (3 grid placements × 9 cusp strengths) that requires a verified O-point at the cusp or a Failed to find magnetic axis error, with the fallback firing at least once. The off-grid error is tested on a smooth grid that excludes the axis.
The fallback has a residual bound: |B_p| ≤ 1e-10·B_scale, with B_scale = ψ_o/(R_o·ΔR). Measured |B_p|/B_scale on the fixtures is 4e-16 to 5e-14.
The fixtures are gone. The two 1.12 MB TJ geqdsks are replaced by an in-test Solovev equilibrium with a cusp at the axis, which reproduces the Newton cycling.
- The test now asserts physics: an O-point, |∇ψ| small against ψ_o/ΔR, Z = 0 to rounding on the symmetric grid, and the axis within a cell of the cusp.
Tidying: bisection uses Roots.find_zero(…, Bisection()), the numerical controls are named, the warnings state which condition stopped Newton, and the docstring's Fortran citation is gone.

NO MERGE WITHOUT HUMAN REVIEW. This PR is a draft and needs a named human reviewer and an assignee before it can be considered for merge. Non-negotiable.

🤖 Generated with Claude Code

github-actions · 2026-09-05T18:13:13Z

This pull request is missing an assignee and a reviewer.

If you are not ready to name them, mark this pull request as a draft.
docs/development/contributors.md suggests lead developers to ask.
Merging is not blocked here, but no pull request may be merged without human review.

…eliable psi-spline Hessian Two failure modes of `direct_position!` on otherwise healthy geqdsks, both caused by the 2-D psi spline's second derivatives near the axis (dB_z/dR wanders between ~5 and ~0 within a few cells, so the Hessian goes near-singular or indefinite): 1. "Jacobian matrix is singular" -- one bad iterate sent the undamped Newton step across the whole box, into a region where det J underflowed. 2. "Failed to find magnetic axis after 200 iterations" -- Newton cycled between two points either side of the axis on an indefinite Hessian. Fix, in two stages that leave the converging path untouched: * Cap each Newton step at one midplane-march step. Inactive for well-behaved iterations (their steps are < dr/2); it only keeps a bad iterate inside the axis neighbourhood until the Hessian recovers. A singular Hessian now returns "not converged" instead of erroring. * If Newton does not converge, locate the axis with first derivatives only: B_z changes sign across the axis along the midplane and B_r along a column, so four alternating 1-D bisections (80 iterations each, widening bracket in Z) converge to grad(psi) = 0 without touching the spline's second derivatives. Newton is retried from that point; if it still fails, the bisection point is used as the axis (it is a converged zero of the first derivatives) with an @info message instead of an error. Healthy equilibria never enter the fallback, so their results are bit-identical (regression harness: diiid_n1 47/47, diiid_slayer_n1 15/15, efit_fixedbdy_separatrix 5/5 unchanged vs 349a0c2). Two 257x257 TJ circular geqdsks reproducing each failure mode are added as fixtures with a test that both now locate the axis where the midplane B_z crossing says it is, and that a healthy geqdsk takes the plain Newton path. Co-Authored-By: Claude Fable 5.1 <noreply@anthropic.com>

…ic-axis search falls back The TJ fixtures fail Newton because their ψ has a curvature spike at the axis, not because the spline's Hessian is unreliable. So the search no longer accepts an unconverged point: - Plain, uncapped Newton runs first, exactly as before, so smooth files are bit-identical by construction rather than because a step cap happened to stay inactive. - On failure, capped Newton restarts from a point found by alternating 1-D bisections on B_z and B_r. The bisection now checks that each bracket straddles a sign change, stays inside the grid, and errors when it cannot find one. - The restart must pass the same step-size test; otherwise "Failed to find magnetic axis" is an error again. A successful fallback warns that ψ is likely not smooth near the axis. - The test asserts the warned fallback on both fixtures and no warning on a healthy file; the fixtures' provenance and defect are recorded in the test-data README. Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>

…t needs no Logging dependency Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>

…ck and report why Newton stopped - After the bisection-seeded Newton restart, evaluate det ∂(B_R,B_Z)/∂(R,Z) at the returned point and error unless it is positive (O-point; an X-point has det < 0). - The Newton helper returns :converged, :singular_jacobian or :iteration_cap, and the warning and error messages state which condition stopped it, plus |B_p| and det at the result, instead of asserting an indefinite Hessian that is never tested. - Replace the hand-written bisection with Roots.find_zero(..., Roots.Bisection()), whose float-bisection stop also terminates for brackets around Z = 0. - Name the singular-Jacobian threshold and the bisection pass count (values unchanged). - Drop the Fortran citation from the docstring and state the returned axis's contract. The plain Newton path is arithmetically unchanged. Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>

… their closure find_separatrix_crossing assigned r and dr, which rebound (and boxed) the enclosing function's variables that the axis search also uses. Declaring them local keeps the closure self-contained; the arithmetic is unchanged. Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>

…est cusped Solovev equilibrium Build the fallback fixture in the test: an up-down symmetric Solovev psi with a conical cusp at the axis, so the axis cells carry a curvature spike and Newton from the midplane march point does not converge. The tests now assert that the returned axis is an O-point of the psi spline, that |grad psi| is below the Newton step tolerance times psio/(rmax-rmin), that Z vanishes to rounding on the symmetric fixture, and that it lies within one grid cell of the cusp apex. A smooth variant and the COCOS_02 geqdsk must take the plain path with no warning. Drops the two 1.1 MB TJ geqdsks. Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>

…e the march bracket, or not O-points On both the plain and the fallback path, the magnetic axis must now lie inside the psi grid and, when the midplane march took a step, inside its R-bracket of the B_z sign change; it must also be an O-point (det d(B_R,B_Z)/d(R,Z) > 0). After the fallback, |B_p| must be at most 1e-10 B_scale with B_scale = |psio| / (R_axis (R_max - R_min)). Each violation is an error that reports where Newton landed and the bounds it broke. Plain Newton on a cusped psi could previously converge silently to a false O-point off the grid, found by spline extrapolation. Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>

…eted axis in the axis-search test The cusped-Solovev fixture takes explicit inboard and outboard grid extents. An inboard-heavy grid makes the midplane march bracket the axis, and cusp strengths 2.25 and 2.5 cells, where plain Newton converges off the grid, must now raise the off-grid error. Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>

Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>

logan-nc · 2026-10-03T22:47:45Z

@d-burg should this be assigned a reviewer?

…tarted axis search The axis of an up-down asymmetric equilibrium need not share its R with the midplane B_z sign change, so the bracket rejected a correct axis found by the direct Newton iteration. It now applies only to the restart. Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>

Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>

…arted axis search On an up-down shifted grid with a cusp at the axis the restart found the correct axis and the bracket rejected it. The restart is still certified by the in-grid, O-point and residual checks. Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>

…path a chaotic Newton takes Newton on the cusp fixture converges, cycles or leaves the grid depending on the machine's rounding, so tests that named the path failed on some hosts. The cusp sweep now asserts a verified axis or a loud failure, and the off-grid error uses a smooth grid that excludes the axis. Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>

matt-pharr · 2026-10-27T12:31:29Z

Merged develop and retained the axis-search changes in Axisymmetric/Direct.jl, leaving ψ-grid construction unchanged. Kept both test-list additions and updated the axis tests to use the renamed SolovevConfig.elongation field.

github-actions Bot added the bugfix Something was wrong and now is not label Sep 5, 2026

d-burg and others added 2 commits September 23, 2026 16:33

d-burg force-pushed the bugfix/equilibrium-axis-search-robustness branch from 7513e20 to a0b33f6 Compare September 23, 2026 22:11

d-burg changed the title ~~Equilibrium - BUGFIX - Make the magnetic-axis search robust to an unreliable psi-spline Hessian~~ Equilibrium - BUGFIX - Recover the magnetic axis when psi is not smooth near it, without relaxing the convergence test Sep 23, 2026

d-burg and others added 7 commits September 23, 2026 18:30

Equilibrium - TEST - Take the warning level from Base so the axis tes…

604fc10

…t needs no Logging dependency Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>

Repo - MINOR - Apply the formatter to the audit fixes

b32ffcb

Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>

d-burg changed the title ~~Equilibrium - BUGFIX - Recover the magnetic axis when psi is not smooth near it, without relaxing the convergence test~~ Equilibrium - BUGFIX! - Recover the magnetic axis when psi is not smooth near it, without relaxing the convergence test Oct 4, 2026

github-actions Bot added the changed-results Results move or an interface breaks - read before upgrading label Oct 4, 2026

d-burg and others added 3 commits October 4, 2026 02:38

Repo - MINOR - Apply the formatter to the axis-bracket change

6d32c85

Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>

d-burg changed the title ~~Equilibrium - BUGFIX! - Recover the magnetic axis when psi is not smooth near it, without relaxing the convergence test~~ Equilibrium - BUGFIX! - ⚠️ Recover the magnetic axis when psi is not smooth near it, without relaxing the convergence test Oct 5, 2026

matt-pharr mentioned this pull request Oct 21, 2026

Equilibrium - REFACTOR - Arrange the Equilibrium module by pipeline stage #551

Merged

Merge 'origin/develop' into bugfix/equilibrium-axis-search-robustness

bf4e298

d-burg added this to the GPEC v2.0.0 milestone Oct 28, 2026

Navigation Menu

Repository navigation

Equilibrium - BUGFIX! - ⚠️ Recover the magnetic axis when psi is not smooth near it, without relaxing the convergence test - #444

Equilibrium - BUGFIX! - ⚠️ Recover the magnetic axis when psi is not smooth near it, without relaxing the convergence test#444
d-burg wants to merge 14 commits into
developfrom
bugfix/equilibrium-axis-search-robustness

d-burg commented Sep 5, 2026 •

edited by matt-pharr

Loading

Uh oh!

github-actions Bot commented Sep 5, 2026

Uh oh!

logan-nc commented Oct 3, 2026

Uh oh!

matt-pharr commented Oct 27, 2026

Uh oh!

This branch has not been deployed

Reviewers

Uh oh!

Assignees

Labels

Projects

Milestone

Development

Uh oh!

3 participants

Conversation

d-burg commented Sep 5, 2026 • edited by matt-pharr Loading Uh oh! There was an error while loading. Please reload this page.

Uh oh!

Release note

Regression report

What changed

Notes for reviewers

Audit round (2026-09-25)

Uh oh!

github-actions Bot commented Sep 5, 2026

Uh oh!

logan-nc commented Oct 3, 2026

Uh oh!

matt-pharr commented Oct 27, 2026

Uh oh!

This branch has not been deployed

Reviewers

Uh oh!

Assignees

Labels

Projects

Milestone

Development

Uh oh!

3 participants

d-burg commented Sep 5, 2026 •

edited by matt-pharr

Loading