PE/KF - BUGFIX! - Disable ideal-singularity regularization in self-consistent kinetic runs

reg_spot smooths the ideal 1/(m-nq) divergence of the displacements before they
drive the NTV integrand. The self-consistent kinetic Euler-Lagrange operator has
no such divergence -- det(F-bar) is complex and nonzero at the rationals (Park &
Logan, Phys. Plasmas 24, 032505 (2017) Eq. 70) -- so regularizing there suppresses
a finite physical response, and inconsistently: xi^psi is never regularized, so
damping the other two breaks their near-resonance cancellation in dB/B.

Measured (DIII-D-like, n=1) against the EL solution's own dissipation: NTV torque
0.1655 vs 0.1322 N*m with reg_spot=0.05, and 0.1324 vs 0.1322 (0.15%) with it off.
Kinetic runs now force reg_spot=0 and log the override; ideal runs are unchanged,
where without it the displacement and torque diverge by four orders of magnitude.

Co-Authored-By: Claude Fable 5 <noreply@anthropic.com>
Claude-Session: https://claude.ai/code/session_01LzbLFQKyuRE5DYZmLokKmk

Author: logan-nc

docs/src/kinetic_forces.md

@@ -89,6 +89,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

src/GeneralizedPerturbedEquilibrium.jl

@@ -603,8 +603,19 @@ function main_from_inputs(
             ft_ctrl = ForcingTerms.ForcingTermsControl()  # Use defaults
         end
 
+        pe_raw_in = inputs["PerturbedEquilibrium"]
+        # reg_spot smooths the ideal 1/(m−nq) divergence of ξ^ψ′ and ξ^α before they reach the NTV
+        # integrand. The self-consistent kinetic Euler-Lagrange operator has no such divergence —
+        # det(F̄) is complex and nonzero at the rationals — so regularizing there suppresses a finite
+        # physical response, and does so inconsistently (ξ^ψ is left unregularized, breaking the
+        # near-resonance cancellation in δB/B). Park & Logan, Phys. Plasmas 24, 032505 (2017) §III D.
+        if ctrl.kinetic_factor > 0 && get(pe_raw_in, "reg_spot", 0.0) != 0
+            @info "Self-consistent kinetic run: overriding reg_spot=$(pe_raw_in["reg_spot"]) with 0 " *
+                  "(the kinetic terms remove the ideal resonant singularity; see docs/src/kinetic_forces.md)"
+            pe_raw_in = merge(pe_raw_in, Dict("reg_spot" => 0.0))
+        end
         pe_ctrl = PerturbedEquilibrium.PerturbedEquilibriumControl(;
-            (Symbol(k) => v for (k, v) in inputs["PerturbedEquilibrium"])...
+            (Symbol(k) => v for (k, v) in pe_raw_in)...
         )
         pe_intr = PerturbedEquilibrium.PerturbedEquilibriumInternal(; dir_path=intr.dir_path)
 

src/PerturbedEquilibrium/FieldReconstruction.jl

@@ -19,6 +19,7 @@ where χ₁ = 2π·Ψ₀ [Park Phys. Plasmas 14, 052110 (2007) eq. 8-10].
 Clebsch displacement components for PENTRC (matches Fortran gpout_xclebsch):
     ξ^ψ         = xsp_mn    (unregularized)
     ∂ξ^ψ/∂ψ    = xmp1_mn   (regularized: xsp1 * singfac²/(singfac² + reg_spot²))
+                           reg_spot = 0 in kinetic runs — no ideal singularity to smooth
     ξ^α         = xms_mn    (regularized: -A⁻¹(B·xmp1 + C·xsp), divided by χ₁ in output)
 
 Contravariant displacement from Jacobian convolution (matches Fortran gpeq_contra):
@@ -404,7 +405,7 @@ function compute_clebsch_displacements(
         # xms = -(A\B)*xmp1 - (A\C)*xsp
         xsp_vec = view(xi_psi_modes, ipsi, :)
         mul!(xms_vec, bmat, xmp1_vec)                     # xms = B*xmp1
-        mul!(xms_vec, cmat_buf, xsp_vec, 1.0+0.0im, 1.0+0.0im)  # xms += C*xsp
+        mul!(xms_vec, cmat_buf, xsp_vec, 1.0 + 0.0im, 1.0 + 0.0im)  # xms += C*xsp
         # amat is positive-definite by construction (Newcomb kinetic-energy form), so cholesky is
         # safe. cholesky! factorizes in place (amat is a per-thread scratch buffer, refilled by
         # ffit.amats each surface), avoiding a fresh factorization allocation per surface.
@@ -980,10 +981,12 @@ function _apply_rzphi_transform(
     # Per-thread scratch (the immutable `ft` functor and `geom` are shared read-only): θ-space
     # transform inputs/outputs (length mtheta) and mode-space forward-DFT outputs (length mpert),
     # so the DFTs run in place with no per-surface allocation.
-    bufs = [(R=zeros(ComplexF64, mtheta), Z=zeros(ComplexF64, mtheta), P=zeros(ComplexF64, mtheta),
-             psi=zeros(ComplexF64, mtheta), th=zeros(ComplexF64, mtheta), ze=zeros(ComplexF64, mtheta),
-             Ro=zeros(ComplexF64, mpert), Zo=zeros(ComplexF64, mpert), Po=zeros(ComplexF64, mpert))
-            for _ in 1:Threads.maxthreadid()]
+    bufs = [
+        (R=zeros(ComplexF64, mtheta), Z=zeros(ComplexF64, mtheta), P=zeros(ComplexF64, mtheta),
+            psi=zeros(ComplexF64, mtheta), th=zeros(ComplexF64, mtheta), ze=zeros(ComplexF64, mtheta),
+            Ro=zeros(ComplexF64, mpert), Zo=zeros(ComplexF64, mpert), Po=zeros(ComplexF64, mpert))
+        for _ in 1:Threads.maxthreadid()
+    ]
 
     Threads.@threads :static for ipsi in 1:npsi
         buf = bufs[Threads.threadid()]

src/PerturbedEquilibrium/PerturbedEquilibriumStructs.jl

@@ -26,8 +26,9 @@ Medium Priority (defer for MWE):
   - `singular_point_method::String` - Method for singular point treatment (default: "standard")
 
 Regularization:
-    # High Priority (MWE)
-  - `reg_spot::Float64` - Regularization width for singular surface smoothing (default: 0.05). Set to 0 to disable. Must be ≥ 0.
+# High Priority (MWE)
+
+  - `reg_spot::Float64` - Regularization width for singular surface smoothing (default: 0.05). Set to 0 to disable. Must be ≥ 0. Forced to 0 in self-consistent kinetic runs, whose Euler-Lagrange operator has no resonant singularity to smooth.
 """
 @kwdef struct PerturbedEquilibriumControl
     # High Priority (MWE)
@@ -121,6 +122,7 @@ Metadata [n_rational] — identifies each (surface, n) row:
 
 Control-surface forcing/response spectra [numpert_total], in the three Pharr (2026) field
 representations (all tesla; no flux/weber is stored):
+
   - `forcing_b`/`response_b` - bare normal field b (Σ⁻¹·b̃)
   - `forcing_b_rootarea`/`response_b_rootarea` - root-area-weighted field b̃ (coordinate-invariant)
   - `forcing_b_area`/`response_b_area` - area-weighted field b̄ (= S·b̃; flux is Φ = A·b̄)
@@ -185,18 +187,18 @@ well-conditioned flux-space inductances L, Λ:
     rational_surface_idx::Vector{Int} = Int[]
 
     # Control-surface forcing/response spectra in the three weightings of field representations [numpert_total], tesla
-    forcing_b::Vector{ComplexF64}           = ComplexF64[]  # bare normal field b (forcing Φ_x)
-    forcing_b_rootarea::Vector{ComplexF64}  = ComplexF64[]  # root-area-weighted field b̃ (coordinate-invariant)
-    forcing_b_area::Vector{ComplexF64}      = ComplexF64[]  # area-weighted field b̄
-    response_b::Vector{ComplexF64}          = ComplexF64[]  # bare normal field b (response Φ_tot = P·Φ_x)
+    forcing_b::Vector{ComplexF64} = ComplexF64[]  # bare normal field b (forcing Φ_x)
+    forcing_b_rootarea::Vector{ComplexF64} = ComplexF64[]  # root-area-weighted field b̃ (coordinate-invariant)
+    forcing_b_area::Vector{ComplexF64} = ComplexF64[]  # area-weighted field b̄
+    response_b::Vector{ComplexF64} = ComplexF64[]  # bare normal field b (response Φ_tot = P·Φ_x)
     response_b_rootarea::Vector{ComplexF64} = ComplexF64[]  # root-area-weighted field b̃
-    response_b_area::Vector{ComplexF64}     = ComplexF64[]  # area-weighted field b̄
+    response_b_area::Vector{ComplexF64} = ComplexF64[]  # area-weighted field b̄
 
     # Control surface matrices [numpert_total × numpert_total], root-area-weighted field (b̃) space
-    plasma_inductance::Matrix{ComplexF64}  = zeros(ComplexF64, 0, 0)  # Λ̃ (field space)
+    plasma_inductance::Matrix{ComplexF64} = zeros(ComplexF64, 0, 0)  # Λ̃ (field space)
     surface_inductance::Matrix{ComplexF64} = zeros(ComplexF64, 0, 0)  # L̃ (field space)
-    permeability::Matrix{ComplexF64}       = zeros(ComplexF64, 0, 0)  # P̃ = R⁻¹·Λ·L⁻¹·R
-    reluctance::Matrix{ComplexF64}         = zeros(ComplexF64, 0, 0)  # ϱ̃ = R†·L⁻¹·(Λ−L)·L⁻¹·R
+    permeability::Matrix{ComplexF64} = zeros(ComplexF64, 0, 0)  # P̃ = R⁻¹·Λ·L⁻¹·R
+    reluctance::Matrix{ComplexF64} = zeros(ComplexF64, 0, 0)  # ϱ̃ = R†·L⁻¹·(Λ−L)·L⁻¹·R
     rootarea_to_area_weight::Matrix{ComplexF64} = zeros(ComplexF64, 0, 0)  # S = Σ/√A at psilim: b̃→b̄ recovery operator
     surface_area::Float64 = 0.0  # scalar control-surface area A = ∫J|∇ψ|dθ (flux: Φ = A·b̄; conform R = S·A)