Repository navigation
Tearing - FEATURE - Classify tearing roots by winding number - #597
matt-pharr merged 5 commits into
Conversation
Add root_classification = "winding" to [SLAYER] and classification = :winding to find_growth_rates. Each contour intersection is classified by the winding number of the dispersion residual on a small circle around it (argument principle): a root, a pole, or a near-miss, instead of the pole threshold and the above-pole/geometry filters. The SLAYER residual is a function of conj(Q), so its zeros wind -1; zero_winding = +1 serves holomorphic residuals. The reported root is the least stable zero, polished inside its circle; winding numbers are stored on GrowthRateResult and written to Tearing/Diagnostics/*Winding. The default stays "threshold". Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
… analysis and harden its circles Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
…e coupled root in develop's normalization Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
… root classification Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
fbdd6eb to
18d88c9
Compare
Describe finite phase sampling as an estimate rather than a proof. Note midpoint and count-convergence checks as possible defensive extensions without changing classification or polishing.
|
Defensive winding-count hardening is deferred; the current classifier remains unchanged. A synthetic counterexample shows that sampled phase increments below pi/2 do not certify the winding count: on the unit circle, f(z) = (z - 0.999 exp(ipi/64)) / (z - 1.001 exp(ipi/64)) has winding +1, but the initial 64 samples alias it to zero. A possible future extension is to require midpoint variation and agreement across successive refinements, returning unresolved when the evaluation budget is exhausted. We tested that prototype on SDCC against this PR's rebased head, 18d88c9, using the DIII-D-like uncoupled case, alternate chi_perp_e transport, coupled case, and a coarse coupled scan. Roots, frequencies, growth rates, poles, and flags were unchanged in all four original/proposed comparisons. Both methods agreed with all 41 independently audited production circles. In 32 radius-stress checks, every finite accepted count agreed with the dense reference; the prototype additionally returned unresolved for one circle where the original and dense reference correctly returned zero. References used phase-offset dense grids and refinement checks; these are numerical evidence, not mathematical certification. The prototype required 4–5 times as many classifier evaluations on production circles, adding roughly 2–4 seconds in the measured classifier work. All 369 focused checks passed on the prototype. This establishes a synthetic failure mode and a possible defensive improvement, but does not demonstrate a classification failure in the tested production cases. The hardening is therefore not being included as a numerical change in this PR. A documentation-only correction removes the claims that the sampled count is an exact proof and records hardening as a possible extension. It is committed locally as 36d54d2 and is not pushed yet. Executable syntax was verified unchanged. Full regression-harness provenance remains the existing fbdd6eb run; the SDCC study is supplemental validation of 18d88c9 plus the separately recorded prototype, not a new full harness stamp. SDCC evidence: array 1300248, report job 1300259, /home/ITER/pharrm/Projects/GPEC_dev/pr597-winding-validation-20261029/report.txt. Prototype patch SHA-256: 1b0b2825f3a77f3c4f1bca5877568a58da9745eddabea3b9fbfcf2e574771d5a. Human review/approval remains with Matthew; no automated approval or merge was performed. |
matt-pharr
left a comment
There was a problem hiding this comment.
This is good besides some wordy documentation. Went through and tested it pretty extensively with claude, I think this is fine as the version for 2.0.0. It came up with special cases where it might not work, but we couldn't produce those so it might be fake lol. Approve
Auto-merge was disabled for another pull request in the stack
Release note
root_classification = "winding"is opt-in (harness @ 36d54d2)Tearing roots can now be classified by the argument principle. With
root_classification = "winding"in[SLAYER](orclassification = :windinginfind_growth_rates), each contour intersection is a root, a pole or a near-miss according to the winding number of the dispersion residual on a small circle around it, instead of the pole threshold and the above-pole/geometry filters. The reported root is the least stable zero, and the winding numbers are written toTearing/Diagnostics/*Winding.Changes since the review
:geom/:gapflags. It is now a pluggable classification step inside_run_analysis(ea316dcfe), so later changes to selection or flags apply to both modes.ea316dcfe,3d8555ba4): the radius follows the local contour resolution instead of one global length; a non-finite circle is retried smaller; the wide recheck may reach 0.9 of the nearest-intersection distance. The never-firing non-integer check is gone, and the docstring describes the argument-principle assumptions; sampled phase increments alone do not certify a count.zero_winding = −1is defined forSLAYERModel{:fitzpatrick}only, so a future analytic variant does not inherit it silently.diiid_slayer_n1_windingruns the SLAYER deck withroot_classification = "winding".abbd767ef); the conflicts were the additions of Tearing - FEATURE - Add critical resonant field calculations to the tearing workflow #415 and InnerLayer.SLAYER - FEATURE! - Build P_perp from the perpendicular particle diffusivity by default and add an ion heat diffusivity profile #501 to the same Runner lines.Regression report
Latest full SDCC comparison (2026-10-29):
aac35ab68(parent) versus36d54d2fa(this PR plus the documentation-only correction), all 22 cases / 44 tasks on sirius, Julia 1.12.6, pinned manifest4c9d54fa, 8 Julia threads / 4 BLAS threads. All 296 comparable quantities across 20 shared cases are unchanged, including the uncoupled, alternate-transport, and coupled SLAYER cases. The newdiiid_slayer_n1_windingcase succeeds on36d54d2fa(203.5 s) and fails on the parent, which does not supportroot_classification; this is not a successful two-ref comparison for that new case.gal_resistive_pestill has eight missing/N/A quantities on both refs and provides no numerical validation. No other task failed. Report job1300276; SDCC report:/home/ITER/pharrm/Projects/GPEC_dev/regress-fanout/runs/pr597_36d54d2fa_full/report.txt.The tested documentation-only correction
36d54d2fais now the PR head. Numerical behavior is unchanged by that commit.Earlier comparison:
regress --refs 5828895a0,fbdd6eb07 --force(all cases, against #577's head) on feynman (SLURM, 8 threads, Julia 1.11.9, manifest pinned): every case is unchanged, includingdiiid_slayer_n1(17) anddiiid_slayer_n1_coupled(7), so threshold mode is untouched. The newdiiid_slayer_n1_windingcase fails on the base, which has noroot_classificationkey, so the harness exits 1;gal_resistive_peis N/A on both refs. Atfbdd6eb07the winding case finds the same three roots as threshold mode, with γ (291.1, −339.4, −1395.3 s⁻¹) agreeing to within 6e-8 relative and ω to within 2e-8, in 277 s.runtests_dispersion_winding.jl(56),runtests_slayer_runner.jl(156),runtests_dispersion_scan.jl(33),runtests_dispersion_amr.jl(50) andruntests_dispersion_polish.jl(65) pass.Supplemental SDCC validation after rebasing
At rebased head
18d88c91c, four DIII-D-like original/proposed comparisons (uncoupled, alternate transport, coupled, coarse coupled) produced identical roots, frequencies, growth rates, poles and flags. Both classifiers agreed with all 41 audited production circles. Of 32 radius-stress checks, the proposed hardening returned unresolved once where the original and dense reference agreed on zero; all finite accepted counts agreed with the reference. The prototype passed 369 focused checks but cost 4–5 times the production-circle evaluations. It is deferred as a possible future extension, not included in this PR's numerical behavior.The documentation-only correction
36d54d2fahas executable syntax identical to18d88c91cand was subsequently tested by the full SDCC harness reported above. The hardening prototype remains deferred.Notes for reviewers
Why. The threshold classifier decides root vs pole from the size of Re Δ along the Im Δ = 0 contour, then prunes with the above-pole and
:geom/:gapheuristics. None of these is a proof, and the threshold is scale-dependent. In coupled SLAYER scans of a current ramp-up, the threshold path listed two intersections as poles and reported no root; both wind −1 on circles of three radii, i.e. they are isolated zeros (that deck used the absolutepole_threshold = 10against |det| ~ 2e8). Near a root–pole crossing the labels are equally unreliable. The argument principle counts zeros minus poles inside the circle; the implementation estimates this from finite samples and can miss rapid phase turns between samples. Classification does not depend on an amplitude threshold. We have used this acceptance rule out of tree (a driver aroundrun_slayer_from_inputs) on every point of the coupling-paper scans since 2026-10-07; this PR moves it into the package so other analyses, and GGJ (stacked PR), can share it.Orientation. The SLAYER residual is a function of conj(Q) (
solve_innerevaluates conj(Δ̂(ĝ)), the same fact #577 fixes the polish for), so its zeros wind −1 and its poles +1. A holomorphic residual, e.g. the GGJ full-matching determinant, has the opposite sign.find_growth_ratestakeszero_winding = ±1; the runner sets −1 forSLAYERModeland refuses other models in winding mode for now. The tests show that the wrong sign swaps roots and poles.Method. The winding classifier (
_winding_classifyinGrowthRateExtraction.jl) is one step inside the shared_run_analysis, in place of the pole threshold; everything after it (polish, the +γ-step flag, root selection,:geom/:gapwarnings, the secondary root, flags) is the same code in both modes.winding_rfrac(default 0.5) × the local contour-segment length at the intersection, so circles follow the cell size on AMR meshes, capped at 0.45 × the distance to the nearest other intersection;:winding_unresolved;filtered_rootswith flag:spurious;The coupled path passes the determinant as
residualwithpolish = false, so in winding mode the determinant is evaluated around intersections only to count; in threshold mode it is not passed and nothing changes.Interface.
GrowthRateResultgainsclassificationand three winding vectors aligned withvalid_roots/poles/filtered_roots(empty in threshold mode); the 13-argument constructor still works.find_growth_ratesgainsclassification,zero_winding,winding_rfracandpolish(defaulttrue, soresidualalone keeps polishing as before). HDF5:Diagnostics/{ValidRoots,Poles,FilteredRoots}Windingare written only when a winding classification ran, so threshold-mode files are unchanged._write_ragged_complex!is now_write_ragged!and takes Float64 rows as well.Cost. Each intersection costs 64 to a few thousand residual evaluations. On
diiid_slayer_n1the winding variant runs in 271 s against 279 s for threshold mode, and its growth rates agree with threshold mode to about 1e-8 relative.Not in this PR. Zeros the scan does not resolve as contour crossings are still missed, in both modes. The default stays
threshold; switching it would move results and belongs in its own PR.Stacked on #577 (
bugfix/slayer-polish-antianalytic).🤖 Generated with Claude Code