Repository navigation
PerturbedEquilibrium - BUGFIX! - Converge the singular coupling with singfac_min = 1e-6 in the forward-integrator decks - #555
Merged
Conversation
…the singular coupling converges At 1e-4 the coupling at the outermost rational surfaces is off by up to 20% and differs between working coordinates; at 1e-6 Hamada, PEST and Boozer agree to 0.5% at no measurable runtime cost. Example decks follow the new default; the test fixtures and Fortran benchmark namelists keep 1e-4. Documents the row and column normalization of the singular coupling matrices. Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com> Claude-Session: https://claude.ai/code/session_01VdRsGGg4YRG6E9fgLGbpdF
…singfac_min = 1e-6 in the forward-integrator decks The Riccati Δ' gains a spurious imaginary part as singfac_min shrinks (on the DIII-D-like deck, surface 4 goes from -2275-9i at 1e-4 to -2278-421i at 1e-5 and -2449-19369i at 1e-6), so the code default and the Riccati, Galerkin and SLAYER decks stay at 1e-4. Only the forward-integrator decks, which feed PerturbedEquilibrium, move to 1e-6, and the singfac_min docstring states the trade-off. Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com> Claude-Session: https://claude.ai/code/session_01VdRsGGg4YRG6E9fgLGbpdF
github-actions
Bot
added
bugfix
This was referenced Oct 22, 2026
d-burg
approved these changes
Oct 26, 2026
d-burg
left a comment
d-burg
left a comment
Collaborator
There was a problem hiding this comment.
Thanks for chasing #237 down to the offset. Showing that the coupling definitions are sound and the scatter lives in the shared solution is useful, and so is the write up of the Riccati trade off. This only moves the jump closer to the surface and loosens nothing, so I'm good with it. Before it goes in, some Claude findings for your consideration:
- Scope: any forward run is affected, not only PE. Your Solovev numbers show the FFS energy
etmoving 32 % (NTV 119 %) between 1e-4 and 1e-6. Could the docstring and Migration line say "forward runs", and could a forward run that keeps 1e-4 warn (or get an integrator-aware default)? Related:diiid_n1_riccatiat 1e-4 giveset= 0.8032 against 0.7979 for converged forward on the same deck, so the Riccati energies look offset-limited too. That seems worth a follow-up issue next to the crossing-match question. - Convergence evidence for the coupling. The Solovev
etscan is convincing. For the DIII-D-like coupling we only see 1e-4 and 1e-6 plus the cross-coordinate check, and the outer surface is still 1.4 % apart. Could you add a 1e-7 point for resonant b^r, island widths and Chirikov ondiiid_n1, and a tighter-rtol control for the Solovev scan (those decks run at rtol 1e-7, so the 1e-8 breakdown may be tolerance-driven)? A small harness case pinning the Hamada/PEST/Boozer agreement would lock in the #237 fix. - Kinetic decks. With
kinetic_factor > 0,singfac_minsets the trapezoidal step across the kinetic singular surfaces, which the evidence here doesn't cover. I'd revert those two decks (suggestions inline) or add a scan for that path.
Smaller items:
- The coupling is read at a fixed 5e-4/(n·q′) offset that doesn't follow
singfac_min; a sentence on why that spot is converged would help (comment fix inline). - The truncation and a10 decks stop before PE, so the PE rationale doesn't cover them; the forward-
etrationale does. - On auto-grid decks
singfac_minalso sets the knot bracketing (covered in suggestion 1). - "Runtime under 4 %" is harness wall time, dominated by compilation; ODE steps go up 35–50 % (
diiid_n1781 → 1174). - The Numerical impact line is missing the kinetic ξ norms (7.8 %), the error-field residual torque (up to 4.8 %) and the truncation
et(0.34 %). - The test fixtures still pin the 1e-4 values. Fine to defer, but worth a tracking issue.
The ! and the regression report look right to me.
| - `numsteps_init::Int` - Initial array size for ODE data storage | ||
| - `numunorms_init::Int` - Initial array size for solution normalization data | ||
| - `singfac_min::Float64` - Fractional distance from rational q at which ideal jump condition is enforced | ||
| - `singfac_min::Float64` - Fractional distance from rational q at which ideal jump condition is enforced. The default 1e-4 suits the Riccati Δ′, whose asymptotic matching gains a spurious imaginary part as it shrinks; `"forward"` runs feeding PerturbedEquilibrium need about 1e-6 for the singular coupling to converge at the outer rational surfaces. |
Collaborator
There was a problem hiding this comment.
Scope the docstring to forward runs in general, mark the Riccati cause as observed, and name the two other places the value acts.
Suggested change
| - `singfac_min::Float64` - Fractional distance from rational q at which ideal jump condition is enforced. The default 1e-4 suits the Riccati Δ′, whose asymptotic matching gains a spurious imaginary part as it shrinks; `"forward"` runs feeding PerturbedEquilibrium need about 1e-6 for the singular coupling to converge at the outer rational surfaces. | |
| - `singfac_min::Float64` - Fractional distance from rational q at which ideal jump condition is enforced. With `kinetic_factor > 0` it instead sets the half-width of the trapezoidal step across each kinetic singular surface, and on the auto ψ grid it also sets the rational-surface knot bracket. `"forward"` runs need about 1e-6 for the free-boundary energies and the PerturbedEquilibrium singular coupling to converge. The default 1e-4 is kept for the Riccati Δ′, whose imaginary part grows as the offset shrinks (observed, cause not yet isolated), so Riccati energies at the default are not offset-converged. |
Comment on lines
+258
to
+260
| - Columns: every stored matrix is right-multiplied by the control-surface `flux_conform`, so it | ||
| acts on the root-area-weighted field b̃ in T and its singular values do not depend on `jac_type`. | ||
| The applied vectors are evaluated before that, on the run's own forcing flux. |
Collaborator
There was a problem hiding this comment.
Say which offset the coupling is actually read at, and limit the jac_type claim to the converged case.
Suggested change
| - Columns: every stored matrix is right-multiplied by the control-surface `flux_conform`, so it | |
| acts on the root-area-weighted field b̃ in T and its singular values do not depend on `jac_type`. | |
| The applied vectors are evaluated before that, on the run's own forcing flux. | |
| - Columns: every stored matrix is right-multiplied by the control-surface `flux_conform`, so it | |
| acts on the root-area-weighted field b̃ in T and, once the solve is converged, its singular values do not depend on `jac_type`. | |
| The applied vectors are evaluated before that, on the run's own forcing flux. | |
| - Offset: the jump in ∂b^ψ/∂ψ is read at ψ_s ± 5e-4/(|n|·|q′|), independent of `singfac_min`. |
| integrator = "forward" # Integration formalism: "forward" (dense ξ, needed for [PerturbedEquilibrium]/kinetic), "riccati" (chunked propagator BVP, unlocks the singular-surface Δ′ matrix), or "galerkin" (RDCON outer-region singular Galerkin Δ′, with the optional RPEC inner-layer match) | ||
| eulerlagrange_tolerance = 1e-7 # Relative tolerance for ODE integration of Euler-Lagrange equations | ||
| singfac_min = 1e-4 # Fractional distance from rational q at which ideal jump enforced | ||
| singfac_min = 1e-6 # Fractional distance from rational q at which ideal jump enforced |
Collaborator
There was a problem hiding this comment.
With kinetic_factor > 0 this sets the kinetic crossing step, not the ideal jump, and the 7.8 % ξ-norm move is unexplained. Revert unless a kinetic-crossing scan is added.
Suggested change
| singfac_min = 1e-6 # Fractional distance from rational q at which ideal jump enforced | |
| singfac_min = 1e-4 # Fractional distance from rational q at which ideal jump enforced |
| integrator = "forward" # Integration formalism: "forward" (dense ξ, needed for [PerturbedEquilibrium]/kinetic), "riccati" (chunked propagator BVP, unlocks the singular-surface Δ′ matrix), or "galerkin" (RDCON outer-region singular Galerkin Δ′, with the optional RPEC inner-layer match) | ||
| eulerlagrange_tolerance = 1e-7 # Relative tolerance for ODE integration of Euler-Lagrange equations | ||
| singfac_min = 1e-4 # Fractional distance from rational q at which ideal jump enforced | ||
| singfac_min = 1e-6 # Fractional distance from rational q at which ideal jump enforced |
Collaborator
There was a problem hiding this comment.
Same mechanism, this deck also stops before PE and has no harness case.
Suggested change
| singfac_min = 1e-6 # Fractional distance from rational q at which ideal jump enforced | |
| singfac_min = 1e-4 # Fractional distance from rational q at which ideal jump enforced |
| shear = abs(nn) * sing_surf.q1 / sing_surf.q | ||
|
|
||
| # Evaluate bwp1_mn = ∂b^ψ/∂ψ at lpsi and rpsi using permeability-weighted eigenstates. | ||
| # Matches Fortran gpout_resp: evaluate bwp1_mn at lpsi/rpsi via gpeq_sol |
Collaborator
There was a problem hiding this comment.
The "matches Fortran default singfac_min" comment no longer matches either value in use. Suggested one-liner:
# Fixed evaluation offset 5e-4/(|n|·|q′|), independent of singfac_min
This was referenced Oct 27, 2026
Collaborator
|
Opened #610 to track further Solovev cutoff-convergence investigation after the 2.0 release, alongside the crossing/conditioning work in #481 and #520. We will not hold this PR up on that additional numerical investigation. The follow-up distinguishes the FFS integration/crossing cutoff from PE's coupling-sampling offset and its field regularization. The ideal-plus-NTV examples use a prescribed regularization; their torque is a downstream diagnostic of the ideal response, not a self-consistent kinetic torque. The calculated-kinetic examples need separate validation. The remaining task is to establish a supported cutoff range through tolerance/resolution checks and focused regression coverage, rather than assume smaller cutoffs are always safer. This does not call for reverting the example changes or adding a solver refactor to #555. The existing human approval by @d-burg is noted. GitHub currently reports merge conflicts, so those still need resolution and the normal review/check requirements still apply before merging. I will do that and get this merged now. |
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.
Closes #237.
Review package: https://claude.ai/artifact/VMcqrrHhyxbnJiMFzKQzLz (figures for every claim below; ask @logan-nc for access if the link doesn't open)
Release note
"forward"into PerturbedEquilibrium should setsingfac_min = 1e-6; the code default stays 1e-4.The resonant coupling matrices (resonant field, current, island width², Δ′) were found not to depend on the working coordinates once the solve is converged. The 8–20% scatter reported in #237 came from the ideal-jump offset
singfac_min = 1e-4, which left the outer rational surfaces unconverged. The forward-integrator example decks now use 1e-6. The code default stays 1e-4, because the Riccati Δ′ degrades at smaller offsets.What was found
singfac_min = 1e-6, Hamada, PEST and Boozer agree on |bʳ| to within 0.5% on four of the five DIII-D-like rational surfaces (1.4% for PEST at ψ = 0.993), and on the singular values of C̃ to within 0.7%. At 1e-4, the q = 4 surface was off by 12% (Hamada) to 20% (PEST, Boozer).mtheta, ODE tolerance, coil forcing (‖b̃ₓ‖ agrees to 0.03%) and the plasma response (leading permeability eigenvalue agrees to about 1%).C_penetrated_area_weighted_fieldcannot be checked yet: runs that match inner-layer solutions skip the singular coupling, so no deck writes it. Its rows are a field in T by construction, now documented.Why only the forward decks
The same offset sets where the Riccati integrator matches to its asymptotic series across each rational surface. Its Δ′ diagonal should be real, but |Im Δ′|/|Re Δ′| on the q = 5 surface grows from 0.4% at 1e-4 to 18% at 1e-5 and 790% at 1e-6. A first version of this branch set 1e-6 globally, and the harness showed the Riccati Δ′ move by 851%. So:
"forward"decks feeding PE now use 1e-6: DIIID-likeerror_field,ideal,ideal_multi_n,truncation; Solovevideal,kinetic_NTV,kinetic_calculated,kinetic_multiion;a10_kinetic.CoreTypes.jldocstring now states the trade-off.compute_singular_coupling_metrics!documents the row and column normalization.Notes for reviewers
sing_get_uainCrossings.jl) loses accuracy. That is my inference and is out of scope here.test/test_data/regression_*) keep 1e-4, so no golden values move. The full test suite passed on the first (global 1e-6) version of the branch; CI covers this one.Regression report
Updated after merging develop: 21 cases / 42 successful runs, zero failed runs, comparing develop
f43b5bee8f1dc1bc772923b86f029c9465f0e3b0with PR head4ecda73be5dd823a65e45dcf505380dd17d2ca10. Harness exit status: 0. Julia 1.12.6, pinned manifest4c9d54fa, 8 Julia threads / 4 BLAS threads on both refs.The cutoff-induced numerical changes remain visible: ideal Solovev NTV torque is 39.47649088015491 → 86.57098121746682 N·m (+119.30%); multi-ion torque changes by +119.18%. The calculated-kinetic Solovev NTV quadrature torque changes by approximately +0.0064%, while displacement norms change by 8.35%. DIII-D-like n=1 resonant-field components change by up to 1.98%, island widths by 1.76%, and Chirikov by 5.32%; the spurious ideal PE torque drops by 99.19%.
Coverage caveat: all eight non-runtime quantities in
gal_resistive_peare missing/N/A on both refs. That case completed but provides no numerical validation of those outputs. Riccati, SLAYER, GGJ, the two Galerkin stability cases, and Solovev multi-n have unchanged tracked non-runtime quantities. Further cutoff-convergence investigation is deferred until after 2.0 in #610; this regression comparison is not a convergence proof.Full 21-case regression report (spacing compacted; values unchanged)
🤖 Generated with Claude Code
https://claude.ai/code/session_01VdRsGGg4YRG6E9fgLGbpdF