KineticForces - FEATURE! - Suppress kinetic terms inside the near-axis validity limit and report validity profiles

docs/src/kinetic_forces.md

@@ -138,6 +138,42 @@ terms respectively, and do not modify the stored kinetic profile splines.
        when computing the `toroidal_rotation_factor` back-solve. Julia uses a
        clean reimplementation with consistent pre-scaling derivatives throughout.
 
+## Regularization: why kinetic runs set `reg_spot = 0`
+
+`[PerturbedEquilibrium] reg_spot` smooths the displacements before they drive the NTV
+integrand, multiplying ``\xi^{\psi\prime}`` and ``\xi^\alpha`` by
+``Q^2/(Q^2 + \mathrm{reg\_spot}^2)`` with ``Q = m - nq``. It exists because **ideal** MHD is
+singular at the rationals: ``\bar F_\mathrm{ideal} = Q F Q`` has ``\det \bar F = 0`` there, so
+those two components diverge as ``1/Q`` and the torque integral does not converge.
+
+The **self-consistent kinetic** workflow has no such singularity. Park & Logan
+([Phys. Plasmas 24, 032505 (2017)](https://doi.org/10.1063/1.4978562), §III D) decompose the
+kinetic composite matrix as ``F_k = Q \bar F_k Q - P_l^\dagger Q - Q P_u + R_1`` where
+``R_1 \neq 0`` at ``Q = 0``; with finite torque ``\det \bar F`` is complex, its zeros leave the
+real ``\psi`` axis, and the singularity is removed from both the solution and the torque
+integral. (Torque-free kinetic energy principles instead *shift and split* the zeros to
+``\psi_r \mp r_{L,R}``, where the singularity is logarithmic and integrable — still not a case
+for smoothing.)
+
+GPEC therefore **forces `reg_spot = 0` whenever `kinetic_factor > 0`**, logging the override.
+Leaving it on suppresses a finite physical response and does so inconsistently — ``\xi^\psi``
+is never regularized, so damping the other two breaks their near-resonance cancellation in
+``\delta B/B`` and leaves a spurious residue driving the NTV integrand.
+
+Measured on the DIII-D-like H-mode case (n = 1, C-coil drive), comparing the NTV torque against
+the Euler–Lagrange solution's own dissipation ``-2n\,\mathrm{Im}\langle \xi, u_2\rangle/4\mu_0``
+— two independent calculations of the same quantity:
+
+| configuration | max ``|\xi^\alpha|`` | NTV torque [N·m] | EL dissipation [N·m] |
+|---|---|---|---|
+| ideal, `reg_spot = 0` | 643.8 | 6074.3 | — |
+| ideal, `reg_spot = 0.05` | 0.058 | 0.554 | — |
+| kinetic, `reg_spot = 0.05` | 0.055 | 0.1655 | 0.1322 |
+| kinetic, `reg_spot = 0` | 0.059 | **0.1324** | **0.1322** |
+
+The ideal rows show why the knob exists; the kinetic rows show why it must be off there — the
+two independent torques agree to 0.15 % with no regularization, and to 20 % with it.
+
 ## HDF5 outputs: complex torque convention and the EnergyIntegrals layout
 
 The method level of `KineticForces/<method>/` reports the two physical scalars a user

examples/DIIID-like_ideal_example/gpec.toml

@@ -87,5 +87,6 @@ f0type = "maxwellian"                   # Equilibrium distribution
 moment = "pressure"                     # Pressure-moment NTV torque
 atol_xlmda = 1e-9                       # Absolute tolerance for inner pitch + energy integrations
 rtol_xlmda = 1e-5                       # Relative tolerance for inner pitch + energy integrations
+axis_validity_suppression = true        # Suppress kinetic terms where the zero-orbit-width ordering fails near the axis (profile-derived boundary, no tuning parameters)
 write_outputs_to_HDF5 = true            # Write outputs to the HDF5 file
 verbose = true                          # Enable verbose logging

examples/Solovev_kinetic_NTV_example/gpec.toml

@@ -87,3 +87,4 @@ f0fac = 1      # Scale toroidal field at constant pressure (β, q change; Φ, p,
 
 [KineticForces]
 kinetic_file = "kinetic.dat"            # Kinetic profile file: psi_n, n_i, n_e, T_i, T_e, omega_E columns
+axis_validity_suppression = true        # Suppress kinetic terms where the zero-orbit-width ordering fails near the axis (profile-derived boundary, no tuning parameters)

examples/a10_kinetic_example/gpec.toml

@@ -59,3 +59,4 @@ nutype = "harmonic"            # Collision operator (zero, small, krook, harmoni
 f0type = "maxwellian"          # Distribution function (maxwellian, jkp, cgl)
 atol_xlmda = 1e-9              # Absolute tolerance for inner pitch + energy integrations
 rtol_xlmda = 1e-5              # Relative tolerance for inner pitch + energy integrations
+axis_validity_suppression = true        # Suppress kinetic terms where the zero-orbit-width ordering fails near the axis (profile-derived boundary, no tuning parameters)

src/Equilibrium/GridRefinement.jl

@@ -418,7 +418,7 @@ Build the refined pass-2 ψ grid from a formed pass-1 equilibrium: measured-curv
 density (`_knot_density`), equidistribution, a global minimum-spacing floor
 (`enforce_min_spacing`), and rational-surface bracketing (`bracket_mandatory_nodes`). `tau` is
 the target interpolation accuracy (`psi_accuracy`); `kin` optionally supplies kinetic profiles
-whose pedestal gradients attract knots; `mandatory` lists rational-surface ψ values to bracket;
+whose pedestal gradients attract knots; `mandatory` lists rational-surface ψ values to bracket; `pinned` lists ψ values inserted as plain knots without a cleared zone (kinetic-resonance surfaces);
 `singfac_min` and `n_min` (smallest |n| in the run) set each surface's matching half-stencil
 `dpsi = singfac_min/(n_min·|q′|)`, and the bracket half-width is `bracket_coef·dpsi` (floored at
 `min_spacing`). Rational surfaces are bracketed, not pinned: a knot on the surface would make the
@@ -428,6 +428,7 @@ function refined_psi_grid(equil::PlasmaEquilibrium;
     tau::Float64,
     kin::Union{Nothing,KineticProfileSplines}=nothing,
     mandatory::Vector{Float64}=Float64[],
+    pinned::Vector{Float64}=Float64[],
     singfac_min::Float64=1e-4,
     n_min::Int=1,
     bracket_coef::Float64=BRACKET_COEF,
@@ -449,6 +450,11 @@ function refined_psi_grid(equil::PlasmaEquilibrium;
     N == N_cap && M_total > N_cap &&
         @warn "refined_psi_grid: knot count capped at $N_cap (density integral wants $(ceil(Int, M_total))); psi_accuracy=$tau may not be attainable"
     grid = enforce_min_spacing(_equidistribute(xs, rho, N), min_spacing)
+    # Pinned knots (e.g. kinetic-resonance surfaces): knot-at-node semantics via
+    # merge_mandatory_nodes — no cleared zone. Inserted before rational bracketing, so a
+    # pinned node inside a rational's bracket zone is cleared by it (the Δ′ clean-interval
+    # requirement wins locally; the rational's own dense floor resolves that neighbourhood).
+    grid = isempty(pinned) ? grid : merge_mandatory_nodes(grid, pinned)
     isempty(mandatory) && return grid
     min_half_widths = [max(bracket_coef * singfac_min / (n_min * abs(equil.profiles.q_deriv(m))), min_spacing) for m in mandatory]
     return bracket_mandatory_nodes(grid, mandatory, min_half_widths, min_spacing)

src/ForceFreeStates/Fourfit.jl

@@ -443,6 +443,7 @@ rational's resolution window, preserving the Δ′-stencil structure the equilib
 (`Equilibrium.RATIONAL_RES_RADIUS`).
 """
 function core_capped_knots(xs::Vector{Float64}, rationals::Vector{Float64})::Vector{Int}
+    length(xs) < 3 && return collect(eachindex(xs))
     cap_edge = 0.1
     isempty(rationals) || (cap_edge = min(cap_edge, minimum(rationals) - Equilibrium.RATIONAL_RES_RADIUS))
     in_rational_window(x) = any(abs(x - r) <= Equilibrium.RATIONAL_RES_RADIUS for r in rationals)

src/ForceFreeStates/Kinetic.jl

@@ -1,3 +1,74 @@
+
+# Knots laid across the envelope's [ψ_c, 2ψ_c] transition. The envelope is a quintic smoothstep,
+# so a cubic spline needs several interior knots to follow it without overshoot; nine keeps the
+# residual below the kernel's own tolerance even on coarse decks.
+const BAND_KNOTS = 9
+
+"""
+    refine_grid_at_fbar_peaks(xs, kw, kt, evaluate, kin, equil, intr, psi_c;
+                              ngrid=1000, relaxed_frac=0.01, target=3, max_add=24) → (xs, kw, kt)
+
+Insert kinetic evaluation knots across near-singular structure of F̄ that the grid does not
+resolve. Scans cond(F̄) (the same operator `find_kinetic_singular_surfaces!` searches — Park &
+Logan Eq. 70, so shifted and split resonances are included), takes peaks between
+`relaxed_frac`·threshold and the singular threshold, measures each peak's FWHM, and adds knots
+only where fewer than `target` knots lie inside it. New points respect `MIN_KNOT_SPACING`, stay
+above the near-axis validity band, and are capped at `max_add`; each costs one kernel evaluation
+and the existing values are reused. A well-resolved grid inserts nothing.
+"""
+function refine_grid_at_fbar_peaks(xs::Vector{Float64}, kw::Array{ComplexF64,3}, kt::Array{ComplexF64,3},
+    evaluate::Function, kin::KineticMatrices, equil::Equilibrium.PlasmaEquilibrium,
+    intr::ForceFreeStatesInternal, psi_c::Float64;
+    ngrid::Int=1000, relaxed_frac::Float64=KINETIC_RELAXED_FRAC, target::Int=3, max_add::Int=24,
+    cond_threshold::Float64=KINETIC_SINGULAR_COND)
+
+    lo, hi = xs[1], xs[end]
+    scan = collect(range(lo, hi; length=ngrid))
+    hint = Ref(1)
+    cond_vals = [
+        try
+            evaluate_fbar_condition(x, kin, equil, intr; hint=hint)
+        catch
+            Inf
+        end for x in scan
+    ]
+
+    add = Float64[]
+    for i in 2:(ngrid-1)
+        c = cond_vals[i]
+        (c > cond_vals[i-1] && c > cond_vals[i+1] && relaxed_frac * cond_threshold < c <= cond_threshold) || continue
+        l = i
+        while l > 1 && cond_vals[l] > c / 2
+            l -= 1
+        end
+        r = i
+        while r < ngrid && cond_vals[r] > c / 2
+            r += 1
+        end
+        inside = count(x -> scan[l] <= x <= scan[r], xs)
+        inside >= target && continue
+        w = (scan[r] - scan[l]) / 3
+        for x in (scan[i], scan[i] - w, scan[i] + w)
+            (lo < x < hi && x > 2 * psi_c) || continue
+            any(y -> abs(y - x) < Equilibrium.MIN_KNOT_SPACING, xs) && continue
+            any(y -> abs(y - x) < Equilibrium.MIN_KNOT_SPACING, add) && continue
+            push!(add, x)
+        end
+    end
+    isempty(add) && return xs, kw, kt
+    length(add) > max_add && (add = sort(add)[1:max_add])
+
+    sort!(add)
+    @info "Kinetic grid: $(length(add)) knot(s) added across unresolved near-singular F̄ structure at " *
+          "ψ=$(round.(add; digits=4)) (cond peaks below the singular threshold)"
+    kw_new, kt_new = evaluate(add)
+    allxs = vcat(xs, add)
+    perm = sortperm(allxs)
+    kw_all = cat(kw, kw_new; dims=1)[perm, :, :]
+    kt_all = cat(kt, kt_new; dims=1)[perm, :, :]
+    return allxs[perm], kw_all, kt_all
+end
+
 """
     build_kinetic_matrix_splines(ctrl, equil, mats, intr, metric;
                         calculated_source=nothing)
@@ -27,11 +98,19 @@ function build_kinetic_matrix_splines(
     mats::MatrixSplines,
     intr::ForceFreeStatesInternal,
     metric::MetricData;
-    calculated_source::Union{Nothing,Function}=nothing
+    calculated_source::Union{Nothing,Function}=nothing,
+    axis_validity_psi_c::Float64=0.0
 )
     xs = metric.xs
     mpsi = length(xs)
 
+    # The near-axis validity envelope (KineticForces) has structure on the scale of the
+    # suppression boundary; coarse equilibrium grids cannot represent env·(increment), and the
+    # spline overshoot can land on a rational surface. Pin the band ends (the smoothstep is
+    # only C² there) and resolve the transition with a fixed set of knots.
+    band_knots(lo, hi) = axis_validity_psi_c > 0 ?
+                         [x for x in range(axis_validity_psi_c, 2 * axis_validity_psi_c; length=BAND_KNOTS) if lo < x < hi] : Float64[]
+
     # Get raw kinetic matrices (scaling is baked into each source)
     if ctrl.kinetic_source == "fixed"
         kw_flat, kt_flat = fixed_kinetic_matrices(intr.mpert, intr.numpert_total, mpsi, ctrl.kinetic_factor, intr.mlow, mats, xs)
@@ -42,15 +121,39 @@ function build_kinetic_matrix_splines(
             "calling build_kinetic_matrix_splines directly, or pass " *
             "`calculated_source=KineticForces.compute_calculated_kinetic_matrices` explicitly."
         )
-        kw_flat, kt_flat = calculated_source(ctrl, equil, intr, metric, mats)
+        band = band_knots(xs[1], xs[end])
+        if isempty(band)
+            kw_flat, kt_flat = calculated_source(ctrl, equil, intr, metric, mats)
+        else
+            xs = sort!(unique!(vcat(collect(xs), band)))
+            mpsi = length(xs)
+            kw_flat, kt_flat = calculated_source(ctrl, equil, intr, metric, mats; psis=xs)
+        end
         kw_flat .*= ctrl.kinetic_factor
         kt_flat .*= ctrl.kinetic_factor
     else
         error("Unknown kinetic_source: $(ctrl.kinetic_source). Must be \"fixed\" or \"calculated\"")
     end
 
     # Pre-compute FKG derived matrices (corresponds to Fortran method=0)
-    return _compute_fkg_matrices(mats, equil, intr, metric, kw_flat, kt_flat)
+    mats = _compute_fkg_matrices(mats, equil, intr, metric, kw_flat, kt_flat; xs=xs)
+
+    # The FKG splines now exist, so F̄ can be scanned: add knots only where near-singular structure
+    # (shifted/split kinetic resonances) falls in an interval that does not resolve it.
+    if ctrl.kinetic_source == "calculated" && calculated_source !== nothing && mats.kinetic !== nothing
+        xs2, kw_flat, kt_flat = refine_grid_at_fbar_peaks(
+            collect(xs), kw_flat, kt_flat,
+            psis -> begin
+                kwn, ktn = calculated_source(ctrl, equil, intr, metric, mats; psis=psis)
+                (kwn .* ctrl.kinetic_factor, ktn .* ctrl.kinetic_factor)
+            end,
+            mats.kinetic, equil, intr, axis_validity_psi_c)
+        if length(xs2) != length(xs)
+            xs = xs2
+            mats = _compute_fkg_matrices(mats, equil, intr, metric, kw_flat, kt_flat; xs=xs)
+        end
+    end
+    return mats
 end
 
 """
@@ -73,9 +176,9 @@ function _compute_fkg_matrices(
     intr::ForceFreeStatesInternal,
     metric::MetricData,
     kw_flat::Array{ComplexF64,3},
-    kt_flat::Array{ComplexF64,3}
+    kt_flat::Array{ComplexF64,3};
+    xs::Vector{Float64}=metric.xs
 )
-    xs = metric.xs
     mpsi = length(xs)
     np = intr.numpert_total
     mpert = intr.mpert

src/ForceFreeStates/Surfaces/Asymptotics.jl

@@ -74,7 +74,7 @@ function compute_sing_asymptotics(
 
     # This is the parameter α but for all modes - α = 0 for non-resonant modes
     power[ipert_res] .= -alpha
-    power[ipert_res .+ intr.numpert_total] .= alpha
+    power[ipert_res.+intr.numpert_total] .= alpha
 
     # Zeroth-order non-resonant solutions
     for ipert in 1:intr.numpert_total
@@ -111,7 +111,7 @@ function compute_sing_asymptotics(
         msg *= @sprintf("  m0mat(1,2)= %+.12e %+.12ei\n", real(m0mat[1, 2]), imag(m0mat[1, 2]))
         msg *= @sprintf("  m0mat(2,1)= %+.12e %+.12ei\n", real(m0mat[2, 1]), imag(m0mat[2, 1]))
         msg *= @sprintf("  m0mat(2,2)= %+.12e %+.12ei\n", real(m0mat[2, 2]), imag(m0mat[2, 2]))
-        di = m0mat[1, 1]*m0mat[2, 2] - m0mat[2, 1]*m0mat[1, 2]
+        di = m0mat[1, 1] * m0mat[2, 2] - m0mat[2, 1] * m0mat[1, 2]
         msg *= @sprintf("  di= %+.12e, alpha= %+.12e %+.12ei\n", real(di), real(alpha[1]), imag(alpha[1]))
         msg *= @sprintf("  psifac= %+.12e, r1=%d, ipert0=%d\n", singp.psifac, r1[1], ipert0)
         msg *= @sprintf("  vmat(ip,ip,2,0)= %+.8e %+.8ei\n", real(vmat[ipert0, ipert0, 2, 1]), imag(vmat[ipert0, ipert0, 2, 1]))
@@ -579,7 +579,7 @@ function sing_get_ua(sing_asymp::SingAsymptotics, dpsi::Float64)
 
     # Restore powers (unshear v→u) — matches Fortran STRIDE sing_get_ua
     for i in eachindex(r1)
-        pfac = pfac_base ^ sing_asymp.alpha[i]  # dpsi^α
+        pfac = pfac_base^sing_asymp.alpha[i]  # dpsi^α
         ua[:, r2[2*i-1], :] ./= pfac  # big solution column: /dpsi^α
         ua[:, r2[2*i], :] .*= pfac    # small solution column: *dpsi^α
         ua[r1[i], :, 1] ./= sqrtfac   # resonant row ξ: /√dpsi

src/ForceFreeStates/Surfaces/Finding.jl

@@ -125,9 +125,11 @@ function sing_lim!(intr::ForceFreeStatesInternal, ctrl::ForceFreeStatesControl,
     # strategy. Multi-n runs are not supported — the "outermost rational + dmlim/n" cutoff depends
     # on which n is used — and fall back to qhigh / psihigh truncation with a warning.
     if ctrl.set_psilim_via_dmlim && intr.nlow <= 0
-        error("sing_lim!: set_psilim_via_dmlim = true requires a resolved toroidal range, but got intr.nlow=$(intr.nlow). " *
-              "Assign intr.nlow / intr.nhigh (from ctrl.nn_low / ctrl.nn_high) before calling sing_lim!, " *
-              "or set set_psilim_via_dmlim = false to truncate via qhigh / psihigh instead.")
+        error(
+            "sing_lim!: set_psilim_via_dmlim = true requires a resolved toroidal range, but got intr.nlow=$(intr.nlow). " *
+            "Assign intr.nlow / intr.nhigh (from ctrl.nn_low / ctrl.nn_high) before calling sing_lim!, " *
+            "or set set_psilim_via_dmlim = false to truncate via qhigh / psihigh instead."
+        )
     elseif ctrl.set_psilim_via_dmlim && intr.nlow != intr.nhigh
         @warn "set_psilim_via_dmlim = true is ignored for multi-n runs (nn_low=$(intr.nlow), nn_high=$(intr.nhigh)); falling back to qhigh / psihigh truncation."
     elseif ctrl.set_psilim_via_dmlim
@@ -239,6 +241,11 @@ function evaluate_fbar_condition(psi::Float64, kin::KineticMatrices, equil::Equi
     return cond(fbar)
 end
 
+# Kinetic F̄ counts as singular above this condition number; structure within
+# KINETIC_RELAXED_FRAC of it is reported as near-singular (the shifted/split resonances).
+const KINETIC_SINGULAR_COND = 1.0e8
+const KINETIC_RELAXED_FRAC = 0.01
+
 """
     find_kinetic_singular_surfaces!(mats, equil, intr; ngrid=2000, cond_threshold=1e8)
 
@@ -257,7 +264,13 @@ Algorithm:
  3. Refine each peak with golden-section minimization of -cond
  4. Filter by threshold and resonance condition
 """
-function find_kinetic_singular_surfaces!(mats::MatrixSplines, equil::Equilibrium.PlasmaEquilibrium, intr::ForceFreeStatesInternal; ngrid::Int=2000, cond_threshold::Float64=1e8)
+function find_kinetic_singular_surfaces!(
+    mats::MatrixSplines,
+    equil::Equilibrium.PlasmaEquilibrium,
+    intr::ForceFreeStatesInternal;
+    ngrid::Int=2000,
+    cond_threshold::Float64=KINETIC_SINGULAR_COND
+)
     kin = mats.kinetic
     kin === nothing && error("find_kinetic_singular_surfaces! requires a kinetic fit; call build_kinetic_matrix_splines first")
     psilow = equil.profiles.xs[1]
@@ -282,11 +295,19 @@ function find_kinetic_singular_surfaces!(mats::MatrixSplines, equil::Equilibrium
     intr.kinsing_scan_threshold = cond_threshold
 
     # Find local maxima of cond(F̄): points where cond increases then decreases
-    peak_indices = Int[]
-    for i in 2:(ngrid-1)
-        if cond_vals[i] > cond_vals[i-1] && cond_vals[i] > cond_vals[i+1] && cond_vals[i] > cond_threshold
-            push!(peak_indices, i)
-        end
+    local_maxima = [i for i in 2:(ngrid-1) if cond_vals[i] > cond_vals[i-1] && cond_vals[i] > cond_vals[i+1]]
+    peak_indices = filter(i -> cond_vals[i] > cond_threshold, local_maxima)
+
+    # Peaks below the threshold are not singular surfaces, but they mark where the kinetic F̄ comes
+    # closest to singular — the shifted/split resonances of Park & Logan Eq. (70). Report the
+    # strongest few so sharp kinetic structure is visible rather than silent (on a DIII-D-like case
+    # these track the NTV torque-density peaks at low collisionality/rotation).
+    subthreshold = filter(i -> KINETIC_RELAXED_FRAC * cond_threshold < cond_vals[i] <= cond_threshold, local_maxima)
+    if !isempty(subthreshold)
+        top = sort(subthreshold; by=i -> cond_vals[i], rev=true)[1:min(3, length(subthreshold))]
+        @info "Kinetic F̄ near-singular structure below the singular threshold at " *
+              join(["ψ=$(round(psi_grid[i]; digits=4)) (cond=$(round(cond_vals[i]; sigdigits=3)))" for i in top], ", ") *
+              " — full scan in SingularSurfaces/Kinetic/scan_cond; check the ψ grid resolves these if results look grid-sensitive"
     end
 
     # Refine each peak to find the precise ψ location