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>
This branch has not been deployed
2 participants
Release note
ucrit = 2, et[1] goes from −30.75 to 0.67734. The shipped decks move only at integration-error level: energies by ≤ 3.5e-9 relative, ODE step counts by ≤ 1.3 %, NTV torques by ≤ 0.09 %. (harness @ 7de3a61)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 97bd878d3,7de3a61df --force(all cases) on feynman (SLURM, 8 threads; develop vs head, both fresh, manifest pinned).solovev_n1_frequent_reduction: new in this PR. It is the Solovev deck atucrit = 2,eulerlagrange_tolerance = 1e-10, where develop removes a non-resonant solution at q = 2.diiid_n1_riccati,diiid_slayer_n1,ggj_reference,ggj_ray_q500i,gal_resistive_diiid,solovev_multi_n,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.solovev_n1diiid_n1solovev_kinetic_ntvsolovev_kinetic_multiiondiiid_error_fieldWhat 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
ucritreduction 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, Solovev n = 1 with
ucrit= 2–5 hits this for most tolerances. Examples (reference et[1] = 0.6773399618, from Vern9 at reltol 1e-12, abstol 1e-14):ucritEvery correct run zeroes m = 2 at q = 2 and m = 3 at q = 3. Across all 18 develop-failing configurations (
ucrit∈ {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 ofucrit(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 aucritreduction.apply_gaussian_reduction!:zeroed_idxrecords their positions fortransform_u!as before.findfirstpreserved is kept.ucritreductions 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 107/107.
regression-harness/cases/solovev_n1_frequent_reduction.toml: the end-to-end reproducer above.Notes for reviewers
3d6acd813) with no conflicts.ucrit = 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. Atucrit= 2–5 it is wrong in 8 of 18 (develop: 12 of 18).ucritthe old choice is usually right.ucritfires 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, loweringucrit, 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.NO MERGE WITHOUT HUMAN REVIEW. This PR is a draft and requires a named human reviewer and an assignee before it can be considered for merge. This is non-negotiable.
🤖 Generated with Claude Code