Tearing - FEATURE - Classify tearing roots by winding number

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/` | `<method>/` (torque/energy profiles, `dTdpsi_ell` resolved by bounce harmonic on the `ell` scale, `KineticMatrices/`); multi-ion runs add `PerSpecies/<species>/<method>/` 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_<k>/`, `CriticalResonantField/` (+ `Scan/Surface_<k>/`) |
+| `Tearing/` | `PerSurface/` (+ `DpMatrix/`), `Roots/`, `LayerWidths/`, `Diagnostics/{ValidRoots,Poles,FilteredRoots}` (+ `*Winding` with `root_classification = winding`), `Scan/Surface_<k>/`, `CriticalResonantField/` (+ `Scan/Surface_<k>/`) |
 | `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`.

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

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)

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)

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{:fitzpatrick}) = -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, " *

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")