Repository navigation
ForceFreeStates - BUGFIX! - ⚠️ Keep both crossings of a rational surface on reverse-shear profiles - #492
Merged
Conversation
… on reverse-shear profiles The singular-surface finder identified a surface by its q value, so the two crossings of q = m/n on either side of an off-axis q minimum were merged into one entry carrying m = [m, m]: the multi-n merge path intended for (2,1)/(4,2) collisions. The outer crossing then had no chunk boundary, no singfac_min standoff and no ideal jump, the sweep crawled through an unregularized 1/(m - nq) (hundreds of thousands of steps), the Δ′ BVP was skipped as a multi-resonance surface, and the free-boundary energy was wrong by a factor of two against STRIDE on DIII-D 153072 (q₀ = 2.44, q_min = 1.86). Identify surfaces by ψ instead, in both sing_find! and the auto-grid knot list, so a genuine multi-(m, n) collision still merges while distinct crossings stay distinct surfaces with their own ψ and q′ sign. Co-Authored-By: Claude Fable 5.1 <noreply@anthropic.com>
…me unresolvable crossings Deduplicate rational-surface roots on both q and psi, so two different rationals can never merge even where steep q puts them within tolerance in psi. Replace the chunk-bounds assert with an error naming the surface when its singfac_min standoff overruns, which happens for crossings near a q extremum. Add a unit test that swaps a synthetic reverse-shear q into the Solovev fixture and checks both q=2 crossings, their opposite-sign shear, and that (2,1) and (4,2) still merge on each. Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
…or and document the surface merge rule The chunk-bounds error no longer asserts a q extremum or recommends reducing singfac_min; it reports the surface, q', the standoff and the chunk start. sing_find! documents the equal-ratio merge and reverse-shear split in its docstring and drops the Fortran citation. Refresh a stale dedupe comment in bracket_mandatory_nodes and tidy the reverse-shear test's isapprox call. Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
…crossing standoffs Two crossings of q = 2 on an analytic parabola resolve when 2ε > singfac_min and raise the named error when their standoffs overlap; a tangent surface (q' = 0) also raises it. Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
…sed-shear profiles When qmin < qlow < q0 the q < qlow region is a ring around q_min rather than a core, and sing_min! (outer leg) and initialize_el_at_axis! (axis) disagree about the lower integration bound. check_qlow_domain raises a named error at domain setup saying what was detected and to set qlow <= qmin or qlow >= q0. Covered on the analytic reverse-shear fixture in runtests_sing.jl. Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
…kets overlap bracket_mandatory_nodes now errors when c2 - w2 <= c1 + w1 for adjacent brackets, or when a half-width is non-finite (tangent surface), instead of interleaving bracket knots and leaving a surface off-center or sharing a spline interval. The collapse_atol dedupe is kept as a safety net. Geometric condition only; no new tolerance. Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
… with an overlapping standoff The kinetic chunking branch still skips a kinetic singular surface whose crossing standoff reaches back past the chunk start, but now warns with its m/n, psi, q' and standoff. No result changes. Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
…ear-double-rational # Conflicts: # src/ForceFreeStates/Surfaces/Finding.jl # src/GeneralizedPerturbedEquilibrium.jl
github-actions
Bot
added
bugfix
d-burg
changed the title
Collaborator
|
Just attaching Matt here since this seems like a big deal |
Collaborator
|
Merged develop and retained the grid changes in |
…ear-double-rational
… the q < qlow region on the grid itself Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
…to their neighbors and the domain edges instead of erroring Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
…ect-ratio equilibrium for reversed shear Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
…with a double q = 2 crossing Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
…the run reports it Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
…ear-double-rational
d-burg
marked this pull request as ready for review
October 29, 2026 04:29
logan-nc
reviewed
Oct 29, 2026
logan-nc
left a comment
logan-nc
left a comment
Collaborator
There was a problem hiding this comment.
Review in progress. Posting now so nobody duplicates this work.
The core fix looks right. Matching surfaces on (q, ψ) keeps both crossings while (2,1)/(4,2) still merge. Both roots come from Brent on the same q-extremum bracket, so they agree far inside the 1e-8 tolerance. The fix also matters on default runs: every shipped deck has qlow = 1.02, below any realistic reversed-shear q_min, so both crossings are in the domain.
In progress, on a local branch (not pushed yet): replacing check_qlow_domain with a single qlow lower bound. Right now three rules disagree:
initialize_el_at_axis!(forward and Riccati) raises the start only whenqlow > q0;sing_min!usesqlow > params.qminand only runs for Galerkin;remove_singular_surfs!filters by q, not ψ.
In the window q_min < qlow ≤ q0 the integration therefore starts at the axis and runs through surfaces that have already been filtered out. The guard errors on that window; the change below removes it.
qlow_start_psiis the outermost crossing of q = qlow, taken from the same spline q extremasing_find!uses and solved with a bracketed root. All integrators start there, andsing_min!runs for every integrator.remove_singular_surfs!drops surfaces by ψ against that bound.check_qlow_domainand its tests are deleted.- This also fixes a latent bug:
params.qminis extrapolated to ψ = 0, so the oldsing_min!Newton solve could start at the first node and converge to a root outside the grid. - On monotonic profiles the start and the filter match develop's, so shipped decks should not move; the harness will confirm. A
qlowscan on the reversed-shear LAR fixture is running and I'll post the table here.
Smaller items, also planned:
- merge develop (it merges cleanly);
maxlogon the bracket-cap warning, which fires per surface on each grid pass;- pin the sign of q′ on each crossing in the new case or a test (the case currently tracks only ψ and q, and its
mpsi = 128bypasses the two-pass grid); - name the remedy (lower
psilowor raiseqlow) in the crossing-standoff error; - re-run the full harness at the merged head. The current stamp
437387af6predates36cbc5e31(bracket cap) and538dcf569(grid-read qlow check), so "no impact on shipped decks" is not yet shown.
Follow-up, not in this PR: downstream code on negative-shear surfaces (resistive E/F/H with q1² and q1³, Δ′, PerturbedEquilibrium coupling, SLAYER) is not exercised; the new case is ideal-only. I'll open an issue.
@d-burg, OK if I push these to your branch, or would you rather take them? @matt-pharr, please hold your review until the push.
…rom the spline q extrema and filter surfaces by psi Replaces check_qlow_domain: qlow_start_psi (outermost crossing of q = qlow, bracketed on the spline q extrema sing_find! uses) is now the lower bound for forward, Riccati and Galerkin alike, and remove_singular_surfs! drops surfaces by psi against it. A qlow between q_min and q0 now excludes the reversed-shear ring and the core instead of starting at the axis or erroring. Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com> Claude-Session: https://claude.ai/code/session_01VdRsGGg4YRG6E9fgLGbpdF
…nd pin dq/dpsi in the reversed-shear case Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com> Claude-Session: https://claude.ai/code/session_01VdRsGGg4YRG6E9fgLGbpdF
…ear-double-rational
Collaborator
|
Pushed
The touched unit tests pass (676). A full 24-case harness run at |
logan-nc
approved these changes
Oct 29, 2026
logan-nc
left a comment
logan-nc
left a comment
Collaborator
There was a problem hiding this comment.
Approving. The reversed-shear fix is correct: both crossings are kept as separate surfaces, and (2,1)/(4,2) still merge. The qlow handling now uses one consistent lower bound for every integrator, confirmed by the qlow scan. The stamp still needs updating from the full harness run in progress before merge.
logan-nc
enabled auto-merge
October 29, 2026 18:41
Collaborator
|
Full 24-case regression report, develop Report |
Collaborator
|
Full 24-case regression report, develop Report |
This file contains hidden or bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment
Add this suggestion to a batch that can be applied as a single commit.This suggestion is invalid because no changes were made to the code.Suggestions cannot be applied while the pull request is closed.Suggestions cannot be applied while viewing a subset of changes.Only one suggestion per line can be applied in a batch.Add this suggestion to a batch that can be applied as a single commit.Applying suggestions on deleted lines is not supported.You must change the existing code in this line in order to create a valid suggestion.Outdated suggestions cannot be applied.This suggestion has been applied or marked resolved.Suggestions cannot be applied from pending reviews.Suggestions cannot be applied on multi-line comments.Suggestions cannot be applied while the pull request is queued to merge.Suggestion cannot be applied right now. Please check back later.
Release note
diiid_truncation_n1, where the edge bracket cap moves the outermost layer-overlap width by 0.13 % and the energies by about 1e-8. Reversed-shear equilibria now give correct results. (harness @ b0f1b06)qlownow sets one lower bound for every integrator: integration starts at the outermost crossing of q = qlow and excludes the q < qlow region inside it, whether a core or, on reversed shear, the ring around q_min and the core within it. On monotonic profiles this is the same start as before. Previously, withq_min < qlow ≤ q0, the forward and Riccati paths started at the axis after the q < qlow surfaces had already been filtered out. Rational surfaces too close for their usual grid brackets get brackets capped at half the gap, with a warning, instead of off-center brackets.On a reversed-shear profile, the two crossings of q = m/n on either side of an off-axis q minimum were merged into one surface. The outer crossing then had no chunk boundary, no standoff and no ideal jump, and the free-boundary energy was wrong by a factor of two against STRIDE on a DIII-D reversed-shear case. Surfaces are now identified by ψ, so distinct crossings stay distinct while genuine multi-(m, n) collisions still merge.
Regression report
Full harness at
b0f1b0685(developb6a28892dmerged in) against developb6a28892d, all 24 cases, on feynman: Julia 1.11.7, CI-pinned manifest, 8 threads / 8 BLAS per case.diiid_n1(53 quantities), every SLAYER deck, the Riccati and Galerkin decks, and the kinetic Solov'ev decks.lar_reversed_shear_n1(new in this PR): develop cannot run it (no example directory there). At this head both q = 2 crossings are separate singular surfaces with opposite-sign dq/dψ (now pinned), et[1] = −0.0226 with no wall, and the forward integration takes 227 saved steps.gal_resistive_pe: eight N/A rows on both refs; pre-existing.diiid_truncation_n1: 3 changed, expected. Its outermost layer-overlap surface sits at ψ = 0.99975, just inside the domain edge. The new edge cap narrows that surface's grid bracket from the 1e-4 stencil to 7.3e-5 (the run logs the "capped" warning for exactly this surface), so the edge grid differs slightly. Only the two outermost of the seven diffusive-resistive widths move (3e-7 relative and 0.13 %); the inner five agree to round-off.ep[1]andet[1]move by about 1e-8 relative andpsilim_overlapby 1e-10. Singular-surface locations, q,psilim, and the overlap index are unchanged.runtests_sing(90),_eulerlagrange(133, 26, 38, 2),_riccati(14),_galerkin(3),_grid_refinement(84),_solve_api(144), and_parallel_integration(3, 82) pass.runtests_equilpasses 337 of 338: oneLoad TokaMaker i-fileassertion compares re-formed and freshly read splines with exact==. It fails the same way on developb6a28892don alfven and passes in CI, so it is a machine-dependent floating-point equality outside this PR.qlowscan onlar_reversed_shear_n1(q₀ = 2.35 on the grid, q_min = 1.82). The "before" column is this branch with the error check removed, which is develop's behaviour on this fixture:diiid_error_fielddiiid_multi_ndiiid_n1_riccati_precompileddiiid_n1_riccatidiiid_n1diiid_slayer_n1_chi_perp_ediiid_slayer_n1_coupleddiiid_slayer_n1diiid_slayer_n1_toroidaldiiid_slayer_n1_windingdiiid_truncation_n1efit_fixedbdy_separatrixgal_resistive_diiid_galinnergal_resistive_diiidgal_resistive_peggj_ray_q500iggj_referencelar_reversed_shear_n1solovev_kinetic_calculatedsolovev_kinetic_multiionsolovev_kinetic_ntvsolovev_kinetic_nuzerosolovev_multi_nsolovev_n1diiid_truncation_n1reportlar_reversed_shear_n1reportThe complete 24-case report is in a comment below (it exceeds the description size limit).
Changes since the review
36cbc5e31). An overlap error was new in this PR and also fired on ordinary profiles develop runs today, e.g. several n on a steep edge, and on rationals outside the integration domain, since the grid brackets every rational up topsihigh. A near-tangent rational could also produce a huge but finite bracket that cleared most of the grid without erroring. Each bracket is now capped at half the gap to its neighboring centers, and its stencil part also at half the distance to the nearer domain edge, so every center keeps its own interval and no bracket can clear the grid. A cap below the Δ′ stencil warns. Cases with non-overlapping brackets are unchanged.psilimtakes the outer crossing (538dcf569).sing_lim!seeded its Newton solve at the node nearest qlim, which on reversed shear can be the inner leg; it now seeds from the edge inward at the outermost node with q ≤ qlim.qlowbound for every integrator (e200765c2). This replacescheck_qlow_domain. The domain had three lower-bound rules:initialize_el_at_axis!raised its start only whenqlow > q0,sing_min!used the extrapolatedparams.qminand ran only for Galerkin, andremove_singular_surfs!filtered by q. Nowqlow_start_psigives the outermost crossing of q = qlow, solved with Brent on the same spline q extremasing_find!uses, so a dip between grid nodes counts and the root cannot leave the grid.sing_min!sets it for every integrator, and surfaces are filtered by ψ against it. The crossing-standoff error now names the remedies.6a6954483). The bracket cap now warns once per grid pass, listing the capped surfaces, instead of once per surface.lar_reversed_shear_n1also tracksSingularSurfaces/dqdpsi, the opposite-sign q′ of the two crossings.a97fe4aa2):sigma_type = "hollow", σ = σ₀(1 + h x²)(1 − x²)^p, which keeps q(0) = q0 and gives reversed shear for large enoughhollow. The three copies of the σ formula are now one helper. The new case and an equilibrium test use it (ca75806a9).Notes for reviewers
b0f1b0685(developb6a28892d, including InnerLayer.SLAYER - FEATURE! - 🚨 Derive the toroidal critical-Δ geometric factor in the r_s reference #441 and Repo - API! - Read kinetic profiles once into the equilibrium and derive every kinetic view from them #581), with no conflicts.🤖 Generated with Claude Code