Repository navigation
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>
|
@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? |
|
@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 |
|
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
left a comment
jhalpern30
left a comment
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.
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[])) |
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 |
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.
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 | ||
|
|
There was a problem hiding this comment.
For clarity/concistency, maybe define
nres = length(ipert_res)here, and then use for loop bounds below?
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 |
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) |
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 |
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) |
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>
|
@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
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 |
|
@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>
|
@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 |
|
Merged develop ( |
|
Sure, go ahead and merge |
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] by 1.6e-8 (1.4e-6 relative), DIII-D-like energies by ≤ 2.3e-8 relative, ODE step counts by ≤ 0.35 %, NTV torques by ≤ 0.01 %. (harness @ 2f436c8)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
regress --refs 9a51446df,2f436c897 --force(all cases) on feynman (SLURM, 8 threads; current develop vs this head with develop merged, both fresh, manifest pinned). The tests inruntests_eulerlagrange.jlandruntests_fullruns.jlpass at this head (126, 26 and 25).solovev_multi_n,diiid_slayer_n1,diiid_n1_riccati,diiid_n1_riccati_precompiled,ggj_reference,ggj_ray_q500i,gal_resistive_diiid,gal_resistive_diiid_galinner,solovev_kinetic_calculated,solovev_kinetic_nuzeroandefit_fixedbdy_separatrix.gal_resistive_pe: extracts nothing on either ref. That is pre-existing on develop and unrelated to this PR.diiid_truncation_n1moves by ≤ 8.8e-9 relative.278057d96/3ef6d15c6,bfc60610e/d50a5e5b6and6fc89d021/50130bafbgave moves of the same size.regress --cases solovev_n1,solovev_multi_n,diiid_n1,diiid_multi_n --refs 1bede6836,d309eebf7and--refs d309eebf7,50130bafb(same job setup, before the Equilibrium - BUGFIX! - 🚨 Tie the equilibrium ODE absolute tolerance to etol #491 merge) compare the merge alone, the merge with the small review edits, and the head with the single-loop pivot selection: every tracked quantity is bit-identical in both comparisons, including both multi-n cases. So all the moves above come from the original fix.solovev_n1diiid_n1diiid_multi_nsolovev_kinetic_ntvsolovev_kinetic_multiiondiiid_error_fielddiiid_truncation_n1What 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.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