Repository navigation
ForceFreeStates - BUGFIX! - ⚠️ Remove the resonant solution at ideal crossings regardless of reduction timing - #489
Merged
Conversation
…sings regardless of reduction timing At an ideal rational-surface crossing the Gaussian reduction ordered columns by growth since the previous reduction and the crossing zeroed the first one. When a ucrit reduction had fired a few steps earlier every column had grown alike, the order was arbitrary, the lead column's pivot fell on a non-resonant harmonic, and the crossing removed the wrong solution. A crossing that landed on the step right after a reduction skipped its own reduction and zeroed a column chosen from the previous one. The ideal crossing now always reduces, leads with the column carrying the most of each resonant harmonic, and pivots it on the resonant row, so that column alone holds the resonant component and is the one zeroed. Kinetic crossings and ucrit reductions are unchanged. Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
… rational-surface crossings ucrit = 2 at eulerlagrange_tolerance 1e-10 makes a ucrit reduction fire two steps before the q = 2 crossing, which on develop leaves the crossing removing a non-resonant solution. Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
Collaborator
|
@d-burg I actually remember what I think is this exact same bug come up in the process of identifying #480. During that, it seemed that this issue had disappeared when using the column wise absolute tolerance. Can you confirm if these results are with that change included, and if not, can you quote the combination of the two fixes and if they are both necessary? |
Collaborator
Author
|
@jhalpern30 my claude answers: the main table in the body is on develop, without the column tolerance. I also ran your #480 head (
So the column tolerance does not remove it. Every wrong run on #480 still zeroes a non-resonant column at a crossing (m = −10, 12, −6, 17, …). At The two fixes do different jobs. #480 controls the error of the small columns; this PR makes the crossing remove the resonant solution whatever the reduction timing. Both are needed, and this branch merges onto #480 without conflicts. To reproduce on your branch: the new |
Collaborator
|
Cool thanks, I'll give this a review at some point in the near future |
…tity Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
d-burg
changed the title
jhalpern30
approved these changes
Oct 6, 2026
jhalpern30
left a comment
jhalpern30
left a comment
Collaborator
There was a problem hiding this comment.
Take what you like, most of my recommendations are for clarity in a section that I think is already confusing and this adds even more logic too.
I'll note that I've found even more bugs in our singular crossing logic and will be documenting them in #481 for further input.
Collaborator
There was a problem hiding this comment.
I am not a fan of our overpopulation of regression harnesses. Is this not a bugfix that can get a targeted unit test for expected behavior? Its a rather niche thing to get an entire regression case, and no one is realistically running with urit = 2 anyway. Is there something covered here that isn't touched by the unit tests below?
|
|
||
| # Fixup solution at singular surface; the ideal jump reduces on the resonant rows so the | ||
| # solution it removes below is the resonant one | ||
| compute_solution_norms!(odet.u, odet, ctrl, intr, true; resonant_rows=(ctrl.kinetic_factor == 0 ? ipert_res : Int[])) |
Collaborator
There was a problem hiding this comment.
This is an ideal crossing - I don't think kinetic factor can ever be nonzero here
| # The reduction above placed the column pivoted on ipert_res[i] at index position i, which also | ||
| # keeps each zeroed column in the same n-block as the resonant mode introduced below; zeroed_idx | ||
| # records the positions for transforming u back to the full solution after integration. | ||
| if ctrl.kinetic_factor == 0 |
Collaborator
There was a problem hiding this comment.
This wasn't changed in this function (and maybe is why the above expression has it too), but I think the same logic applies - this crossing is only called when kinetic factor is > 0.
Collaborator
There was a problem hiding this comment.
Same here- I think remove unless you can see a reason to keep
| compute_solution_norms!(odet.u, odet, ctrl, intr, true) | ||
| singp = intr.sing[ising] | ||
| ipert_res = 1 .+ singp.m .- intr.mlow .+ (singp.n .- intr.nlow) .* intr.mpert | ||
|
|
Collaborator
There was a problem hiding this comment.
For clarity/concistency, maybe define
nres = length(ipert_res)here, and then use for loop bounds below?
Collaborator
There was a problem hiding this comment.
Claude says with your changes, this could be reduced to
odet.u[:, odet.index[i, odet.ifix], :] .= ua[:, ipert_res[i]+intr.numpert_total, :]with the loop for i in 1:nres
|
|
||
| - u: Current solution vector array, updated in-place if fixfac is called | ||
| - sing_flag: Indicates if normalization is occuring at a singular surface or not | ||
| - resonant_rows: Rows of the resonant harmonics at an ideal crossing. When given, the reduction |
Collaborator
There was a problem hiding this comment.
"resonant_rows: Resonant rows at an ideal crossing; when given, the reduction always runs and pivots on them." is shorter and gets the same point across
|
|
||
| # Normalize unorm and perform Gaussian reduction if required | ||
| if odet.new | ||
| if odet.new && isempty(resonant_rows) |
Collaborator
There was a problem hiding this comment.
Completely up to you on this, but I had to read it a few times to understand the control flow here. An option I think is a bit more clear:
# The first call after a reduction records the reference norms; an ideal crossing reduces regardless
at_ideal_crossing = !isempty(resonant_rows)
if odet.new && !at_ideal_crossing
odet.new = false
odet.unorm0 .= odet.unorm
return
end
# Growth since the last reduction (raw norms if a crossing comes right after one), then reduce if required
if !odet.new
odet.unorm ./= odet.unorm0
end
uratio = maximum(odet.unorm) / minimum(odet.unorm)
# Perform Gaussian reduction if ucrit ratio is reached or at singular surface
if uratio > ctrl.ucrit || sing_flag
# …
apply_gaussian_reduction!(u, odet, intr, sing_flag; resonant_rows)
odet.new = true
end| description of `compute_solution_norms!` for more details on the benefits of in-place `u` | ||
| updates. | ||
|
|
||
| Columns are triangularized in order of their growth since the last reduction. At an ideal |
Collaborator
There was a problem hiding this comment.
Shorter option: "At an ideal crossing, the first pivots are the resonant_rows, each on the remaining column largest there. Growth order cannot pick the resonant solution if a ucrit reduction fired just before the crossing, since every column has then grown alike."
| end | ||
| # Sort unorm in descending order (since we triangularize from largest to smallest) | ||
| odet.index[:, ifix] = sortperm(odet.unorm; rev=true) | ||
| # Resonant columns first, then the rest in descending unorm (we triangularize from largest to smallest) |
Collaborator
There was a problem hiding this comment.
Claude's suggested cleanup for this section that I actually kinda like:
# Triangularize primary solutions, resonant rows first, then columns from largest to smallest growth
growth_order = sortperm(odet.unorm; rev=true)
row_free = trues(intr.numpert_total)
col_free = trues(intr.numpert_total)
for isol in 1:intr.numpert_total
if isol <= length(resonant_rows)
kpert = resonant_rows[isol]
ksol = argmax(j -> col_free[j] ? abs(u[kpert, j, 1]) : -Inf, 1:intr.numpert_total)
else
ksol = growth_order[findfirst(j -> col_free[j], growth_order)]
@views kpert = argmax(abs.(u[:, ksol, 1]) .* row_free)
end
odet.index[isol, ifix] = ksol
col_free[ksol] = false
row_free[kpert] = false
# Eliminate other solution vectors below the pivot
for jsol in 1:intr.numpert_total
if col_free[jsol]
# … elimination body unchanged …Note that some of it (i.e. renaming mask to row_free/col_free is outside the scope of your PR but can be done here if you agree with it). This fix put the resonant logic in the same loops and removes some of this complicated lead/filter logic
…g-resonant-elimination
…deal crossings and flatten the reduction control flow The ideal crossing is only reached with kinetic_factor = 0, so its guards never took the other branch. Review suggestions; no change in behavior intended. Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
…duction loop One loop picks each resonant row's column from the columns still free at that step, then continues in growth order. Identical for a single resonant row; with several, later resonant columns are chosen after the earlier eliminations. Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
…tion frequency in the full-run tests Runs the coarse Solovev deck at ucrit = 2, 3, 5 and 10 against its default and replaces the solovev_n1_frequent_reduction harness case. Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
d-burg
marked this pull request as ready for review
October 9, 2026 22:16
Collaborator
Author
|
@jhalpern30 thanks, all taken. Pushed as d309eeb, c9cb55d and 50130ba on a merge of current develop, and the body has the new regression report. From claude: "Harness case: dropped. The end-to-end check is now in runtests_fullruns.jl: the coarse Solovev deck at ucrit = 2, 3, 5, 10 must give the same et[1] as its default (rtol 1e-6). On develop two of the four fail on my machine (13.74 and 16.397 against 16.4037); on this branch all agree to about 1e-9. One thing to look at in the new report: against current develop, solovev_n1 et[1] moves 3.8e-6 relative, where it moved 3.5e-9 against the September develop. The review edits contribute nothing to that, so it is the original fix on post-#480 develop, and it is the size of the scatter in #490." |
d-burg
changed the title
github-actions
Bot
added
the
changed-results
…g-resonant-elimination
Collaborator
I mean I guess, but this still feels like not a good fix. We really need 4 full runthroughs of the code in the unit tests suite to test a relatively niche behavior? Up to you if you think its worth it to change it but I'm not a fan |
Collaborator
Author
|
@jhalpern30 The testset costs about 6 s in CI (five coarse Solovev runs, 23:34:15 → 23:34:21 on the last 1.11 job), and it is the only check that fails on develop for this bug; #520 keeps it with its own Before this merges I'll bring in current develop (#503 moves the vacuum results, so the report needs rerunning). Once that's pushed, could you re-approve the head? Your approval is on |
…g-resonant-elimination
Carries the TOML rename of ucrit to renorm_threshold into the crossing-reduction logic and its tests. Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
Collaborator
|
@d-burg @jhalpern30 — can we wrap this up by next week? It's an approved What's left:
This should land before #520, since #520 keeps this PR's |
Collaborator
Author
|
Merged develop ( |
Collaborator
|
Sure, go ahead and merge |
Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
…g-resonant-elimination
…g-resonant-elimination
Collaborator
|
@d-burg is this waiting on anything? |
matt-pharr
approved these changes
Oct 28, 2026
matt-pharr
left a comment
matt-pharr
left a comment
Collaborator
There was a problem hiding this comment.
Changing this to approved on behalf of Jake's comment
Collaborator
|
@d-burg — completed the requested full SDCC fanout comparison against develop Result: seven cases have changed quantities; harness exit code 0 (no run failures).
These are numerical differences, not an unchanged-results claim; their magnitudes agree with the changes already described in this PR's release note. Environment: Julia 1.12.6, x86_64 Linux, CI manifest pin ./regress_fanout.sh --refs 057c4af0d,a319ed295 --harness-ref a319ed295 --name pr489_a319ed295_vs_057c4af0d_20261028SDCC report: Full regression report (21 cases)Third-party human review and approval are required before merge. |
…g-resonant-elimination
d-burg
enabled stacked auto-merge
October 28, 2026 20:18
…g-resonant-elimination
Collaborator
Author
|
Restamp for the record, since the develop merges after @matt-pharr's run left the stamp behind. |
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
renorm_threshold = 5, et[1] went from 13.74 to 16.4037, its value at the defaultrenorm_threshold(Apple Silicon, develop6fc89d021). The shipped decks move little: Solovev n = 1 et[1] moves by 1.153e-8 absolute (1.06e-6 relative); DIII-D-like headline energies move by up to 1.710e-8 absolute; saved ODE step counts move by up to 0.32%; some kinetic/error-field quantities move by approximately 0.02%. (harness @ a319ed2)In the standard ForceFreeStates (FFS) integration, an ideal rational-surface crossing could remove the wrong solution when a Gaussian reduction had just fired. It did so silently, and the free-boundary energies came out arbitrarily wrong. The crossing now always removes the resonant solution.
Regression report
Latest SDCC fanout rerun (2026-10-28): develop
057c4af0dvsa319ed295— 21 cases; seven cases have numerical changes; no run failures. Full report and limitations. This comparison also includes the existinggal_resistive_pecase with all quantities missing. Earlier reports are retained below.Full SDCC fanout run on 2026-10-28:
regress_fanout.sh --refs 057c4af0dd3c3a834140fad4645b995552944a56,a319ed2955ed3603a75bf074a4a66e4b201d5e46 --name pr489_a319ed2. Baseline: develop057c4af0d; tested head:a319ed295. All 21 cases ran on sirius with Julia 1.12.6 and CI manifest pin4c9d54fa. Cases used 2 Julia threads / 1 BLAS thread, exceptdiiid_error_fieldat 16 / 8; sizing matched between refs. All 42 tasks and report job 1298694 completed with exit code zero.solovev_n1,diiid_n1,diiid_multi_n,diiid_truncation_n1,solovev_kinetic_multiion,solovev_kinetic_ntv,diiid_error_field.gal_resistive_pe: all eight quantities are missing on both refs, as in the previous report; this case provides no numerical validation.The targeted unit/full-run test results documented below are historical results, not reruns at this merge head. Current validation is the full regression harness above.
Full SDCC regression report (21 cases)
What was wrong
At a crossing,
cross_ideal_singular_surf!reduces the solution columns and then zeroes the one it takes to be the large resonant solution. The reduction orders columns by their growth since the previous reduction (unorm ./ unorm0) and zeroes the first.That order carries no information when a
renorm_thresholdreduction fired only a few steps earlier: every column has grown by about 1. The lead column's pivot then lands on an arbitrary harmonic, and the crossing removes a non-resonant solution. A crossing on the very step after a reduction was worse:compute_solution_norms!skipped the crossing's reduction altogether, and the crossing zeroed a column chosen by the previous reduction's ordering.On develop as of
97bd878d3(before #480), Solovev n = 1 withrenorm_threshold= 2–5 hits this for most tolerances. Examples (reference et[1] = 0.6773399618, from Vern9 at reltol 1e-12, abstol 1e-14):renorm_thresholdEvery correct run zeroes m = 2 at q = 2 and m = 3 at q = 3. Across all 18 develop-failing configurations (
renorm_threshold∈ {2, 3, 5} at six tolerances from 1e-7 to 1.4e-10), this PR gives et[1] within 1.8e-7 of the reference at tolerance 1e-7, and within 7.5e-10 at 1e-9 and tighter. At the deck values ofrenorm_threshold(1e3 and 1e4) reductions are rare enough that this was seldom hit, but nothing prevented it.The fix
cross_ideal_singular_surf!: computes the resonant rows before its reduction, and passes them in asresonant_rows.compute_solution_norms!: withresonant_rows, always reduces, including right after arenorm_thresholdreduction.apply_gaussian_reduction!:zeroed_idxrecords their positions (1:nres) fortransform_u!as before. Thekinetic_factor == 0guards in the ideal crossing are removed: that function is only reached on the ideal branch.findfirstpreserved is kept.renorm_thresholdreductions are unchanged. The Riccati path has its own crossing, which zeroes byipert_resdirectly; its harness case is bit-identical.Tests
test/runtests_eulerlagrange.jl, new testset "ideal-crossing reduction isolates the resonant harmonic": a 3 × 3 case in which the resonant component is not in the largest column. It checks three things:The file passes 126/126 and 26/26.
test/runtests_fullruns.jl, new testset "et[1] is independent of the reduction frequency": the coarse Solovev test deck atrenorm_threshold= 2, 3, 5 and 10 must give the same et[1] as at its defaultrenorm_threshold= 1e3, tortol = 1e-6.6fc89d021(Apple Silicon) two of the four fail: 13.74 atrenorm_threshold= 5 and 16.397 atrenorm_threshold= 10, against 16.4037. Which values fail on develop depends on where the steps land, so it can differ by platform.solovev_n1_frequent_reductionharness case from the first version of this PR.Changes since the first review
d309eebf7: deadkinetic_factorguards removed,nres,zeroed_idx = collect(1:nres), shorter docstrings, flattened control flow incompute_solution_norms!.c9cb55d91: single-loop pivot selection withrow_free/col_free. Identical for one resonant row; with several, later resonant columns are chosen after the earlier eliminations. Bit-identical on both multi-n harness cases.50130bafb: the harness case is replaced by the full-run test above.d50a5e5b6: merge of develop after Equilibrium - BUGFIX! - 🚨 Tie the equilibrium ODE absolute tolerance to etol #491; no conflicts.3ef6d15c6: merge of develop after Equilibrium - BUGFIX - Support original and rho grids for inverse equilibria #523, Vacuum/ForceFreeStates - BUGFIX! - Put the vacuum and norm matrices in the plasma's Fourier frame #503, Equilibrium - BUGFIX! - Place read_imas profiles on the stored psi grid #499 and ForceFreeStates - FEATURE! - Cap the integration domain where resistive layers overlap #298; no conflicts.a319ed295: conflict-free merge of develop057c4af0d; the crossing fix and its tests are unchanged. Full SDCC report above.2f436c897: merge of develop9a51446df. One conflict, inEulerLagrange.jl, where develop renameducrittorenorm_threshold: resolved to this PR's structure with the new name. The tests use the new name too. The text of this description uses the new name throughout, including for runs made before the rename.Notes for reviewers
3d6acd813) with no conflicts.renorm_threshold = 10gives a wrong et[1] in 6 of 9 tolerances (1.75, −3.64, −3.84, −0.19, 0.6762, 1.75), where develop gives 0 of 9. Atrenorm_threshold= 2–5 it is wrong in 8 of 18 (develop: 12 of 18).renorm_thresholdthe old choice is usually right.renorm_thresholdfires on norm spread and misses near-parallel columns, giving step-dependent errors of ~1e-5 at the deck values. It is filed as ForceFreeStates - Gaussian reduction trigger misses near-parallel solutions #490. Its temporary remedy, loweringrenorm_threshold, is only safe once this PR lands, since frequent reductions are exactly what trigger the wrong-solution removal fixed here.ode_unorm/ode_fixup. I have not checked whether the Fortran code has the same exposure.🤖 Generated with Claude Code