From 8b5dd913eed956d9d1ce0eacbddf27691da0dc94 Mon Sep 17 00:00:00 2001 From: d-burg Date: Mon, 26 Oct 2026 16:40:09 -0400 Subject: [PATCH 1/5] Tearing - FEATURE - Classify tearing roots by winding number 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 --- docs/development/hdf5-conventions.md | 2 +- .../Dispersion/GrowthRateExtraction.jl | 301 ++++++++++++++++-- src/Tearing/Runner/Control.jl | 15 +- src/Tearing/Runner/HDF5Output.jl | 33 +- src/Tearing/Runner/run_slayer.jl | 27 +- test/runtests.jl | 1 + test/runtests_dispersion_winding.jl | 119 +++++++ test/runtests_slayer_runner.jl | 62 ++++ 8 files changed, 511 insertions(+), 49 deletions(-) create mode 100644 test/runtests_dispersion_winding.jl diff --git a/docs/development/hdf5-conventions.md b/docs/development/hdf5-conventions.md index 40f381ebc..e106a9f30 100644 --- a/docs/development/hdf5-conventions.md +++ b/docs/development/hdf5-conventions.md @@ -46,7 +46,7 @@ Top level (11 groups): | `PerturbedEquilibrium/` | `ForcingModes/`, `Response/`, `ResponseMatrices/`, `SingularCoupling/`, `Energies/`, `TorqueResponse/` (kinetic forward solves: ψ-resolved `T_xe`, `T_source`), control-surface spectra | | `KineticForces/` | `/` (torque/energy profiles, `dTdpsi_ell` resolved by bounce harmonic on the `ell` scale, `KineticMatrices/`); multi-ion runs add `PerSpecies///` with the same per-method layout, summing to the top-level total | | `ErrorFields/` | `CoilSensitivities/` (per-coil-set control-surface spectra and their rigid shift/tilt derivatives, `DominantMode/` full-window projection); `MonteCarlo/` (intrinsic and corrected `\|δ\|` histograms over the sampled tolerances, per batch and averaged); `Risk/` (threshold density, `P(lock\|δ)`, locking probabilities, `ToleranceScan/`); `NTV/` (correction-coil overlap and NTV torque couplings per kAt) | -| `Tearing/` | `PerSurface/` (+ `DpMatrix/`), `Roots/`, `LayerWidths/`, `Diagnostics/{ValidRoots,Poles,FilteredRoots}`, `Scan/Surface_/`, `CriticalResonantField/` (+ `Scan/Surface_/`) | +| `Tearing/` | `PerSurface/` (+ `DpMatrix/`), `Roots/`, `LayerWidths/`, `Diagnostics/{ValidRoots,Poles,FilteredRoots}` (+ `*Winding` with `root_classification = winding`), `Scan/Surface_/`, `CriticalResonantField/` (+ `Scan/Surface_/`) | | `SurfaceGeometries/` | `{Plasma,Wall}/{x,y,z}` point clouds | Reserved (documented, not yet written): `ForceFreeStates/Solutions/RiccatiIntegration/` — the third integrator backend slot alongside `ForwardIntegration` and `GalerkinIntegration`. diff --git a/src/Tearing/Dispersion/GrowthRateExtraction.jl b/src/Tearing/Dispersion/GrowthRateExtraction.jl index d16667889..db8582972 100644 --- a/src/Tearing/Dispersion/GrowthRateExtraction.jl +++ b/src/Tearing/Dispersion/GrowthRateExtraction.jl @@ -19,6 +19,11 @@ # Roots whose +γ step exits the loop AND that lie above the highest # pole are filtered out (spurious upper branches). # 5. Return the highest-γ surviving root in physical units. +# +# With `classification = :winding`, steps 3-4 are replaced by the argument +# principle: each intersection is a root, a pole or neither according to the +# winding number of the residual on a small circle around it (see +# `_run_analysis_winding`). using Contour using DelaunayTriangulation @@ -33,30 +38,36 @@ using Random: Xoshiro Output of `find_growth_rates`. -| field | meaning | -|:-------------------- |:----------------------------------------------------- | -| `Q_root` | Best (highest-γ surviving) root, normalized | -| `omega_Hz` | `Re(Q_root) / tauk` — physical rotation frequency | -| `gamma_Hz` | `Im(Q_root) / tauk` — physical growth rate | -| `Q_root_secondary` | Second-most-unstable root flagged for ambiguity, or | -| | `NaN+NaNim` if the primary root was unambiguous. | -| `omega_Hz_secondary` | physical ω of the secondary root, or 0 if none | -| `gamma_Hz_secondary` | physical γ of the secondary root, or 0 if none | -| `warning_flags` | `Vector{Symbol}` of warnings raised on `Q_root`: | -| | `:geom`, `:gap` (root accepted with caveat); | -| | `:spurious` when ≥1 contour near-miss was dropped by | -| | the validity gate (parked in `filtered_roots`); or | -| | `:no_root` when NO usable root was found (`Q_root` is | -| | `NaN`; `omega_Hz`/`gamma_Hz` fall back to 0 — check | -| | this flag to tell that apart from a true γ≈0 result). | -| | Empty if root is clean. | -| `valid_roots` | All non-pole intersections that survived pole filter | -| `poles` | Intersections classified as poles | -| `filtered_roots` | Intersections rejected by the above-pole/outside-Re | -| | filter or the new geom+gap recursion | -| `re_contours` | Extracted Re(Δ)=`re_target` polylines | -| `im_contours` | Extracted Im(Δ)=`im_target` polylines | -| `pole_threshold` | Threshold used for pole classification | +| field | meaning | +|:------------------------ |:----------------------------------------------------- | +| `Q_root` | Best (highest-γ surviving) root, normalized | +| `omega_Hz` | `Re(Q_root) / tauk` — physical rotation frequency | +| `gamma_Hz` | `Im(Q_root) / tauk` — physical growth rate | +| `Q_root_secondary` | Second-most-unstable root flagged for ambiguity, or | +| | `NaN+NaNim` if the primary root was unambiguous. | +| `omega_Hz_secondary` | physical ω of the secondary root, or 0 if none | +| `gamma_Hz_secondary` | physical γ of the secondary root, or 0 if none | +| `warning_flags` | `Vector{Symbol}` of warnings raised on `Q_root`: | +| | `:geom`, `:gap` (root accepted with caveat); | +| | `:spurious` when ≥1 contour near-miss was dropped by | +| | the validity gate (parked in `filtered_roots`); or | +| | `:no_root` when NO usable root was found (`Q_root` is | +| | `NaN`; `omega_Hz`/`gamma_Hz` fall back to 0 — check | +| | this flag to tell that apart from a true γ≈0 result). | +| | Empty if root is clean. | +| `valid_roots` | All non-pole intersections that survived pole filter | +| `poles` | Intersections classified as poles | +| `filtered_roots` | Intersections rejected by the above-pole/outside-Re | +| | filter or the new geom+gap recursion | +| `re_contours` | Extracted Re(Δ)=`re_target` polylines | +| `im_contours` | Extracted Im(Δ)=`im_target` polylines | +| `pole_threshold` | Threshold used for pole classification | +| | (`:threshold` mode only) | +| `classification` | `:threshold` or `:winding` — how intersections were | +| | classified (see `find_growth_rates`) | +| `valid_roots_winding` | winding number of each entry of `valid_roots`, | +| `poles_winding` | `poles` and `filtered_roots` (`NaN` where the circle | +| `filtered_roots_winding` | was unresolved); empty in `:threshold` mode | """ struct GrowthRateResult Q_root::ComplexF64 @@ -72,8 +83,22 @@ struct GrowthRateResult re_contours::Vector{Vector{ComplexF64}} im_contours::Vector{Vector{ComplexF64}} pole_threshold::Float64 + classification::Symbol + valid_roots_winding::Vector{Float64} + poles_winding::Vector{Float64} + filtered_roots_winding::Vector{Float64} end +# Threshold-mode constructor (no winding numbers), the form used before +# winding classification existed. +GrowthRateResult(Q_root, omega_Hz, gamma_Hz, Q_root_secondary, omega_Hz_secondary, + gamma_Hz_secondary, warning_flags, valid_roots, poles, filtered_roots, + re_contours, im_contours, pole_threshold) = + GrowthRateResult(Q_root, omega_Hz, gamma_Hz, Q_root_secondary, omega_Hz_secondary, + gamma_Hz_secondary, warning_flags, valid_roots, poles, filtered_roots, + re_contours, im_contours, pole_threshold, :threshold, + Float64[], Float64[], Float64[]) + """ find_growth_rates(scan::ScanResult, tauk::Real; re_target=0.0, im_target=0.0, @@ -122,6 +147,40 @@ single-surface scans; `mc.surfaces[mc.ref_idx].tauk` for coupled scans). Δ' BVP failed → `|Δ'|` huge and uncancellable) are dropped from `valid_roots`, parked in `filtered_roots`, and flagged `:spurious`. If every candidate is spurious, the result carries `:no_root`. + - `polish` -- polish roots when a `residual` is supplied (default `true`). + `false` keeps the raw contour estimate, e.g. to classify by winding number + without polishing. + - `classification` -- `:threshold` (default) or `:winding`; see below. + - `zero_winding` -- `:winding` mode only: the winding number of an isolated + zero of `residual`, `+1` for a function holomorphic in Q and `-1` for one + that is anti-holomorphic (a function of conj(Q)), such as the SLAYER + residual. Default `-1`. + - `winding_rfrac` -- `:winding` mode only: radius of the classifying circle + in units of the median contour-segment length (default 0.5), capped at + 0.45 × the distance to the nearest other intersection. + +# Winding classification (`classification = :winding`) + +Each distinct contour intersection is classified by the argument principle +instead of the pole-threshold heuristics. The counter-clockwise winding number +of `residual` on a small circle around the intersection counts the zeros minus +the poles inside it (for a holomorphic function; the sign is reversed for an +anti-holomorphic one, hence `zero_winding`). The circle is resampled until no +phase step between neighboring samples exceeds π/2, so the count is exact +rather than estimated. An intersection is a root when the count has the sign of +`zero_winding`, a pole when it has the opposite sign, and is moved to +`filtered_roots` (flag `:spurious`) when the count is zero, e.g. a contour +near-miss or an unresolved pole-zero pair. A zero count is checked once more on +a circle eight times wider (still clear of every other intersection), because +the contour estimate of a pole can lie outside the small circle. A circle that +stays unresolved (winding `NaN`, a zero or pole on it) is filtered too and flags +`:winding_unresolved`. `pole_threshold`, `filter_above_poles`, +`filter_outside_re`, `gap_kHz_threshold` and `validity_rtol` decide nothing in +this mode, and `Q_root` is the least stable root, with no secondary root. +With `polish = true` a root is polished inside its classifying circle, which is +proven to contain it. `residual` is required. Each intersection costs between +64 and a few hundred residual evaluations. Only zeros that the scan resolves +as contour crossings are found, as in `:threshold` mode. # Spurious-root recursion @@ -148,7 +207,12 @@ function find_growth_rates(scan::ScanResult, tauk::Real; gap_kHz_threshold::Real=1.0, residual=nothing, polish_maxit::Integer=20, - validity_rtol::Real=1e-3) + validity_rtol::Real=1e-3, + polish::Bool=true, + classification::Symbol=:threshold, + zero_winding::Integer=-1, + winding_rfrac::Real=0.5) + _check_classification(classification, residual, zero_winding) return _extract_growth_rates(scan.re_axis, scan.im_axis, scan.Δ, Float64(tauk); re_target=Float64(re_target), @@ -160,7 +224,11 @@ function find_growth_rates(scan::ScanResult, tauk::Real; residual=residual, polish_maxit=Int(polish_maxit), validity_scale=_residual_scale(scan.Δ), - validity_rtol=Float64(validity_rtol)) + validity_rtol=Float64(validity_rtol), + polish=polish, + classification=classification, + zero_winding=Int(zero_winding), + winding_rfrac=Float64(winding_rfrac)) end """ @@ -186,7 +254,12 @@ function find_growth_rates(amr::AMRResult, tauk::Real; gap_kHz_threshold::Real=1.0, residual=nothing, polish_maxit::Integer=20, - validity_rtol::Real=1e-3) + validity_rtol::Real=1e-3, + polish::Bool=true, + classification::Symbol=:threshold, + zero_winding::Integer=-1, + winding_rfrac::Real=0.5) + _check_classification(classification, residual, zero_winding) return _extract_growth_rates_amr(amr.Q, amr.Δ, Float64(tauk); re_target=Float64(re_target), im_target=Float64(im_target), @@ -197,7 +270,11 @@ function find_growth_rates(amr::AMRResult, tauk::Real; residual=residual, polish_maxit=Int(polish_maxit), validity_scale=_residual_scale(amr.Δ), - validity_rtol=Float64(validity_rtol)) + validity_rtol=Float64(validity_rtol), + polish=polish, + classification=classification, + zero_winding=Int(zero_winding), + winding_rfrac=Float64(winding_rfrac)) end # --------------------------------------------------------------------- @@ -584,7 +661,15 @@ function _run_analysis(re_paths::Vector{Vector{ComplexF64}}, residual=nothing, polish_maxit::Int=20, validity_scale::Float64=0.0, - validity_rtol::Float64=1e-3) + validity_rtol::Float64=1e-3, + polish::Bool=true, + classification::Symbol=:threshold, + zero_winding::Int=-1, + winding_rfrac::Float64=0.5) + classification === :winding && + return _run_analysis_winding(re_paths, im_paths, tauk, residual; + zero_winding=zero_winding, rfrac=winding_rfrac, + polish=polish, polish_maxit=polish_maxit, pole_threshold=pole_threshold) raw_intersections = _all_intersections(re_paths, im_paths) poles = ComplexF64[] @@ -595,7 +680,7 @@ function _run_analysis(re_paths::Vector{Vector{ComplexF64}}, # (non-pole) contour intersection to the true zero of `f`. Trust radii are # capped by the spacing to neighbouring intersections so closely-spaced # coupled roots stay coherently bound to their own contour crossing. - do_polish = residual !== nothing + do_polish = polish && residual !== nothing h_grid = do_polish ? _median_segment_length(re_paths, im_paths) : 0.0 # Validity gate: a genuine root drives |residual| ≈ 0; a spurious contour # near-miss (e.g. a surface whose Δ' BVP failed → |Δ'| huge and @@ -768,6 +853,138 @@ function _run_analysis(re_paths::Vector{Vector{ComplexF64}}, re_paths, im_paths, pole_threshold) end +# --------------------------------------------------------------------- +# Winding-number classification (`classification = :winding`). +# +# The pole-threshold path above decides root vs pole from the local size of +# Re(Δ) along the Im=0 contour, then prunes with geometric heuristics. Both +# can misjudge a root that sits close to a pole, and none of them is a proof. +# The argument principle is: the winding number of f on a closed curve counts +# zeros minus poles inside it (holomorphic f; the sign flips for a function of +# conj(Q)). Evaluated on a small circle around each intersection with the +# phase resolved, it classifies the intersection exactly. +# --------------------------------------------------------------------- + +const _VALID_CLASSIFICATIONS = (:threshold, :winding) + +function _check_classification(classification::Symbol, residual, zero_winding::Integer) + classification in _VALID_CLASSIFICATIONS || + throw(ArgumentError("find_growth_rates: classification=$(classification) " * + "not in $(_VALID_CLASSIFICATIONS)")) + classification === :winding || return nothing + residual === nothing && + throw(ArgumentError("find_growth_rates: classification=:winding needs a residual callable")) + abs(zero_winding) == 1 || + throw(ArgumentError("find_growth_rates: zero_winding must be +1 (holomorphic residual) " * + "or -1 (anti-holomorphic), got $(zero_winding)")) + return nothing +end + +# Counter-clockwise winding number of `f` on the circle |Q − Q0| = r. The +# circle is sampled at `nmin` points and the sampling doubled (reusing the +# previous samples) until no phase step between neighbors exceeds π/2, up to +# `nmax` points. Returns (W, max_step); W is NaN when `f` is non-finite or zero +# on the circle or the phase stays unresolved (a zero or pole on the circle). +function _winding_number(f, Q0::ComplexF64, r::Float64; nmin::Int=64, nmax::Int=4096) + N = nmin + vals = ComplexF64[f(Q0 + r * cis(2π * j / N)) for j in 0:N-1] + max_step = Inf + while true + all(z -> isfinite(z) && !iszero(z), vals) || return (NaN, Inf) + total = 0.0 + max_step = 0.0 + @inbounds for j in 1:N + d = angle(vals[mod1(j + 1, N)] / vals[j]) + total += d + max_step = max(max_step, abs(d)) + end + max_step < π / 2 && return (total / 2π, max_step) + 2N > nmax && return (NaN, max_step) + # double the sampling: old samples at even positions, new ones between + new = Vector{ComplexF64}(undef, 2N) + @inbounds for j in 1:N + new[2j-1] = vals[j] + new[2j] = f(Q0 + r * cis(2π * (2j - 1) / (2N))) + end + vals = new + N *= 2 + end +end + +# Merge intersections that coincide to within `tol` (contours crossing at a +# shared vertex are reported once per adjacent segment pair). +function _dedupe_points(pts::Vector{ComplexF64}, tol::Float64) + out = ComplexF64[] + for p in pts + any(q -> abs(q - p) <= tol, out) || push!(out, p) + end + return out +end + +function _run_analysis_winding(re_paths::Vector{Vector{ComplexF64}}, + im_paths::Vector{Vector{ComplexF64}}, + tauk::Float64, residual; + zero_winding::Int, rfrac::Float64, polish::Bool, polish_maxit::Int, + pole_threshold::Float64) + h_grid = _median_segment_length(re_paths, im_paths) + pts = _dedupe_points(_all_intersections(re_paths, im_paths), 1e-9 * max(h_grid, eps())) + + valid_roots, poles, filtered = ComplexF64[], ComplexF64[], ComplexF64[] + w_valid, w_poles, w_filtered = Float64[], Float64[], Float64[] + warning_flags = Symbol[] + + for (i, pt) in enumerate(pts) + # circle radius: a fraction of the local resolution, and less than half + # the distance to any other intersection so circles never overlap + r = _polish_trust_radius(pts, i, h_grid; k_cell=rfrac, nn_frac=0.45) + W, _ = _winding_number(residual, pt, r) + # The contour estimate of a zero is accurate (the field is nearly linear + # there); that of a pole is not, and can fall outside the small circle. + # If it counts nothing, count once more on a wider circle, still clear + # of every other intersection. + if !(isfinite(W) && round(Int, W) != 0) + r_wide = _polish_trust_radius(pts, i, h_grid; k_cell=8rfrac, nn_frac=0.45) + if r_wide > r + W_wide, _ = _winding_number(residual, pt, r_wide) + if isfinite(W_wide) && round(Int, W_wide) != 0 + W, r = W_wide, r_wide + end + end + end + n = isfinite(W) ? round(Int, W) : 0 + if !isfinite(W) || abs(W - n) > 0.1 + push!(filtered, pt) + push!(w_filtered, NaN) + :winding_unresolved in warning_flags || push!(warning_flags, :winding_unresolved) + elseif n * zero_winding >= 1 # net zeros inside the circle + if polish + pt, _, _ = _polish_root(residual, pt, r; maxit=polish_maxit) + end + push!(valid_roots, pt) + push!(w_valid, W) + elseif n * zero_winding <= -1 # net poles inside + push!(poles, pt) + push!(w_poles, W) + else # zero count: near-miss or pole-zero pair + push!(filtered, pt) + push!(w_filtered, W) + :spurious in warning_flags || push!(warning_flags, :spurious) + end + end + + Q_root = isempty(valid_roots) ? ComplexF64(NaN, NaN) : valid_roots[argmax(imag.(valid_roots))] + isnan(real(Q_root)) && push!(warning_flags, :no_root) + omega_Hz = isnan(real(Q_root)) ? 0.0 : real(Q_root) / tauk + gamma_Hz = isnan(imag(Q_root)) ? 0.0 : imag(Q_root) / tauk + + return GrowthRateResult(Q_root, omega_Hz, gamma_Hz, + ComplexF64(NaN, NaN), 0.0, 0.0, + warning_flags, + valid_roots, poles, filtered, + re_paths, im_paths, pole_threshold, + :winding, w_valid, w_poles, w_filtered) +end + # Regular-grid path: extract contours via Contour.jl, compute im_re_vals by # bilinear interpolation on the grid, then run the shared analysis. function _extract_growth_rates(re_axis::Vector{Float64}, @@ -783,7 +1000,11 @@ function _extract_growth_rates(re_axis::Vector{Float64}, residual=nothing, polish_maxit::Int=20, validity_scale::Float64=0.0, - validity_rtol::Float64=1e-3) + validity_rtol::Float64=1e-3, + polish::Bool=true, + classification::Symbol=:threshold, + zero_winding::Int=-1, + winding_rfrac::Float64=0.5) re_field = real.(Δ_grid) im_field = imag.(Δ_grid) @@ -803,7 +1024,11 @@ function _extract_growth_rates(re_axis::Vector{Float64}, residual=residual, polish_maxit=polish_maxit, validity_scale=validity_scale, - validity_rtol=validity_rtol) + validity_rtol=validity_rtol, + polish=polish, + classification=classification, + zero_winding=zero_winding, + winding_rfrac=winding_rfrac) end # --------------------------------------------------------------------- @@ -958,7 +1183,11 @@ function _extract_growth_rates_amr(Q::Vector{ComplexF64}, residual=nothing, polish_maxit::Int=20, validity_scale::Float64=0.0, - validity_rtol::Float64=1e-3) + validity_rtol::Float64=1e-3, + polish::Bool=true, + classification::Symbol=:threshold, + zero_winding::Int=-1, + winding_rfrac::Float64=0.5) length(Q) == length(Δ) || throw(ArgumentError("_extract_growth_rates_amr: length(Q) ≠ length(Δ)")) length(Q) >= 3 || @@ -1000,5 +1229,9 @@ function _extract_growth_rates_amr(Q::Vector{ComplexF64}, residual=residual, polish_maxit=polish_maxit, validity_scale=validity_scale, - validity_rtol=validity_rtol) + validity_rtol=validity_rtol, + polish=polish, + classification=classification, + zero_winding=zero_winding, + winding_rfrac=winding_rfrac) end diff --git a/src/Tearing/Runner/Control.jl b/src/Tearing/Runner/Control.jl index a35f3ef9c..5dd2c5ccc 100644 --- a/src/Tearing/Runner/Control.jl +++ b/src/Tearing/Runner/Control.jl @@ -132,6 +132,11 @@ constructor. - `filter_above_poles` -- discard roots above the highest pole γ - `filter_outside_re` -- condition the above-pole filter on the +γ step exiting the Re(Δ)=0 contour loop + - `root_classification` -- `:threshold` (default) classifies contour + intersections with the pole threshold and the filters above; `:winding` + classifies each by the winding number of the dispersion residual on a + small circle around it (argument principle), and the filters above decide + nothing. Applies to both coupling modes; see `find_growth_rates`. # Kinetic-profile source @@ -217,6 +222,9 @@ there is one consistent interface for resistive and kinetic profiles. # |residual| is not below validity_rtol × median(|Δ|) — a contour near-miss # / failed-Δ'-BVP surface, not a real root. Flagged `:spurious`. validity_rtol::Float64 = 1e-3 + # How contour intersections are classified as roots or poles: :threshold + # (pole_threshold + filters) or :winding (argument principle). + root_classification::Symbol = :threshold omega_E_kHz::Dict{String,Float64} = Dict{String,Float64}() @@ -236,6 +244,7 @@ const _VALID_D_C_TYPES = (:none, :lar, :fitzpatrick, :toroidal) const _VALID_P_PERP_MODELS = (:chi_perp_e, :chi_perp_i, :P_phi, :D_perp, :c_beta, :tau_E) const _VALID_RESISTIVITY_MODELS = (:sauter, :redl, :spitzer, :spitzer_harm) const _VALID_LNLAMBDA_FORMS = (:nrl, :sauter, :wesson) +const _VALID_ROOT_CLASSIFICATIONS = (:threshold, :winding) function validate(ctrl::SLAYERControl) ctrl.inner_model in _VALID_INNER_MODELS || ctrl.inner_model in _UNSUPPORTED_INNER_MODELS || @@ -269,6 +278,9 @@ function validate(ctrl::SLAYERControl) ctrl.lnLambda_form in _VALID_LNLAMBDA_FORMS || throw(ArgumentError("SLAYERControl: lnLambda_form=$(ctrl.lnLambda_form) " * "not in $(_VALID_LNLAMBDA_FORMS)")) + ctrl.root_classification in _VALID_ROOT_CLASSIFICATIONS || + throw(ArgumentError("SLAYERControl: root_classification=$(ctrl.root_classification) " * + "not in $(_VALID_ROOT_CLASSIFICATIONS)")) ctrl.impurity_charge > 1 || throw(ArgumentError("SLAYERControl: impurity_charge=$(ctrl.impurity_charge) must exceed the main-ion charge 1")) ctrl.msing_max >= 1 || @@ -317,6 +329,7 @@ function slayer_control_from_toml(section::AbstractDict) haskey(v, "pole_threshold") && (flat["pole_threshold"] = v["pole_threshold"]) haskey(v, "filter_above_poles") && (flat["filter_above_poles"] = v["filter_above_poles"]) haskey(v, "filter_outside_re") && (flat["filter_outside_re"] = v["filter_outside_re"]) + haskey(v, "root_classification") && (flat["root_classification"] = v["root_classification"]) elseif k == "CriticalResonantField" && v isa AbstractDict flat["critical_resonant_field"] = critical_resonant_field_control_from_toml(v) else @@ -338,7 +351,7 @@ function slayer_control_from_toml(section::AbstractDict) for (k, v) in flat sym = Symbol(k) if sym in (:inner_model, :scan_mode, :coupling_mode, :D_c_type, - :resistivity_model, :lnLambda_form, :P_perp_model) + :resistivity_model, :lnLambda_form, :P_perp_model, :root_classification) kwargs[sym] = v isa Symbol ? v : Symbol(String(v)) elseif sym in (:Q_re_range, :Q_im_range) kwargs[sym] = _as_range(v) diff --git a/src/Tearing/Runner/HDF5Output.jl b/src/Tearing/Runner/HDF5Output.jl index ab7914191..21bf90a06 100644 --- a/src/Tearing/Runner/HDF5Output.jl +++ b/src/Tearing/Runner/HDF5Output.jl @@ -13,7 +13,8 @@ # │ ├── psi, q, q1, ... # │ └── ... # ├── Roots/ -- complex Q_root, omega_Hz, gamma_Hz -# ├── Diagnostics/ -- ValidRoots, Poles, FilteredRoots (uncoupled: each surface's E×B frame) +# ├── Diagnostics/ -- ValidRoots, Poles, FilteredRoots (uncoupled: each surface's E×B frame; +# │ + *Winding when root_classification = :winding) # │ (flat-plus-offsets ragged encoding) # ├── Scan/ -- optional: full Q/Δ scan data # └── CriticalResonantField/ -- optional: torque-balance b_r/B_φ crit per surface @@ -184,6 +185,11 @@ const TEARING_RAGGED_H5_ANNOTATIONS = [ "offsets" => (; long_name="ragged-array offsets: row k spans offsets[k]+1:offsets[k+1]") ] +const TEARING_RAGGED_WINDING_H5_ANNOTATIONS = [ + "flat" => (; long_name="winding number of the dispersion residual around each entry of the matching root list (NaN: circle unresolved)"), + "offsets" => (; long_name="ragged-array offsets: row k spans offsets[k]+1:offsets[k+1]") +] + const TEARING_SCAN_H5_ANNOTATIONS = [ "kind" => (; long_name="scan kind: brute_force or amr"), "Q" => (; long_name="sampled complex normalized frequency Q"), @@ -200,7 +206,8 @@ function _annotate_tearing!(g) ann.annotate!(g, TEARING_H5_ANNOTATIONS) if haskey(g, "Diagnostics") for sub in keys(g["Diagnostics"]) - ann.annotate!(g["Diagnostics"][sub], TEARING_RAGGED_H5_ANNOTATIONS) + ann.annotate!(g["Diagnostics"][sub], + endswith(sub, "Winding") ? TEARING_RAGGED_WINDING_H5_ANNOTATIONS : TEARING_RAGGED_H5_ANNOTATIONS) end end if haskey(g, "CriticalResonantField/Scan") @@ -319,22 +326,30 @@ function _write_diagnostics!(g, r::SLAYERResult) r.per_surface_extraction end - _write_ragged_complex!(diag, "ValidRoots", + _write_ragged!(diag, "ValidRoots", [gr.valid_roots for gr in extractions]) - _write_ragged_complex!(diag, "Poles", + _write_ragged!(diag, "Poles", [gr.poles for gr in extractions]) - _write_ragged_complex!(diag, "FilteredRoots", + _write_ragged!(diag, "FilteredRoots", [gr.filtered_roots for gr in extractions]) + if any(gr -> gr.classification === :winding, extractions) + _write_ragged!(diag, "ValidRootsWinding", + [gr.valid_roots_winding for gr in extractions]) + _write_ragged!(diag, "PolesWinding", + [gr.poles_winding for gr in extractions]) + _write_ragged!(diag, "FilteredRootsWinding", + [gr.filtered_roots_winding for gr in extractions]) + end return nothing end -# Write a ragged vector-of-vectors of ComplexF64 as (flat, offsets) — +# Write a ragged vector-of-vectors (ComplexF64 or Float64) as (flat, offsets) — # `offsets[k+1] - offsets[k]` is the length of row `k`. This avoids HDF5 VLEN # types, which have patchy cross-language support. -function _write_ragged_complex!(parent, name::String, - data::Vector{Vector{ComplexF64}}) +function _write_ragged!(parent, name::String, + data::Vector{Vector{T}}) where {T<:Union{ComplexF64,Float64}} g = create_group(parent, name) - flat = ComplexF64[] + flat = T[] offsets = Int[0] for v in data append!(flat, v) diff --git a/src/Tearing/Runner/run_slayer.jl b/src/Tearing/Runner/run_slayer.jl index cf07e5954..80baab4a7 100644 --- a/src/Tearing/Runner/run_slayer.jl +++ b/src/Tearing/Runner/run_slayer.jl @@ -133,6 +133,13 @@ _surface_label(p::InnerLayerParameters) = "ising=$(getfield(p, :ising))" _is_ggj(::GGJModel) = true _is_ggj(::Any) = false +# Winding number of an isolated zero of the dispersion residual, for +# root_classification = :winding. The SLAYER residual is a function of conj(Q) +# (solve_inner evaluates conj(Δ̂(ĝ)); see Riccati.jl), so a zero winds −1. +_zero_winding(::SLAYERModel) = -1 +_zero_winding(model) = + throw(ArgumentError("root_classification=:winding is not available for inner model $(typeof(model))")) + # --------------------------------------------------------------------- # Scan dispatch # --------------------------------------------------------------------- @@ -425,6 +432,9 @@ function run_slayer_from_inputs(model::InnerLayer.InnerLayerModel, return 10.0 * median(abs.(finite)) end + winding = control.root_classification === :winding + zw = winding ? _zero_winding(model) : -1 + if control.coupling_mode === :uncoupled for (k, sc) in enumerate(scs) scan = _run_scan(sc, control) @@ -434,8 +444,11 @@ function run_slayer_from_inputs(model::InnerLayer.InnerLayerModel, filter_above_poles=control.filter_above_poles, filter_outside_re=control.filter_outside_re, gap_kHz_threshold=control.gap_kHz_threshold, - residual=control.polish_roots ? sc : nothing, - validity_rtol=control.validity_rtol) + residual=(control.polish_roots || winding) ? sc : nothing, + polish=control.polish_roots, + validity_rtol=control.validity_rtol, + classification=control.root_classification, + zero_winding=zw) :no_root in gr.warning_flags && @warn( "SLAYER: no usable growth-rate root found for surface " * "$(_surface_label(params[k])); reported γ=0 is a " * @@ -460,12 +473,18 @@ function run_slayer_from_inputs(model::InnerLayer.InnerLayerModel, # determinant det(D'−D(Q)) is ill-conditioned (its magnitude floors well # above zero), so |det|-based polishing is a no-op and the residual-scale # gate is unreliable. The coupled root is therefore the raw contour - # intersection and carries the triangulation's interpolation error. + # intersection and carries the triangulation's interpolation error. With + # root_classification = :winding the determinant is evaluated around + # each intersection only to count its zeros and poles. gr = find_growth_rates(scan, ref_tauk; pole_threshold=pthr, filter_above_poles=control.filter_above_poles, filter_outside_re=control.filter_outside_re, - gap_kHz_threshold=control.gap_kHz_threshold) + gap_kHz_threshold=control.gap_kHz_threshold, + residual=winding ? mc : nothing, + polish=false, + classification=control.root_classification, + zero_winding=zw) :no_root in gr.warning_flags && @warn( "SLAYER: no usable growth-rate root found for the coupled " * "$(m_use)-surface determinant; reported γ=0 is a placeholder, " * diff --git a/test/runtests.jl b/test/runtests.jl index c52bbf8ab..8f93b8c23 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -64,6 +64,7 @@ else include("./runtests_dispersion_amr.jl") include("./runtests_dispersion_polish.jl") include("./runtests_slayer_profile_transport.jl") + include("./runtests_dispersion_winding.jl") include("./runtests_slayer_runner.jl") include("./runtests_critical_resonant_field.jl") include("./runtests_kinetic.jl") diff --git a/test/runtests_dispersion_winding.jl b/test/runtests_dispersion_winding.jl new file mode 100644 index 000000000..2a961bd0b --- /dev/null +++ b/test/runtests_dispersion_winding.jl @@ -0,0 +1,119 @@ +@testset "Dispersion winding-number root classification" begin + using GeneralizedPerturbedEquilibrium.Dispersion + D = GeneralizedPerturbedEquilibrium.Dispersion + + # ---------------------------------------------------------------- + # _winding_number: the argument principle on synthetic functions + # ---------------------------------------------------------------- + @testset "_winding_number: zeros, poles, both orientations" begin + z = 0.3 - 0.2im + @test D._winding_number(Q -> Q - z, z, 0.1)[1] ≈ 1 atol = 1e-10 + @test D._winding_number(Q -> 1 / (Q - z), z, 0.1)[1] ≈ -1 atol = 1e-10 + # a function of conj(Q): its zero winds -1, its pole +1 + @test D._winding_number(Q -> conj(Q - z), z, 0.1)[1] ≈ -1 atol = 1e-10 + @test D._winding_number(Q -> 1 / conj(Q - z), z, 0.1)[1] ≈ 1 atol = 1e-10 + # nothing inside the circle + @test D._winding_number(Q -> Q - z, z + 1.0, 0.1)[1] ≈ 0 atol = 1e-10 + # two zeros inside, and a zero-pole pair inside + @test D._winding_number(Q -> (Q - z) * (Q - z - 0.02), z, 0.1)[1] ≈ 2 atol = 1e-10 + @test D._winding_number(Q -> (Q - z) / (Q - z - 0.02), z, 0.1)[1] ≈ 0 atol = 1e-10 + end + + @testset "_winding_number: refines the sampling near the circle" begin + # zero just inside the circle: the phase turns fast where the circle + # passes it, so the default 64 samples are not enough + z = 0.0 + 0.0im + W, step = D._winding_number(Q -> Q - z, z + 0.0998, 0.1) + @test W ≈ 1 atol = 1e-10 + @test step < π / 2 + # zero exactly on the circle → unresolved + W, _ = D._winding_number(Q -> Q - z, z + 0.1, 0.1) + @test isnan(W) + # non-finite residual → unresolved + W, _ = D._winding_number(Q -> ComplexF64(NaN, NaN), z, 0.1) + @test isnan(W) + end + + # ---------------------------------------------------------------- + # find_growth_rates(...; classification=:winding) + # ---------------------------------------------------------------- + # Two zeros and a pole between them, off the scan nodes; the upper zero lies + # ABOVE the pole. + z_lo, z_hi, p = -0.4137 - 0.3021im, 0.3519 + 0.4473im, 0.0531 + 0.1077im + g(Q) = (ComplexF64(Q) - z_lo) * (ComplexF64(Q) - z_hi) / (ComplexF64(Q) - p) + fanti(Q) = conj(g(Q)) # SLAYER-like: a function of conj(Q) + + @testset "brute-force scan: roots and pole by winding number (anti-holomorphic)" begin + scan = brute_force_scan(fanti, (-1.0, 1.0), (-1.0, 1.0); nre=81, nim=81, threaded=false) + r = find_growth_rates(scan, 1.0; residual=fanti, classification=:winding, zero_winding=-1) + @test r.classification === :winding + @test length(r.valid_roots) == 2 + @test length(r.poles) == 1 + @test all(w -> w ≈ -1, r.valid_roots_winding) + @test all(w -> w ≈ 1, r.poles_winding) + @test length(r.valid_roots_winding) == length(r.valid_roots) + @test length(r.poles_winding) == length(r.poles) + @test minimum(abs.(r.valid_roots .- z_hi)) < 1e-8 # polished + @test minimum(abs.(r.valid_roots .- z_lo)) < 1e-8 + @test abs(r.poles[1] - p) < 0.05 + # the least stable zero is reported even though it lies above the pole + @test abs(r.Q_root - z_hi) < 1e-8 + @test r.gamma_Hz ≈ imag(z_hi) atol = 1e-8 + @test isnan(real(r.Q_root_secondary)) + @test !(:no_root in r.warning_flags) + end + + @testset "holomorphic residual: zero_winding = +1" begin + scan = brute_force_scan(g, (-1.0, 1.0), (-1.0, 1.0); nre=81, nim=81, threaded=false) + r = find_growth_rates(scan, 1.0; residual=g, classification=:winding, zero_winding=1) + @test length(r.valid_roots) == 2 && length(r.poles) == 1 + @test abs(r.Q_root - z_hi) < 1e-8 + # the wrong orientation swaps roots and poles: the convention matters + w = find_growth_rates(scan, 1.0; residual=g, classification=:winding, zero_winding=-1) + @test length(w.valid_roots) == 1 && length(w.poles) == 2 + end + + @testset "polish = false keeps the contour estimate" begin + scan = brute_force_scan(fanti, (-1.0, 1.0), (-1.0, 1.0); nre=41, nim=41, threaded=false) + r = find_growth_rates(scan, 1.0; residual=fanti, classification=:winding, polish=false) + @test length(r.valid_roots) == 2 + @test 0 < abs(r.Q_root - z_hi) < 0.05 + end + + @testset "AMR scan" begin + amr = amr_scan(fanti, (-1.0, 1.0), (-1.0, 1.0); nre0=12, nim0=12, passes=3) + r = find_growth_rates(amr, 1.0; residual=fanti, classification=:winding) + @test r.classification === :winding + @test length(r.valid_roots) == 2 + @test length(r.poles) == 1 + @test abs(r.Q_root - z_hi) < 1e-8 + end + + @testset "an intersection with no zero inside is filtered as spurious" begin + re_paths = [ComplexF64[-1.0, 1.0]] + im_paths = [ComplexF64[-1.0im, 1.0im]] + r = D._run_analysis(re_paths, im_paths, [Float64[0.0, 0.0]], 1.0; + pole_threshold=10.0, filter_above_poles=true, filter_outside_re=true, + residual=Q -> ComplexF64(1.0), classification=:winding, zero_winding=-1) + @test isempty(r.valid_roots) && isempty(r.poles) + @test r.filtered_roots == [0.0 + 0.0im] + @test r.filtered_roots_winding ≈ [0.0] + @test :spurious in r.warning_flags + @test :no_root in r.warning_flags + @test r.gamma_Hz == 0.0 + end + + @testset "threshold mode is unchanged and carries no winding numbers" begin + scan = brute_force_scan(fanti, (-1.0, 1.0), (-1.0, 1.0); nre=41, nim=41, threaded=false) + r = find_growth_rates(scan, 1.0) + @test r.classification === :threshold + @test isempty(r.valid_roots_winding) && isempty(r.poles_winding) && isempty(r.filtered_roots_winding) + end + + @testset "argument validation" begin + scan = brute_force_scan(fanti, (-1.0, 1.0), (-1.0, 1.0); nre=11, nim=11, threaded=false) + @test_throws ArgumentError find_growth_rates(scan, 1.0; classification=:winding) + @test_throws ArgumentError find_growth_rates(scan, 1.0; residual=fanti, classification=:winding, zero_winding=0) + @test_throws ArgumentError find_growth_rates(scan, 1.0; classification=:bogus) + end +end diff --git a/test/runtests_slayer_runner.jl b/test/runtests_slayer_runner.jl index 2d9afd6e9..10d58bc87 100644 --- a/test/runtests_slayer_runner.jl +++ b/test/runtests_slayer_runner.jl @@ -45,6 +45,12 @@ SLAYERControl(; msing_max=0)) @test_throws ArgumentError Runner.validate( SLAYERControl(; nre=1)) + @test c.root_classification === :threshold + @test_throws ArgumentError Runner.validate( + SLAYERControl(; root_classification=:bogus)) + @test slayer_control_from_toml(Dict{String,Any}("root_classification" => "winding")).root_classification === :winding + @test slayer_control_from_toml(Dict{String,Any}("growth_rate_filter" => + Dict{String,Any}("root_classification" => "winding"))).root_classification === :winding end @testset "slayer_control_from_toml: nested sections flatten" begin @@ -398,6 +404,7 @@ # Ragged diagnostics use flat+offsets encoding @test haskey(g["Diagnostics/ValidRoots"], "flat") @test haskey(g["Diagnostics/ValidRoots"], "offsets") + @test !haskey(g["Diagnostics"], "ValidRootsWinding") # threshold mode # Scan group present (store_scan=true) @test haskey(g, "Scan/Surface_1") @@ -406,6 +413,61 @@ end end + @testset "root_classification = :winding: uncoupled + coupled, HDF5" begin + p1 = _mk_params(; rs=0.5, lu=1.0e7, tauk=1.0e-4, m=2, ising=1) + p2 = _mk_params(; rs=0.6, lu=2.0e7, tauk=1.2e-4, m=3, ising=2) + model = SLAYERModel(; variant=:fitzpatrick) + # Diagonal Δ' placing surface 1's zero at Q1 and surface 2's at Q2 (each + # in its own normalization), so the coupled determinant has two + # distinct zeros. Both points are off the scan nodes and away from Q = 0. + Q1, Q2 = -0.1213 + 0.0537im, 0.1517 + 0.1029im + Δ1 = InnerLayer.solve_inner(model, p1, Q1).tearing * p1.lu^(1 / 3) + Δ2 = InnerLayer.solve_inner(model, p2, Q2).tearing * p2.lu^(1 / 3) + dp = ComplexF64[Δ1 0.0; 0.0 Δ2] + for mode in (:uncoupled, :coupled) + ctrl(cls) = SLAYERControl(; enabled=true, scan_mode=:brute_force, coupling_mode=mode, + Q_re_range=(-0.5, 0.5), Q_im_range=(-0.3, 0.3), nre=40, nim=40, + root_classification=cls, pole_threshold_adaptive=true) # adaptive: threshold-mode reference + r = run_slayer_from_inputs(model, [p1, p2], dp, ctrl(:winding); + rational_psi=[0.45, 0.72], rational_q=[2.0, 3.0]) + ref = run_slayer_from_inputs(model, [p1, p2], dp, ctrl(:threshold)) + exts = mode === :coupled ? [r.coupled_extraction] : r.per_surface_extraction + @test all(gr -> gr.classification === :winding, exts) + @test all(gr -> all(w -> w ≈ -1, gr.valid_roots_winding), exts) + if mode === :uncoupled + # Polished. The SLAYER residual is only smooth to |r| ~ 1e-2 within + # 1e-4 of these zeros, which bounds any polish; the grid step is 0.026. + @test abs(r.Q_root[1] - Q1) < 1e-3 + @test abs(r.Q_root[2] - Q2) < 1e-3 + @test maximum(abs.(r.Q_root .- ref.Q_root)) < 1e-3 + else + # zeros at Q1 and Q2·τ_k2/τ_k1 in the reference normalization + # (surface k sees Q·τ_k1/τ_kk); the least stable is reported, + # unpolished since the coupled path does not polish + @test length(r.coupled_extraction.valid_roots) == 2 + @test abs(r.Q_root[1] - Q2 * p2.tauk / p1.tauk) < 2e-2 # grid-resolution limited + @test r.Q_root == ref.Q_root + end + mktemp() do path, io + close(io) + h5open(path, "w") do f + write_slayer_hdf5!(f, r) + end + h5open(path, "r") do f + viol = _collect_metadata_violations(f) + isempty(viol) || @error "metadata contract violations in Tearing/" viol + @test isempty(viol) + g = f["Tearing/Diagnostics"] + for k in ("ValidRootsWinding", "PolesWinding", "FilteredRootsWinding") + @test haskey(g, k) + end + @test read(g["ValidRootsWinding/offsets"]) == read(g["ValidRoots/offsets"]) + @test all(w -> w ≈ -1, read(g["ValidRootsWinding/flat"])) + end + end + end + end + @testset "GGJ per-surface writer: metadata contract" begin # The GGJ branch of _write_per_surface! is unreachable from the SLAYER # round-trip above; enforce its annotation-table rows on a synthetic write. From c7057033d78b85c5126f74874e996041acc594a6 Mon Sep 17 00:00:00 2001 From: d-burg Date: Wed, 28 Oct 2026 21:58:37 -0400 Subject: [PATCH 2/5] Tearing - REFACTOR - Plug winding classification into the shared root analysis and harden its circles Co-Authored-By: Claude Opus 5.5 --- .../Dispersion/GrowthRateExtraction.jl | 191 ++++++++++-------- src/Tearing/Runner/run_slayer.jl | 2 +- test/runtests_dispersion_winding.jl | 46 +++++ test/runtests_slayer_runner.jl | 10 + 4 files changed, 161 insertions(+), 88 deletions(-) diff --git a/src/Tearing/Dispersion/GrowthRateExtraction.jl b/src/Tearing/Dispersion/GrowthRateExtraction.jl index db8582972..9d5da441a 100644 --- a/src/Tearing/Dispersion/GrowthRateExtraction.jl +++ b/src/Tearing/Dispersion/GrowthRateExtraction.jl @@ -20,10 +20,10 @@ # pole are filtered out (spurious upper branches). # 5. Return the highest-γ surviving root in physical units. # -# With `classification = :winding`, steps 3-4 are replaced by the argument +# With `classification = :winding`, step 3 is replaced by the argument # principle: each intersection is a root, a pole or neither according to the # winding number of the residual on a small circle around it (see -# `_run_analysis_winding`). +# `_winding_classify`); the rest of the analysis is shared. using Contour using DelaunayTriangulation @@ -55,10 +55,14 @@ Output of `find_growth_rates`. | | `NaN`; `omega_Hz`/`gamma_Hz` fall back to 0 — check | | | this flag to tell that apart from a true γ≈0 result). | | | Empty if root is clean. | -| `valid_roots` | All non-pole intersections that survived pole filter | +| `valid_roots` | Intersections classified as roots (threshold: not a | +| | pole and passed the validity gate; winding: a zero | +| | count of the sign of `zero_winding`) | | `poles` | Intersections classified as poles | -| `filtered_roots` | Intersections rejected by the above-pole/outside-Re | -| | filter or the new geom+gap recursion | +| `filtered_roots` | Gate-rejected near-misses and roots rejected by the | +| | above-pole/outside-Re filter (threshold), or | +| | intersections with a zero or unresolved count | +| | (winding) | | `re_contours` | Extracted Re(Δ)=`re_target` polylines | | `im_contours` | Extracted Im(Δ)=`im_target` polylines | | `pole_threshold` | Threshold used for pole classification | @@ -156,8 +160,8 @@ single-surface scans; `mc.surfaces[mc.ref_idx].tauk` for coupled scans). that is anti-holomorphic (a function of conj(Q)), such as the SLAYER residual. Default `-1`. - `winding_rfrac` -- `:winding` mode only: radius of the classifying circle - in units of the median contour-segment length (default 0.5), capped at - 0.45 × the distance to the nearest other intersection. + in units of the local contour-segment length at the intersection (default + 0.5), capped at 0.45 × the distance to the nearest other intersection. # Winding classification (`classification = :winding`) @@ -172,14 +176,20 @@ rather than estimated. An intersection is a root when the count has the sign of `filtered_roots` (flag `:spurious`) when the count is zero, e.g. a contour near-miss or an unresolved pole-zero pair. A zero count is checked once more on a circle eight times wider (still clear of every other intersection), because -the contour estimate of a pole can lie outside the small circle. A circle that -stays unresolved (winding `NaN`, a zero or pole on it) is filtered too and flags -`:winding_unresolved`. `pole_threshold`, `filter_above_poles`, -`filter_outside_re`, `gap_kHz_threshold` and `validity_rtol` decide nothing in -this mode, and `Q_root` is the least stable root, with no secondary root. -With `polish = true` a root is polished inside its classifying circle, which is -proven to contain it. `residual` is required. Each intersection costs between -64 and a few hundred residual evaluations. Only zeros that the scan resolves +the contour estimate of a pole can lie outside the small circle. A non-finite +circle (the residual failed on it) is retried at half and a quarter of the +radius, and a zero count must agree on the half-radius circle; a circle that +stays unresolved is filtered too and flags `:winding_unresolved`. The count is +exact only for a residual that is single-valued and meromorphic on the circle: +a branch cut crossing it, or a jump in the residual, can change the count. +The classification replaces the pole threshold, the validity gate and the +above-pole/outside-Re filter, so `pole_threshold`, `validity_rtol`, +`filter_above_poles` and `filter_outside_re` decide nothing in this mode; the +selection of `Q_root`, the `:geom`/`:gap` warnings and the secondary root are +the same as in `:threshold` mode. With `polish = true` a root is polished inside +its classifying circle, which is proven to contain it. `residual` is required. +Each intersection costs between about 130 and a few thousand residual +evaluations (two circles of 64–4096 samples). Only zeros that the scan resolves as contour crossings are found, as in `:threshold` mode. # Spurious-root recursion @@ -666,35 +676,55 @@ function _run_analysis(re_paths::Vector{Vector{ComplexF64}}, classification::Symbol=:threshold, zero_winding::Int=-1, winding_rfrac::Float64=0.5) - classification === :winding && - return _run_analysis_winding(re_paths, im_paths, tauk, residual; - zero_winding=zero_winding, rfrac=winding_rfrac, - polish=polish, polish_maxit=polish_maxit, pole_threshold=pole_threshold) + winding = classification === :winding raw_intersections = _all_intersections(re_paths, im_paths) poles = ComplexF64[] candidates = Tuple{ComplexF64,Bool}[] # (pt, on_top_half_re_flag) - spurious_roots = ComplexF64[] # polished, but |residual| ≫ 0 + spurious_roots = ComplexF64[] # polished, but |residual| ≫ 0, or a zero winding count + w_valid, w_poles, w_filtered = Float64[], Float64[], Float64[] + warning_flags = Symbol[] # Isolate-then-polish: when a residual callable is supplied, refine each # (non-pole) contour intersection to the true zero of `f`. Trust radii are # capped by the spacing to neighbouring intersections so closely-spaced # coupled roots stay coherently bound to their own contour crossing. do_polish = polish && residual !== nothing - h_grid = do_polish ? _median_segment_length(re_paths, im_paths) : 0.0 + h_grid = (do_polish || winding) ? _median_segment_length(re_paths, im_paths) : 0.0 + # Contours crossing at a shared vertex are reported once per adjacent segment pair; a winding count needs each once. + winding && (raw_intersections = _dedupe_points(raw_intersections, 1e-9 * max(h_grid, eps()))) # Validity gate: a genuine root drives |residual| ≈ 0; a spurious contour # near-miss (e.g. a surface whose Δ' BVP failed → |Δ'| huge and # uncancellable by the inner layer) leaves |residual| stuck near the typical # scan magnitude. Only active when polishing AND a scale is supplied. - do_gate = do_polish && validity_scale > 0 + # A winding count already proves a zero, so the gate applies to the threshold classification only. + do_gate = do_polish && validity_scale > 0 && !winding polish_evals = 0 for (i_pt, pt) in enumerate(raw_intersections) - # --- 1. classify as pole or root via local Re-magnitude on Im contour best_im_path_idx, best_im_vert_idx, _ = _closest_polyline_vertex(im_paths, pt) + # --- 1. classify as root or pole: by winding number, or by the local Re-magnitude on the Im contour + R_root = 0.0 # trust radius for the polish; the classifying circle in :winding mode + if winding + kind, W, r = _winding_classify(residual, raw_intersections, i_pt, re_paths, im_paths, h_grid; + zero_winding=zero_winding, rfrac=winding_rfrac) + if kind === :pole + push!(poles, pt) + push!(w_poles, W) + continue + elseif kind !== :root + push!(spurious_roots, pt) + push!(w_filtered, W) + flag = kind === :unresolved ? :winding_unresolved : :spurious + flag in warning_flags || push!(warning_flags, flag) + continue + end + push!(w_valid, W) + R_root = r + end is_pole = false - if best_im_path_idx > 0 + if !winding && best_im_path_idx > 0 re_vals = im_re_vals[best_im_path_idx] n = length(re_vals) i_prev = max(1, best_im_vert_idx - 1) @@ -713,7 +743,7 @@ function _run_analysis(re_paths::Vector{Vector{ComplexF64}}, # --- 1b. polish this (non-pole) intersection to the true zero of `f`, # bounded to a neighbour-aware trust region (no-op on failure). if do_polish - R = _polish_trust_radius(raw_intersections, i_pt, h_grid) + R = winding ? R_root : _polish_trust_radius(raw_intersections, i_pt, h_grid) pt, nev, _ = _polish_root(residual, pt, R; maxit=polish_maxit) polish_evals += nev # --- 1c. validity gate: drop intersections that don't polish to a @@ -777,8 +807,7 @@ function _run_analysis(re_paths::Vector{Vector{ComplexF64}}, filtered_roots = copy(spurious_roots) # gate-rejected near-misses kept here Q_root = ComplexF64(NaN, NaN) Q_root_2nd = ComplexF64(NaN, NaN) - warning_flags = Symbol[] - isempty(spurious_roots) || push!(warning_flags, :spurious) + !winding && !isempty(spurious_roots) && push!(warning_flags, :spurious) if !isempty(valid_roots) order = sortperm(valid_roots; by=q -> -imag(q)) @@ -791,8 +820,8 @@ function _run_analysis(re_paths::Vector{Vector{ComplexF64}}, for k in 1:length(sorted_pts) cand = sorted_pts[k] top_re = sorted_top[k] - # Legacy filter: above-pole + closed-loop outside-Re - legacy_reject = filter_above_poles && imag(cand) > max_pole_gamma && + # Legacy filter: above-pole + closed-loop outside-Re (a winding-classified root is a proven zero) + legacy_reject = !winding && filter_above_poles && imag(cand) > max_pole_gamma && (!filter_outside_re || top_re) if legacy_reject push!(filtered_roots, cand) @@ -850,7 +879,8 @@ function _run_analysis(re_paths::Vector{Vector{ComplexF64}}, Q_root_2nd, omega_Hz_2nd, gamma_Hz_2nd, warning_flags, valid_roots, poles, filtered_roots, - re_paths, im_paths, pole_threshold) + re_paths, im_paths, pole_threshold, + classification, w_valid, w_poles, w_filtered) end # --------------------------------------------------------------------- @@ -921,68 +951,55 @@ function _dedupe_points(pts::Vector{ComplexF64}, tol::Float64) return out end -function _run_analysis_winding(re_paths::Vector{Vector{ComplexF64}}, - im_paths::Vector{Vector{ComplexF64}}, - tauk::Float64, residual; - zero_winding::Int, rfrac::Float64, polish::Bool, polish_maxit::Int, - pole_threshold::Float64) - h_grid = _median_segment_length(re_paths, im_paths) - pts = _dedupe_points(_all_intersections(re_paths, im_paths), 1e-9 * max(h_grid, eps())) - - valid_roots, poles, filtered = ComplexF64[], ComplexF64[], ComplexF64[] - w_valid, w_poles, w_filtered = Float64[], Float64[], Float64[] - warning_flags = Symbol[] +# Length of the contour segments that meet at the vertices nearest `pt`, the local resolution there; on an AMR mesh it +# follows the cell size, unlike the median over the whole scan. +function _local_segment_length(re_paths::Vector{Vector{ComplexF64}}, im_paths::Vector{Vector{ComplexF64}}, pt::ComplexF64, + fallback::Float64) + h = 0.0 + for paths in (re_paths, im_paths) + ip, iv, _ = _closest_polyline_vertex(paths, pt) + ip > 0 || continue + path = paths[ip] + iv > 1 && (h = max(h, abs(path[iv] - path[iv-1]))) + iv < length(path) && (h = max(h, abs(path[iv+1] - path[iv]))) + end + return h > 0 ? h : fallback +end - for (i, pt) in enumerate(pts) - # circle radius: a fraction of the local resolution, and less than half - # the distance to any other intersection so circles never overlap - r = _polish_trust_radius(pts, i, h_grid; k_cell=rfrac, nn_frac=0.45) +# Classify intersection `i` of `pts` by the winding number of `residual` on a circle around it: `(kind, W, r)` with kind +# `:root`, `:pole`, `:none` (zero count) or `:unresolved`, the count W and the circle radius r. The circle is a fraction +# `rfrac` of the local contour resolution and under half the distance to every other intersection; a non-finite circle is +# retried at r/2 and r/4. A zero count must agree on the half-radius circle (the contour estimate of a zero is accurate, so a +# disagreement means a second zero or a pole near the circle), and a zero count is checked once on a circle eight times +# wider, because the contour estimate of a pole can lie outside the small one. +function _winding_classify(residual, pts::Vector{ComplexF64}, i::Int, re_paths, im_paths, h_grid::Float64; + zero_winding::Int, rfrac::Float64) + pt = pts[i] + h_loc = _local_segment_length(re_paths, im_paths, pt, h_grid) + r = _polish_trust_radius(pts, i, h_loc; k_cell=rfrac, nn_frac=0.45) + W = NaN + for _ in 1:3 W, _ = _winding_number(residual, pt, r) - # The contour estimate of a zero is accurate (the field is nearly linear - # there); that of a pole is not, and can fall outside the small circle. - # If it counts nothing, count once more on a wider circle, still clear - # of every other intersection. - if !(isfinite(W) && round(Int, W) != 0) - r_wide = _polish_trust_radius(pts, i, h_grid; k_cell=8rfrac, nn_frac=0.45) - if r_wide > r - W_wide, _ = _winding_number(residual, pt, r_wide) - if isfinite(W_wide) && round(Int, W_wide) != 0 - W, r = W_wide, r_wide - end + isfinite(W) && break + r /= 2 + end + isfinite(W) || return (:unresolved, NaN, r) + n = round(Int, W) + if n * zero_winding >= 1 + W_half, _ = _winding_number(residual, pt, r / 2) + (isfinite(W_half) && round(Int, W_half) == n) || return (:unresolved, NaN, r) + elseif n == 0 + r_wide = _polish_trust_radius(pts, i, h_loc; k_cell=8rfrac, nn_frac=0.45) + if r_wide > r + W_wide, _ = _winding_number(residual, pt, r_wide) + if isfinite(W_wide) && round(Int, W_wide) != 0 + W, r, n = W_wide, r_wide, round(Int, W_wide) end end - n = isfinite(W) ? round(Int, W) : 0 - if !isfinite(W) || abs(W - n) > 0.1 - push!(filtered, pt) - push!(w_filtered, NaN) - :winding_unresolved in warning_flags || push!(warning_flags, :winding_unresolved) - elseif n * zero_winding >= 1 # net zeros inside the circle - if polish - pt, _, _ = _polish_root(residual, pt, r; maxit=polish_maxit) - end - push!(valid_roots, pt) - push!(w_valid, W) - elseif n * zero_winding <= -1 # net poles inside - push!(poles, pt) - push!(w_poles, W) - else # zero count: near-miss or pole-zero pair - push!(filtered, pt) - push!(w_filtered, W) - :spurious in warning_flags || push!(warning_flags, :spurious) - end end - - Q_root = isempty(valid_roots) ? ComplexF64(NaN, NaN) : valid_roots[argmax(imag.(valid_roots))] - isnan(real(Q_root)) && push!(warning_flags, :no_root) - omega_Hz = isnan(real(Q_root)) ? 0.0 : real(Q_root) / tauk - gamma_Hz = isnan(imag(Q_root)) ? 0.0 : imag(Q_root) / tauk - - return GrowthRateResult(Q_root, omega_Hz, gamma_Hz, - ComplexF64(NaN, NaN), 0.0, 0.0, - warning_flags, - valid_roots, poles, filtered, - re_paths, im_paths, pole_threshold, - :winding, w_valid, w_poles, w_filtered) + n * zero_winding >= 1 && return (:root, W, r) + n * zero_winding <= -1 && return (:pole, W, r) + return (:none, W, r) end # Regular-grid path: extract contours via Contour.jl, compute im_re_vals by diff --git a/src/Tearing/Runner/run_slayer.jl b/src/Tearing/Runner/run_slayer.jl index 80baab4a7..df9d76be7 100644 --- a/src/Tearing/Runner/run_slayer.jl +++ b/src/Tearing/Runner/run_slayer.jl @@ -136,7 +136,7 @@ _is_ggj(::Any) = false # Winding number of an isolated zero of the dispersion residual, for # root_classification = :winding. The SLAYER residual is a function of conj(Q) # (solve_inner evaluates conj(Δ̂(ĝ)); see Riccati.jl), so a zero winds −1. -_zero_winding(::SLAYERModel) = -1 +_zero_winding(::SLAYERModel{:fitzpatrick}) = -1 _zero_winding(model) = throw(ArgumentError("root_classification=:winding is not available for inner model $(typeof(model))")) diff --git a/test/runtests_dispersion_winding.jl b/test/runtests_dispersion_winding.jl index 2a961bd0b..b05ba6803 100644 --- a/test/runtests_dispersion_winding.jl +++ b/test/runtests_dispersion_winding.jl @@ -103,6 +103,52 @@ @test r.gamma_Hz == 0.0 end + @testset "a root about one cell from a pole, with the validity gate on" begin + # Zero 0.04 from a pole on a 0.025 grid: the pole's field leaks into the root's neighborhood. + zr, pp = 0.2131 + 0.1517im, 0.2131 + 0.1917im + h(Q) = conj((ComplexF64(Q) - zr) / (ComplexF64(Q) - pp)) + scan = brute_force_scan(h, (-0.5, 0.5), (-0.5, 0.5); nre=41, nim=41, threaded=false) + r = find_growth_rates(scan, 1.0; residual=h, classification=:winding, zero_winding=-1) + @test length(r.valid_roots) == 1 && length(r.poles) == 1 + @test abs(r.Q_root - zr) < 1e-8 + @test abs(r.poles[1] - pp) < 0.03 + end + + @testset "a residual that fails near a pole: the circle is retried smaller" begin + # The residual is NaN within 0.06 of the pole, as a failed inner-layer solve would be; the zero sits just outside. + zr, pp = -0.1173 + 0.0811im, -0.1173 + 0.1611im + h(Q) = abs(ComplexF64(Q) - pp) < 0.06 ? ComplexF64(NaN, NaN) : conj((ComplexF64(Q) - zr) / (ComplexF64(Q) - pp)) + scan = brute_force_scan(h, (-0.5, 0.5), (-0.5, 0.5); nre=21, nim=21, threaded=false) + r = find_growth_rates(scan, 1.0; residual=h, classification=:winding, zero_winding=-1) + @test any(z -> abs(z - zr) < 1e-8, r.valid_roots) + @test abs(r.Q_root - zr) < 1e-8 + end + + @testset "a coupled two-surface determinant with an off-diagonal term" begin + # det of [d1 − a(Q) c; c d2 − b(Q)] with a, b functions of conj(Q) and real d1, d2, c is itself a function of + # conj(Q); its zeros solve (d1 − Q̄)(d2 − 2Q̄) = c² and are not those of either diagonal entry alone. + d1, d2, c = 0.31 + 0.12im, 0.18 - 0.05im, 0.07 + det2(Q) = conj((d1 - ComplexF64(Q)) * (d2 - 2ComplexF64(Q)) - c^2) + # (d1 − x)(d2 − 2x) − c² = 0 for x = Q̄ + A, B, C = 2.0, -(d2 + 2d1), d1 * d2 - c^2 + xs = [(-B + sqrt(B^2 - 4A * C)) / 2A, (-B - sqrt(B^2 - 4A * C)) / 2A] + zs = conj.(xs) + scan = brute_force_scan(det2, (-0.5, 0.8), (-0.6, 0.6); nre=53, nim=49, threaded=false) + r = find_growth_rates(scan, 1.0; residual=det2, classification=:winding, zero_winding=-1) + @test length(r.valid_roots) == 2 + for z in zs + @test any(v -> abs(v - z) < 1e-8, r.valid_roots) + end + @test all(w -> w ≈ -1, r.valid_roots_winding) + end + + @testset "coarse AMR first pass: circles follow the local resolution" begin + amr = amr_scan(fanti, (-1.0, 1.0), (-1.0, 1.0); nre0=6, nim0=6, passes=2) + r = find_growth_rates(amr, 1.0; residual=fanti, classification=:winding) + @test length(r.valid_roots) == 2 && length(r.poles) == 1 + @test abs(r.Q_root - z_hi) < 1e-8 + end + @testset "threshold mode is unchanged and carries no winding numbers" begin scan = brute_force_scan(fanti, (-1.0, 1.0), (-1.0, 1.0); nre=41, nim=41, threaded=false) r = find_growth_rates(scan, 1.0) diff --git a/test/runtests_slayer_runner.jl b/test/runtests_slayer_runner.jl index 10d58bc87..3945de941 100644 --- a/test/runtests_slayer_runner.jl +++ b/test/runtests_slayer_runner.jl @@ -448,6 +448,16 @@ @test abs(r.Q_root[1] - Q2 * p2.tauk / p1.tauk) < 2e-2 # grid-resolution limited @test r.Q_root == ref.Q_root end + if mode === :coupled + # A Hermitian off-diagonal coupling moves the zeros off the diagonal entries' own; both classifications must + # still agree on the reported root, and every winding root counts −1. + dpc = ComplexF64[Δ1 0.3*sqrt(abs(Δ1 * Δ2)); 0.3*sqrt(abs(Δ1 * Δ2)) Δ2] + rc = run_slayer_from_inputs(model, [p1, p2], dpc, ctrl(:winding); rational_psi=[0.45, 0.72], rational_q=[2.0, 3.0]) + refc = run_slayer_from_inputs(model, [p1, p2], dpc, ctrl(:threshold)) + @test !isempty(rc.coupled_extraction.valid_roots) + @test all(w -> w ≈ -1, rc.coupled_extraction.valid_roots_winding) + @test rc.Q_root == refc.Q_root + end mktemp() do path, io close(io) h5open(path, "w") do f From b06477d0afec59db2fb01b583ea3605d64803752 Mon Sep 17 00:00:00 2001 From: d-burg Date: Wed, 28 Oct 2026 22:18:55 -0400 Subject: [PATCH 3/5] Tearing - TEST - Widen only the winding recheck circle, and expect the coupled root in develop's normalization Co-Authored-By: Claude Opus 5.5 --- .../Dispersion/GrowthRateExtraction.jl | 26 +++++++++---------- test/runtests_dispersion_winding.jl | 8 +++--- test/runtests_slayer_runner.jl | 6 ++--- 3 files changed, 18 insertions(+), 22 deletions(-) diff --git a/src/Tearing/Dispersion/GrowthRateExtraction.jl b/src/Tearing/Dispersion/GrowthRateExtraction.jl index 9d5da441a..db9210e31 100644 --- a/src/Tearing/Dispersion/GrowthRateExtraction.jl +++ b/src/Tearing/Dispersion/GrowthRateExtraction.jl @@ -175,11 +175,12 @@ rather than estimated. An intersection is a root when the count has the sign of `zero_winding`, a pole when it has the opposite sign, and is moved to `filtered_roots` (flag `:spurious`) when the count is zero, e.g. a contour near-miss or an unresolved pole-zero pair. A zero count is checked once more on -a circle eight times wider (still clear of every other intersection), because -the contour estimate of a pole can lie outside the small circle. A non-finite +a circle eight times wider, reaching up to 0.9 of the distance to the nearest +other intersection, because the contour estimate of a pole can lie outside the +small circle. A non-finite circle (the residual failed on it) is retried at half and a quarter of the -radius, and a zero count must agree on the half-radius circle; a circle that -stays unresolved is filtered too and flags `:winding_unresolved`. The count is +radius; a circle that stays unresolved is filtered too and flags +`:winding_unresolved`. The count is exact only for a residual that is single-valued and meromorphic on the circle: a branch cut crossing it, or a jump in the residual, can change the count. The classification replaces the pole threshold, the validity gate and the @@ -188,8 +189,8 @@ above-pole/outside-Re filter, so `pole_threshold`, `validity_rtol`, selection of `Q_root`, the `:geom`/`:gap` warnings and the secondary root are the same as in `:threshold` mode. With `polish = true` a root is polished inside its classifying circle, which is proven to contain it. `residual` is required. -Each intersection costs between about 130 and a few thousand residual -evaluations (two circles of 64–4096 samples). Only zeros that the scan resolves +Each intersection costs between 64 and a few thousand residual evaluations +(one or two circles of 64–4096 samples). Only zeros that the scan resolves as contour crossings are found, as in `:threshold` mode. # Spurious-root recursion @@ -969,9 +970,9 @@ end # Classify intersection `i` of `pts` by the winding number of `residual` on a circle around it: `(kind, W, r)` with kind # `:root`, `:pole`, `:none` (zero count) or `:unresolved`, the count W and the circle radius r. The circle is a fraction # `rfrac` of the local contour resolution and under half the distance to every other intersection; a non-finite circle is -# retried at r/2 and r/4. A zero count must agree on the half-radius circle (the contour estimate of a zero is accurate, so a -# disagreement means a second zero or a pole near the circle), and a zero count is checked once on a circle eight times -# wider, because the contour estimate of a pole can lie outside the small one. +# retried at r/2 and r/4, and a zero count is checked once on a circle eight times wider, because the contour estimate of a +# pole can lie outside the small one; that circle only has to exclude the other intersections, so it may reach 0.9 of the +# distance to the nearest. function _winding_classify(residual, pts::Vector{ComplexF64}, i::Int, re_paths, im_paths, h_grid::Float64; zero_winding::Int, rfrac::Float64) pt = pts[i] @@ -985,11 +986,8 @@ function _winding_classify(residual, pts::Vector{ComplexF64}, i::Int, re_paths, end isfinite(W) || return (:unresolved, NaN, r) n = round(Int, W) - if n * zero_winding >= 1 - W_half, _ = _winding_number(residual, pt, r / 2) - (isfinite(W_half) && round(Int, W_half) == n) || return (:unresolved, NaN, r) - elseif n == 0 - r_wide = _polish_trust_radius(pts, i, h_loc; k_cell=8rfrac, nn_frac=0.45) + if n == 0 + r_wide = _polish_trust_radius(pts, i, h_loc; k_cell=8rfrac, nn_frac=0.9) if r_wide > r W_wide, _ = _winding_number(residual, pt, r_wide) if isfinite(W_wide) && round(Int, W_wide) != 0 diff --git a/test/runtests_dispersion_winding.jl b/test/runtests_dispersion_winding.jl index b05ba6803..015fcc2cd 100644 --- a/test/runtests_dispersion_winding.jl +++ b/test/runtests_dispersion_winding.jl @@ -125,14 +125,12 @@ end @testset "a coupled two-surface determinant with an off-diagonal term" begin - # det of [d1 − a(Q) c; c d2 − b(Q)] with a, b functions of conj(Q) and real d1, d2, c is itself a function of - # conj(Q); its zeros solve (d1 − Q̄)(d2 − 2Q̄) = c² and are not those of either diagonal entry alone. + # The conjugate of det [d1 − Q c; c d2 − 2Q], a function of conj(Q) like the SLAYER coupled determinant; its zeros + # solve (d1 − Q)(d2 − 2Q) = c² and are not those of either diagonal entry alone. d1, d2, c = 0.31 + 0.12im, 0.18 - 0.05im, 0.07 det2(Q) = conj((d1 - ComplexF64(Q)) * (d2 - 2ComplexF64(Q)) - c^2) - # (d1 − x)(d2 − 2x) − c² = 0 for x = Q̄ A, B, C = 2.0, -(d2 + 2d1), d1 * d2 - c^2 - xs = [(-B + sqrt(B^2 - 4A * C)) / 2A, (-B - sqrt(B^2 - 4A * C)) / 2A] - zs = conj.(xs) + zs = [(-B + sqrt(B^2 - 4A * C)) / 2A, (-B - sqrt(B^2 - 4A * C)) / 2A] scan = brute_force_scan(det2, (-0.5, 0.8), (-0.6, 0.6); nre=53, nim=49, threaded=false) r = find_growth_rates(scan, 1.0; residual=det2, classification=:winding, zero_winding=-1) @test length(r.valid_roots) == 2 diff --git a/test/runtests_slayer_runner.jl b/test/runtests_slayer_runner.jl index 3945de941..4f7a4b8f3 100644 --- a/test/runtests_slayer_runner.jl +++ b/test/runtests_slayer_runner.jl @@ -441,11 +441,11 @@ @test abs(r.Q_root[2] - Q2) < 1e-3 @test maximum(abs.(r.Q_root .- ref.Q_root)) < 1e-3 else - # zeros at Q1 and Q2·τ_k2/τ_k1 in the reference normalization - # (surface k sees Q·τ_k1/τ_kk); the least stable is reported, + # zeros at Q1 and Q2·τ_k1/τ_k2 in the reference normalization + # (surface k sees Q·τ_kk/τ_k1); the least stable is reported, # unpolished since the coupled path does not polish @test length(r.coupled_extraction.valid_roots) == 2 - @test abs(r.Q_root[1] - Q2 * p2.tauk / p1.tauk) < 2e-2 # grid-resolution limited + @test abs(r.Q_root[1] - Q2 * p1.tauk / p2.tauk) < 2e-2 # grid-resolution limited @test r.Q_root == ref.Q_root end if mode === :coupled From 18d88c91c719c3c5913761e3ed5b546b1868f90c Mon Sep 17 00:00:00 2001 From: d-burg Date: Wed, 28 Oct 2026 22:19:15 -0400 Subject: [PATCH 4/5] Regression - TEST - Add a diiid_slayer_n1 variant with winding-number root classification Co-Authored-By: Claude Opus 5.5 --- regression-harness/cases/diiid_slayer_n1_winding.toml | 10 ++++++++++ 1 file changed, 10 insertions(+) create mode 100644 regression-harness/cases/diiid_slayer_n1_winding.toml diff --git a/regression-harness/cases/diiid_slayer_n1_winding.toml b/regression-harness/cases/diiid_slayer_n1_winding.toml new file mode 100644 index 000000000..b5c3a3963 --- /dev/null +++ b/regression-harness/cases/diiid_slayer_n1_winding.toml @@ -0,0 +1,10 @@ +# Regression case: the diiid_slayer_n1 deck with SLAYER roots classified by winding number instead of the +# pole threshold. Pins the winding path's roots and growth rates on a real DIII-D-like layer, and its +# runtime row tracks the cost of the classifying circles against the threshold case. +[case] +name = "diiid_slayer_n1_winding" +base = "diiid_slayer_n1" # Inherit example_dir, [quantities.*] and [overrides] from this case +description = "DIII-D-like H-mode equilibrium, n=1, SLAYER tearing-mode analysis with root_classification = winding" + +[overrides] +"SLAYER.root_classification" = "winding" # Classify contour intersections by the argument principle From 36d54d2fac446caa5fe5ccf24b1b8b5968afce69 Mon Sep 17 00:00:00 2001 From: Matthew Pharr Date: Thu, 29 Oct 2026 11:10:25 +0100 Subject: [PATCH 5/5] Tearing - DOCS - Clarify sampled winding-count limitations 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. --- .../Dispersion/GrowthRateExtraction.jl | 25 ++++++++----------- 1 file changed, 11 insertions(+), 14 deletions(-) diff --git a/src/Tearing/Dispersion/GrowthRateExtraction.jl b/src/Tearing/Dispersion/GrowthRateExtraction.jl index db9210e31..c21408965 100644 --- a/src/Tearing/Dispersion/GrowthRateExtraction.jl +++ b/src/Tearing/Dispersion/GrowthRateExtraction.jl @@ -170,8 +170,10 @@ instead of the pole-threshold heuristics. The counter-clockwise winding number of `residual` on a small circle around the intersection counts the zeros minus the poles inside it (for a holomorphic function; the sign is reversed for an anti-holomorphic one, hence `zero_winding`). The circle is resampled until no -phase step between neighboring samples exceeds π/2, so the count is exact -rather than estimated. An intersection is a root when the count has the sign of +sampled phase step exceeds π/2. This numerical criterion can miss rapid phase +turns between samples and does not certify the count. Additional midpoint and +count-convergence checks are possible defensive extensions. An intersection is +classified as a root when the estimated count has the sign of `zero_winding`, a pole when it has the opposite sign, and is moved to `filtered_roots` (flag `:spurious`) when the count is zero, e.g. a contour near-miss or an unresolved pole-zero pair. A zero count is checked once more on @@ -180,15 +182,15 @@ other intersection, because the contour estimate of a pole can lie outside the small circle. A non-finite circle (the residual failed on it) is retried at half and a quarter of the radius; a circle that stays unresolved is filtered too and flags -`:winding_unresolved`. The count is -exact only for a residual that is single-valued and meromorphic on the circle: +`:winding_unresolved`. The argument principle requires a residual that is +single-valued and meromorphic on and inside the circle: a branch cut crossing it, or a jump in the residual, can change the count. The classification replaces the pole threshold, the validity gate and the above-pole/outside-Re filter, so `pole_threshold`, `validity_rtol`, `filter_above_poles` and `filter_outside_re` decide nothing in this mode; the selection of `Q_root`, the `:geom`/`:gap` warnings and the secondary root are the same as in `:threshold` mode. With `polish = true` a root is polished inside -its classifying circle, which is proven to contain it. `residual` is required. +its classifying circle. `residual` is required. Each intersection costs between 64 and a few thousand residual evaluations (one or two circles of 64–4096 samples). Only zeros that the scan resolves as contour crossings are found, as in `:threshold` mode. @@ -698,7 +700,7 @@ function _run_analysis(re_paths::Vector{Vector{ComplexF64}}, # near-miss (e.g. a surface whose Δ' BVP failed → |Δ'| huge and # uncancellable by the inner layer) leaves |residual| stuck near the typical # scan magnitude. Only active when polishing AND a scale is supplied. - # A winding count already proves a zero, so the gate applies to the threshold classification only. + # Winding classification uses the sampled count instead of this residual-magnitude gate. do_gate = do_polish && validity_scale > 0 && !winding polish_evals = 0 @@ -821,7 +823,7 @@ function _run_analysis(re_paths::Vector{Vector{ComplexF64}}, for k in 1:length(sorted_pts) cand = sorted_pts[k] top_re = sorted_top[k] - # Legacy filter: above-pole + closed-loop outside-Re (a winding-classified root is a proven zero) + # Apply the above-pole and outside-Re heuristics only to threshold classification. legacy_reject = !winding && filter_above_poles && imag(cand) > max_pole_gamma && (!filter_outside_re || top_re) if legacy_reject @@ -892,8 +894,7 @@ end # can misjudge a root that sits close to a pole, and none of them is a proof. # The argument principle is: the winding number of f on a closed curve counts # zeros minus poles inside it (holomorphic f; the sign flips for a function of -# conj(Q)). Evaluated on a small circle around each intersection with the -# phase resolved, it classifies the intersection exactly. +# conj(Q)). Finite sampling can miss phase turns between points and is not a proof of the count. # --------------------------------------------------------------------- const _VALID_CLASSIFICATIONS = (:threshold, :winding) @@ -911,11 +912,7 @@ function _check_classification(classification::Symbol, residual, zero_winding::I return nothing end -# Counter-clockwise winding number of `f` on the circle |Q − Q0| = r. The -# circle is sampled at `nmin` points and the sampling doubled (reusing the -# previous samples) until no phase step between neighbors exceeds π/2, up to -# `nmax` points. Returns (W, max_step); W is NaN when `f` is non-finite or zero -# on the circle or the phase stays unresolved (a zero or pole on the circle). +# Estimate the counter-clockwise winding count by refining sampled phase steps; NaN means unresolved. function _winding_number(f, Q0::ComplexF64, r::Float64; nmin::Int=64, nmax::Int=4096) N = nmin vals = ComplexF64[f(Q0 + r * cis(2π * j / N)) for j in 0:N-1]