Repository navigation
ForceFreeStates.Galerkin - BUGFIX! - Jacobi-scale the banded Galerkin LU solve - #512
Merged
Merged
Conversation
… LU solve The assembled Galerkin matrix diagonal spans ~2e-15..2e10, and unscaled banded LU (gbtrf/gbtrs) lost the small resonant coefficients that carry Delta': on DIII-D g147131 n=1 the 4/1 row varied erratically with gal_nx and disagreed with Riccati. Symmetrically scale A <- D A D, d_j = |A_jj|^(-1/2) (1 for a zero diagonal), before factoring, and un-scale every solution column (incl. rpec coil columns). Add test/runtests_galerkin.jl: a badly scaled banded system that unscaled LU solves to ~1e-4 relative error, now ~1e-15. Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
github-actions
Bot
added
bugfix
Contributor
|
This pull request is missing an assignee. If you are not ready to name them, mark this pull request as a draft. |
matt-pharr
approved these changes
Oct 8, 2026
matt-pharr
left a comment
matt-pharr
left a comment
Collaborator
There was a problem hiding this comment.
@StuartBenjamin this is YUGE and I wish I had the numerical linear algebra knowledge to have thought of this myself. I am very excited to do a broader sensitivity study to see how much this improves RDCON's convergence/accuracy. This is a nice atomic addition, thank you.
matt-pharr
enabled auto-merge
October 8, 2026 13:15
matt-pharr
disabled auto-merge
October 8, 2026 13:37
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
2 participants
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
integrator = "galerkin") results move. Ingal_resistive_diiid: Δ′ diagonal 0.05%, ‖Δ′‖ 0.01%, ‖Δ_coil‖ 1.85%, ‖match cout‖ 0.40%. Forward, Riccati and ideal results are bit-identical (diiid_n1,diiid_n1_riccati). (harness @ 8a8a3cf)integrator = "galerkin"cases and discard previously computed Galerkin Δ′, which could be wrong by O(1), including sign, on outer surfaces.The Galerkin outer-region Δ′ solve lost accuracy in the small resonant coefficients that carry Δ′. As a result, some Δ′ entries changed erratically with
gal_nxand could even flip sign, for example the 4/1 row of DIII-D 147131 n=1. The banded LU system is now Jacobi-scaled before factorization. Δ′ now converges ingal_nx, and the 2/1–4/1 diagonal agrees with Riccati to within 0.4%.Regression report
Notes for reviewers
Harness:
gal_resistive_peis N/A on both refs. This is pre-existing ondevelopand not caused by this PR. With the Galerkin integrator, PerturbedEquilibrium skips the plasma-response and singular-coupling stages ("free_boundarywas not produced by the galerkin integrator"). None of that case's quantities are written on either commit. The Galerkin Δ′ it shares withgal_resistive_diiidis covered there.Cause. The assembled Galerkin matrix diagonal spans about 25 decades: resonant-cell entries are around 1e-14, sitting next to extension-cell entries around 1e-3. Unscaled partial-pivoting banded LU (
gbtrf!/gbtrs!) loses the small resonant coefficients, and Δ′ is extracted from exactly those.Fix (
src/ForceFreeStates/Galerkin/GalerkinSolve.jl, LU branch only):gal_scaled_lu_solve!applies symmetric Jacobi scaling,A ← D·A·Dwithd_j = |A_jj|^(-1/2)(1 for a zero diagonal), ingbtrf!band storage.gal_solver = "cholesky"branch is untouched. It still fails ondevelopwithUndefVarError: pbtrf!. A separate feature PR (feature/galerkin-cholesky) fixes that. The two PRs edit disjoint lines and merge cleanly in either order.Evidence: DIII-D g147131 n=1, no wall, Re Δ′ diagonal,
gal_solver = "LU"galerkin_pr_package.zip
All runs on
develop@6553fe205("before") and this PR @8a8a3cf17("after"), same decks. Riccati is identical on both commits, which serves as a control.gal_nx = 128gal_nx = 256gal_nx = 384gal_nx = 512gal_nx = 128gal_nx = 256gal_nx = 384gal_nx = 512gal_nx = 512was[-29.2, -184.7, 125.8, 2482, 5474]before. After, it is[3.192, 5.214, -60.50, 57.84, 68.17], against Riccati's[3.184, 5.196, -60.48, 57.90, 77.06].gal_nxbut stay 7–10% from Riccati. That is a separate edge/vacuum-dominated discrepancy and is not addressed here.vac_flag = false,gal_nx = 256). The diagonal moves slightly (2.033 / -14.48 / -76.32 / -207.0 / -1968 → 2.041 / -14.33 / -78.38 / -207.5 / -1933). The 4/1 off-diagonals change substantially ([-6.86, -110.1]→[18.4, 0.06]in the 5/1 and 6/1 columns).In-repo example decks.
DIIID-like_gal_resistive{,_pe}_exampleback harness casesgal_resistive_diiidandgal_resistive_pe. They are well enough conditioned that they move only slightly:Tests. The new
test/runtests_galerkin.jlbuilds a banded system with an alternating 4e-8 / 1e5 diagonal scaling. Unscaled LU solves it to about 1e-4 relative error; the scaled solve reaches about 1e-15. The test also covers multiple-RHS un-scaling and a zero diagonal entry. At8a8a3cf17(1 CPU job):Galerkin scaled banded LUsolve API(includes the Galerkinsolvetestset)ForceFreeStatesResultRunner: Control + run_slayer + HDF5 outputgpec.h5 schema namingReproducing. All commands run from a GPEC checkout.
Unit test (on this branch):
Regression report:
julia --project=regression-harness regression-harness/regress.jl \ --cases diiid_n1,gal_resistive_diiid,gal_resistive_pe,diiid_n1_riccati \ --refs origin/develop,bugfix/galerkin-scaled-lu-pr --forceBefore/after Δ′ tables. The attached
galerkin_pr_package.zip
contains
reproduce/:validate.jldriver;run_validation.sh, run once against a checkout at6553fe205and once at8a8a3cf17;summarize.py, which prints the tables above;The expected output and the logs of the runs cited here are in
validation_8a8a3cf/.Δ′ test coverage found while preparing this (not addressed here):
gal_resistive_diiidandgal_resistive_peare the only value-level Galerkin Δ′ trackers.diiid_n1_riccatitracks Riccati Δ′.docs/src/developer_notes.mdsays the Δ′ diagonal is pinned intest/runtests_parallel_integration.jl. That file now only has a comment saying the per-surface assertions were removed.diiid_n1sample report in the docs shows adelta primerow that the currentdiiid_n1.tomlno longer tracks; it tracks only the ballooning Δ′ checksum.