Repository navigation
ForceFreeStates - BUGFIX! - Trigger the Gaussian reduction on the conditioning of U₁ and factor U₁ in the plasma response matrix - #520
Conversation
…ditioning of U₁ and factor U₁ in the plasma response matrix The forward integrator reduced only when the U₁ column-norm spread exceeded ucrit, which misses columns that grow alike and turn near-parallel. It now reduces when the 1-norm condition estimate of the column-normalized U₁ exceeds ucrit (ucrit keeps its entry-size meaning in the Riccati path). The reduction storage grows on demand instead of overflowing. wp = U₂ / U₁ factored U₁ᵀ, so partial pivoting ran across U₁'s columns, whose norms span ~1e110 right after a reduction; it now factors U₁ itself. Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
github-actions
Bot
added
bugfix
…reduction check Drop the column_condition helper and its linear-algebra unit tests. Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
…reduction check Drop the column_condition helper and its linear-algebra unit tests. Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
…m:OpenFUSIONToolkit/GPEC into bugfix/ffs-reduction-conditioning-trigger
|
cc-ing @logan-nc @matt-pharr since this would basically fundamentally change the decades-old approach of how fixups are done in ForceFreeStates. After a lot of investigation/tests, it seems that the growth oriented reduction does not properly capture when a fixup is needed. Instead, the condition number of the displacement columns is more accurate, producing good numerics. Specifically, this fixes the Solov'ev case, which has a higher alpha (Frobenius coefficient in the expansion) at the singular surfaces, and was giving garbage as And on a related note, I've verified and tested that the 3D code also has this same issue (but worse) and is fixed by this method |
|
I'm happy to implement tested improvemlike this! this! The fixup is the magic sauce, so having our own recipe is a good result of the rewrite effort |
|
@jhalpern30 be sure to post some plots highlighting the results, a table is nice but a plot is a much better record! |
|
@matt-pharr figures are in the artifact linked in in the description |
d-burg
left a comment
d-burg
left a comment
There was a problem hiding this comment.
The mechanism is convincing to me: columns that grow alike and turn parallel are invisible to a spread test, and the new unit tests separate the two triggers nicely. The storage overflow and the U₂ / lu(U₁) fix are clear bug fixes. I'm less sure the evidence is there yet for replacing the reduction trigger itself.
My and claude's brains merged into a disconcerting synthetic hybrid network and yielded the following items:
Blocking
- Migration. The title carries
!and Migration starts with "none", which is what fails the title check, and substantively a forward-path deck with a hand-setucritis reinterpreted without a warning. Suggested text: "Migration: forward-integrator decks that setucritby hand should return it to the default 1e3. The value is now compared with a condition estimate of the column-normalized U₁, not with the spread of column-norm growth; small values such as 10 now reduce very often, and values above about 1e10 lose digits. Riccati decks need no change." - Regression report. The code block is empty and the stamp (
1dc6138aa) predates the last two commits. This changes when every forward run reduces and how W is formed for every integrator, so I'd like the full harness now,d50a5e5b6against the current head as two git refs, especially the Riccati, perturbed-equilibrium, kinetic Solovev and multi-n cases. Rerunning at merge time is cheap; reviewing blind is not.
Evidence I'd want for "changes the decades-old approach"
- Every error in the table is measured against the same build with
ucrit= 1, which shows self-consistency, not correctness. One independent reference per headline row would help: the Riccati integrator on the same Solovev decks atsingfac_min1e-8, or a tolerance-converged or Fortran value where available. The figures sit behind a link "shared on request"; the numbers the argument rests on should be in the PR. - Per case: accepted steps, number of reductions, the largest estimate at the moment of reduction, wall time and peak memory against #489's head. Without reduction counts nobody can judge 1e3 against 1e6, and with unbounded storage a run that reduces every other step is a memory problem the old 100-slot cap used to hide.
- Cases not yet covered: an ideal-unstable equilibrium where det U₁ passes through zero inside the plasma (U₁ is singular there for a physical reason, so the trigger fires whatever its threshold), a kinetic run, and the largest multi-n deck.
On your five decisions
- 1e3 vs 1e6: undecided until the reduction counts are shown; leaning 1e6 for the forward path. The stated reason for 1e3 is about Riccati sharing the key, not about the forward path.
- One key, two meanings: I'd add a forward-path key. Decks are reinterpreted silently, switching
integratorchanges what the same number means, and the shared key is already forcing a default neither path would pick on its own. The!in the title makes a new key cheap now. - cond(U₁) rather than cond([U₁; U₂]): agree; the elimination pivots on U₁ and W inverts U₁.
- Elimination rather than an (S, I) or QR reset: agree for this PR. Worth a follow-up issue: the elimination never rescales columns, which is why the norms span 1e110 in the first place.
- No test for the
lufix: disagree; a 4×4 matrix reproduces it (inline).
Suggestion: the storage fix and the lu fix stand on their own and are easy to verify; the trigger change is a design decision that needs the evidence above. Splitting them would let the two bug fixes go in quickly on #489.
Small, outside the diff
- The
OdeStatefield docstrings inEulerLagrange.jlstill give the reduction arrays fixed shapes ((numpert_total, numunorms_init), "lengthnumunorms_init"); they are initial capacities now.docs/src/stability.mdstill callsnumunorms_inita "renorm checkpoint budget". - The comment at
Riccati/Driver.jl:77still quoteswp = u[:,:,2] / u[:,:,1]. - The description's notes refer to
column_condition, which the last two commits inlined.
| Increase `numunorms_init` if needed. Automatic resizing will be added in a future version." | ||
| end | ||
| # Reduce at a singular surface or once the U₁ columns lose independence (1-norm condition estimate of the column-normalized U₁) | ||
| @views if sing_flag || cond(UpperTriangular(qr!(u[:, :, 1] ./ norm.(eachcol(u[:, :, 1]))').R), 1) > ctrl.ucrit |
There was a problem hiding this comment.
Two claude things about this trigger.
A 1-norm condition estimate of a column-normalised matrix has a floor that grows with the matrix size. Applying this elimination to random complex matrices leaves the estimate around 3e2 at N = 34 and above 1e3 from about N = 64, so for large numpert_total (multi-n, high n with a wide m range) the default 1e3 could fire on every second step. Real solution matrices are structured and may sit well below that, so this is a question: could you post odet.ifix and odet.step for each harness case, in particular the multi-n deck and DIIID n = 6/8/10? If ifix is close to step/2 anywhere, I'd compare against the value right after the last reduction, or scale the threshold with N, rather than raise the default.
It also does a QR plus about five N×N allocations on every accepted step, where the old test reused norms already computed. Standalone that is 0.05 ms at N = 34 and 7.7 ms at N = 400; probably small next to the right-hand side, but could you post benchmarks/benchmark_git_branches.jl numbers for the DIIID ideal and multi-n decks against #489's head? I'd also bring back the named helper the description mentions: it is the natural place for a reused work matrix and a direct unit test of the estimate.
| # Compute plasma response matrix W = U₂ * U₁⁻¹ | ||
| wp = zeros(ComplexF64, numpert_total, numpert_total) | ||
| @views wp .= (odet.u[:, :, 2] / odet.u[:, :, 1]) ./ equil.psio^2 | ||
| @views wp .= (odet.u[:, :, 2] / lu(odet.u[:, :, 1])) ./ equil.psio^2 |
There was a problem hiding this comment.
The claim checks out, and it is reproducible at 4×4: take L[i,j] = 1 on the diagonal and 0.8cis(i + 2j) below it, U1 = (L .* exp10.(range(110, 0; length=4))')[4:-1:1, :], A = [cis(0.3i*j)/(i+j) for i in 1:4, j in 1:4], W = A + A', and U2 = W * U1. Then U2 / U1 is off by 0.15 relative and U2 / lu(U1) by 2e-16 (other Hermitian W give 0.1 to 1). With random dense U₁ of the same column scales both forms are fine, which is probably why the search came up empty.
Since the expression now appears three times (here, line 279 and Result.jl:243), I'd put it in one small helper used by free_run, free_compute_total and build_result, and unit-test it with that matrix at rtol = 1e-12.
This path is shared with Riccati, so "the Riccati path is untouched" should be backed by the Riccati rows of the harness table. I expect rounding-level moves only.
| - `ksing::Int` - Singular surface handling parameter | ||
| - `eulerlagrange_tolerance::Float64` - Relative tolerance for ODE integration of Euler-Lagrange equations | ||
| - `ucrit::Float64` - Critical value of unorm ratio to trigger solution normalization. In the standard path it triggers Gaussian reduction; in the Riccati path it triggers `renormalize_riccati_inplace!`. Default `1e3`. The trigger compares column-norm growth, not conditioning, so near-parallel columns can leave an error in the free-boundary energies that does not shrink with `eulerlagrange_tolerance`. | ||
| - `ucrit::Float64` - Normalization threshold, default `1e3`. Forward path: the condition number of the column-normalized U₁ above which a Gaussian reduction is applied; the reduction pivots on U₁, so its rounding error and that of the reconstructed `u_store` scale as eps·`ucrit` (values ≥ 1e10 lose digits). Riccati path: the largest entry of U₁ or U₂ above which `renormalize_riccati_inplace!` runs. |
There was a problem hiding this comment.
The docstring says the rounding error scales as eps·ucrit, but the check runs once per accepted step and is skipped on the first step after a reduction, so the conditioning when the elimination actually runs is at least ucrit and can be well above it where growth per step is large, which is the small-singfac_min regime this PR targets. Do you have the distribution of the estimate at the moment each reduction fires for the Solovev singfac_min = 1e-8 case? If the overshoot is routinely a few decades, I'd say "at least", and also that this is a 1-norm estimate whose healthy baseline grows with the number of modes.
| Doubles the capacity of the per-reduction arrays (`fixfac`, `index`, `sing_flag`, `fixstep`, `zeroed_idx`). | ||
| """ | ||
| function grow_reduction_storage!(odet::OdeState) | ||
| nfix = size(odet.fixfac, 3) |
There was a problem hiding this comment.
Doubling zero is still zero, so numunorms_init = 0 throws a BoundsError on the first reduction now that the docstring says the storage grows as needed.
| nfix = size(odet.fixfac, 3) | |
| nfix = max(size(odet.fixfac, 3), 1) |
| end | ||
| et_ref = total_energy(1e3) | ||
| for ucrit in (2.0, 3.0, 5.0, 10.0) | ||
| for ucrit in (3.0, 10.0, 1e2, 1e6) |
There was a problem hiding this comment.
Why were 2.0 and 5.0 dropped from this ladder? The rtol is unchanged, but removing the smallest values narrows what the test accepts, and 5.0 was the value that fails on develop for the #489 bug. The reference for your whole results table is the ucrit = 1 run, and the description claims insensitivity from 10 to 1e8. If 2.0 still passes at 1e-6, please keep it and add 1.0 and 1e8 so the test covers the claimed range; if it doesn't, we should understand that before trusting the reference.
| end | ||
|
|
||
| # A ucrit reduction landing just before an ideal crossing must not change which solution the crossing removes. | ||
| # A reduction landing just before an ideal crossing must not change which solution the crossing removes. |
There was a problem hiding this comment.
Nothing pins the results that motivate this PR. This testset uses the stock Solovev deck (singfac_min = 1e-4), where the old and new triggers agree. Since it already rewrites the deck per run, a second loop at singfac_min = 1e-8 checking et₁ across the same ucrit ladder is a few lines, and a harness case for one high-n perturbed-equilibrium quantity (an island width at n = 6) would pin the storage fix end to end.
logan-nc
left a comment
logan-nc
left a comment
There was a problem hiding this comment.
One item not in the six open threads, on the kinetic crossing. compute_solution_norms! (head 2baa20b, EulerLagrange.jl:1007) returns early when odet.new && !at_ideal_crossing, resetting the baseline. The kinetic-crossing call at EulerLagrange.jl:847 passes sing_flag = true with no resonant_rows, so when the step before the surface already reduced (odet.new == true), that crossing resets and skips the reduction the flag is supposed to force. Fortran GPEC closed the same hole in PrincetonUniversity/GPEC#332 (fix_required = sing_flag ahead of the branch). The early return is pre-existing on develop, but this PR owns the trigger now, so if odet.new && !at_ideal_crossing && !sing_flag (or passing the surface's rows from the kinetic path too) belongs here.
…nditioning trigger
|
@jhalpern30 what is the status on this? Can this make it in for 2.0.0? |
|
@logan-nc thoughts? This seems like it would be valuable. If Jake doesn't have time before APS then I guess this could be 2.1? Or someone else could get into the weeds on this? |
Probably not. Daniel raised some valid points that I don't have the time/effort to look into before APS |
|
@jhalpern30 @matt-pharr it's time to face reality on v2.0 scope and it seems reasonable that such a fundamental change could come in a minor version update with clean validation. That lets us release the "decades-old approach" for warm familiarity and creates a clean historical record of the improvement. I am not afraid of releasing 2.1 in Dec. so please do keep this on your todo list for immediately post APS. |
|
@jhalpern30 The kinetic-crossing hole I flagged above (Fortran GPEC PR 332) is fixed separately in #626, so it doesn't wait on this PR's trigger change. In This PR rewrites the trigger in the same function, so whichever of the two lands second needs a small rebase there. The changes compose: #626 only replaces the early return, and your conditioning test slots in where the |
Closes #490 and enables #481. Stacked on #489 (base
bugfix/ffs-crossing-resonant-elimination); retarget todeveloponce #489 merges.Accurate ideal energies close to rational surfaces, and trustworthy perturbed equilibria at high n.
singfac_minsingfac_min1e-8. Before, anything below ~1e-5 was unusable on this equilibrium; now accurate down to 1e-8. Prerequisite for #481, whose crossing defect scales withsingfac_min.singfac_min1e-6renorm_thresholdnow insensitive from 10 to 1e8The current full harness finds small DIII-D-like headline energy changes (n = 1 et₁: 4.741e-9 absolute); both Riccati cases are unchanged. Kinetic-field norm changes require further interpretation, as detailed below. Errors are relative to the same build reducing at every step (
renorm_threshold= 1).Review package with the figures, the #490 rerun, the threshold ladder (called ucrit in the original study) and the decision table: https://claude.ai/artifact/JqVHgVrx3xHbZLRJmqNAPe
Release note
solovev_kinetic_calculated, raw displacement/derivative matrix norms change by approximately 5% and the regularized magnetic-field norm by 1.64%; these require a common-grid comparison before physical interpretation. Both Riccati cases are unchanged. (harness @ b8e1e71)renorm_thresholdkeeps its default of 1e3, but in the forward path it now means the U₁ condition-number limit (Riccati keeps the entry-size meaning).You can now move the ideal jump much closer to rational surfaces (
singfac_mindown to 1e-8) and get accurate energies, and high-n perturbed-equilibrium runs no longer silently corrupt their output. The forward solve now applies a Gaussian reduction when the U₁ solution columns lose independence, not when their norms spread apart; the plasma response matrix W = U₂U₁⁻¹ is also formed stably.Regression report
Full SDCC fanout run on 2026-10-28, after merging #489 into this branch:
Baseline: #489 at
a319ed295. Tested head:b8e1e7103. This stacked comparison isolates #520's changes. All 21 cases ran on SDCC sirius using Julia 1.12.6 and CI manifest pin4c9d54fa. Cases used 2 Julia threads / 1 BLAS thread, exceptdiiid_error_fieldat 16 / 8, with matching settings between refs. All 42 tasks, prep, and report job1298807completed with exit code zero.diiid_error_field,diiid_multi_n,diiid_n1,diiid_truncation_n1,solovev_kinetic_calculated,solovev_kinetic_multiion,solovev_kinetic_ntv,solovev_kinetic_nuzero, andsolovev_n1.gal_resistive_peextracts none of its eight quantities on either ref, as before; it provides no numerical validation.solovev_kinetic_calculated,||xi_clebsch_alpha||changes by 5.02%,||dxi_clebsch_psi/dpsi||by 5.01%, and||Jb_theta_reg||by 1.64%, while Re(et₁) changes by 1.935e-7 absolute. These are raw matrix norms; the saved ODE grid also changes from 608 to 595 points. A common-grid comparison is needed to distinguish sampling effects from changes in the physical fields. The report alone does not establish the cause.solovev_kinetic_nuzeroRe(et₁) changes by 7.647e-7 absolute, with saved ODE steps changing from 608 to 597.The earlier parameter-scan results above are historical evidence; this run covers the shipped harness cases and does not rerun those scans or the targeted unit/full-run tests.
SDCC output:
/home/ITER/pharrm/Projects/GPEC_dev/regress-fanout/runs/pr520_b8e1e7103/report.txt.Full regression report (21 cases)
Notes for reviewers
compute_solution_norms!reduces whencolumn_condition(U₁)(1-norm estimate of the column-normalized U₁, from QR viacond(UpperTriangular(R), 1)) exceedsrenorm_threshold. Singular surfaces still always reduce.grow_reduction_storage!). The old 100-slot cap overflowed silently at n ≥ 6 and corrupted PE output (island widths off by 1e73).free_run,free_compute_totalandbuild_result:U₂ / U₁factors U₁ᵀ, so partial pivoting runs across U₁'s column norms, which span ~1e110 just after a reduction. When the last event before the edge is a reduction, W was garbage (DIIID n = 4, δm = 24: et₁ −0.026 instead of 0.157). Older than this PR; the new trigger makes it more likely.renorm_threshold1e3 rather than 1e6 (1e6 is equally accurate in the forward path; 1e3 keeps Riccati, which shares the key, on its tested value); redefiningrenorm_thresholdrather than adding a key; cond(U₁) rather than cond[U₁;U₂]; keeping elimination rather than an (S, I) or QR reset; no dedicated test for thelufix (no small matrix reproduces it).