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

.gitignore

@@ -38,3 +38,6 @@ scratch/
 
 # Local profiling scratch (one-off study scripts, not part of the package)
 profiling/
+
+# Override-case rundirs materialized beside their example by the regression harness
+.regress-override-*

docs/development/hdf5-conventions.md

@@ -25,7 +25,7 @@ These rules govern `gpec.h5` (and any future GPEC-produced HDF5 output); harness
 - **"rational" over "singular"** in dataset names (`rational_psi`, `rational_q`, `rational_m`, `rational_n`, `rational_index`, `rational_count`) — kinetic/resistive runs are not singular at the rationals. Specifier order is standardized specifier-first (`rational_psi`, never `psi_rational`).
 - **Vector components** follow `[d]<variable>[_<representation>]_<coordinate>[dpsi]` — the variable always comes first and the coordinate is always the trailing subscript. A bare coordinate suffix is the **contravariant** component (`xi_psi` = ξ^ψ), `_cov_` marks the **covariant** one (`b_cov_theta` = b_θ), a leading `J` marks a **Jacobian-weighted** component (`Jxi_theta` = J·ξ^θ), and other representations sit in the same slot (`xi_clebsch_psi`, `dxi_clebsch_psidpsi`). Never drop the variable: `clebsch_psi` is wrong because it does not say *what* is being represented. There is no HDF5/netCDF standard for super- vs subscripts — flat `_` names are universal — so the typeset form always appears in the dataset's `long_name`.
 - **Coordinates**: the radial abscissa is `psi` (normalized poloidal flux ψ_N) and the poloidal one is `theta` in every group; never `psi_n`, `xs`, or `ys`.
-- **One name per physical quantity**: a quantity written in several groups carries the identical leaf name everywhere (`rational_psi` in `SingularSurfaces/`, `Solutions/GalerkinIntegration/`, and `SingularCoupling/`; `Delta_prime_matrix` in `SingularSurfaces/` and `Tearing/PerSurface/`; `dVdpsi` in `Profiles/`, `SingularSurfaces/`, and `KineticForces/<method>/`) — the group supplies the context, the leaf supplies the identity.
+- **One name per physical quantity**: a quantity written in several groups carries the identical leaf name everywhere (`rational_psi` in `SingularSurfaces/`, `Solutions/GalerkinIntegration/`, and `SingularCoupling/`; `Delta_prime_matrix` in `SingularSurfaces/` and, on the GGJ path, `Tearing/PerSurface/` — the r_s-referenced matrix the slab layer matches is a different quantity and so carries its own name, `Tearing/PerSurface/Delta_prime_matrix_rs`; `dVdpsi` in `Profiles/`, `SingularSurfaces/`, and `KineticForces/<method>/`) — the group supplies the context, the leaf supplies the identity.
 
 ## Inputs live only under `Input/`
 

docs/src/kinetic_forces.md

@@ -138,6 +138,35 @@ 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.) Logan, *Electromagnetic Torque in Tokamaks with Toroidal Asymmetries* (PhD thesis,
+Princeton, 2015), Ch. 7 Eq. 7.46 states the same resonance shifting/splitting result in GPEC's own
+notation.
+
+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.
+
+The check that this is right is that the NTV torque and 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 — agree to a fraction of a percent with the regularization off, and disagree
+by tens of percent with it on.
+
 ## 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_SLAYER_example/gpec.toml

@@ -61,8 +61,11 @@ dmlim                 = 0.2    # Truncate integration at (last_rational_q + dmli
 # Runs in the force_termination path (PE skipped). Kinetic profiles (incl. the
 # χ⊥/χ_φ heat/momentum diffusivities) are read from the standardized HDF5
 # kinetic file shared with the kinetic/NTV physics. Per-surface (uncoupled)
-# analysis: the headline is the unstable 2/1; the validity gate drops surfaces
-# with no real root (e.g. the 5/1, whose Δ' BVP yields a huge uncancellable Δ').
+# analysis: track the 2/1, 3/1, and 4/1 (the headline is the unstable 2/1).
+# The 5/1 and 6/1 sit near the edge (ψ_N ≥ 0.97), far from the core where the
+# slab-layer approximation holds; their Δ' are excessively large and likely
+# invalid, no Re/Im contour crossing polishes to a true zero, and the validity
+# gate correctly reports them as no_root rather than a spurious growth rate.
 enabled                 = true       # Run the SLAYER tearing-mode analysis
 inner_model             = "slayer_fitzpatrick" # Inner-layer Δ(Q) model: "slayer_fitzpatrick", "ggj_shooting", or "ggj_galerkin"
 scan_mode               = "amr"      # Q-plane scan strategy: "amr" (adaptive refinement) or "brute_force"

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 inside the near-axis orbit-width validity limit
 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 inside the near-axis orbit-width validity limit

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 inside the near-axis orbit-width validity limit

regression-harness/cases/diiid_n1.toml

@@ -411,3 +411,14 @@ extract = "value"
 label = "Runtime (s)"
 noise_threshold = 0.0
 order = 999
+
+# Root-area-weighted (coordinate-invariant) forcing field b̃ on the control surface: the
+# applied coil spectrum pushed through the √(J|∇ψ|) weight operator. Tracks the field-space
+# conversion that every ResponseMatrices/ field-space output depends on.
+[quantities.forcing_b_rootarea_norm]
+h5path = "PerturbedEquilibrium/forcing_b_root_area"
+type = "complex_vector"
+extract = "norm"
+label = "||forcing b~|| (root-area-weighted)"
+noise_threshold = 1e-10
+order = 1000

regression-harness/cases/diiid_slayer_n1.toml

@@ -1,6 +1,8 @@
 # Regression case: DIII-D-like H-mode n=1 SLAYER tearing-mode growth rates.
 # Uncoupled per-surface analysis with the AMR Q-plane scan; tracks the per-surface
-# layer parameters and the extracted roots (Q, ω, γ) for the inner three rationals.
+# layer parameters, the diamagnetic drive Q_e/Q_i, and the extracted roots (Q, ω, γ)
+# for the inner three rationals. The n-scaling of Q_e/Q_i, which no single-n case can
+# constrain, is covered by a unit test in test/runtests_slayer_inputs.jl instead.
 # Each [quantities.*] block names an HDF5 path in the run output, how to extract it,
 # and the noise floor below which a difference is treated as zero.
 [case]
@@ -89,6 +91,24 @@ label = "SLAYER iota_e"
 noise_threshold = 1e-12
 order = 19
 
+# Normalized diamagnetic frequencies. iota_e above is their ratio, so it is blind
+# to a common factor on both; these pin the drive itself.
+[quantities.slayer_Q_e]
+h5path = "Tearing/PerSurface/Q_e"
+type = "real_vector"
+extract = "all_real"
+label = "SLAYER Q_e"
+noise_threshold = 1e-12
+order = 20
+
+[quantities.slayer_Q_i]
+h5path = "Tearing/PerSurface/Q_i"
+type = "real_vector"
+extract = "all_real"
+label = "SLAYER Q_i"
+noise_threshold = 1e-12
+order = 21
+
 # Tearing eigenvalue (coupled mode → length 1). The headline deliverable:
 # a real, nonzero growth rate on a realistic equilibrium. Root-extraction
 # is sensitive to the AMR cell topology and ODE solver, so the Q/ω/γ

regression-harness/src/runner.jl

@@ -27,16 +27,18 @@ end
 Materialize the directory GPEC will actually run in.
 
 With no `overrides`, the example deck runs in place (returns it untouched). With overrides,
-the deck is copied to a throwaway temp dir and the named gpec.toml keys are patched there,
-so one shared example can serve several cases (e.g. a collisionless variant via
-`"KineticForces.nutype" => "zero"`). Override keys are dotted paths into the TOML; missing
-intermediate tables are created. Returns `(rundir, is_temp)`; the caller removes the temp
-tree when `is_temp`.
+the deck is copied to a uniquely named throwaway sibling of the example dir and the named gpec.toml keys are
+patched there, so one shared example can serve several cases (e.g. a collisionless variant via
+`"KineticForces.nutype" => "zero"`). The copy sits beside the original so relative file
+references in the deck (e.g. `SLAYER.profile_file = "../<other_example>/..."`) still resolve.
+Override keys are dotted paths into the TOML; missing intermediate tables are created.
+Returns `(rundir, is_temp)`; the caller removes `rundir` when `is_temp`.
 """
 function _materialize_rundir(example_path::String, overrides::Dict{String,Any})
     isempty(overrides) && return (example_path, false)
-    rundir = joinpath(mktempdir(), "deck")
-    cp(example_path, rundir)
+    # Unique per invocation, so concurrent runs from one checkout cannot share (and race on) a rundir.
+    rundir = mktempdir(dirname(example_path); prefix=".regress-override-" * basename(example_path) * "-", cleanup=false)
+    cp(example_path, rundir; force=true)
     rm(joinpath(rundir, "gpec.h5"); force=true)   # drop any stale output copied along
     cfg = TOML.parsefile(joinpath(rundir, "gpec.toml"))
     for (dotted, val) in overrides
@@ -523,7 +525,8 @@ function run_local(db::SQLite.DB, case_spec::CaseSpec, repo_root::String;
             rm(runinfo_file; force=true)
         end
         if rundir_is_temp && rundir !== nothing
-            rm(dirname(rundir); recursive=true, force=true)
+            # rundir sits beside the example; removing its parent would delete examples/ itself
+            rm(rundir; recursive=true, force=true)
         end
     end
 end
@@ -647,7 +650,8 @@ function run_at_commit(db::SQLite.DB, commit_hash::String, ref_name::String,
             rm(runinfo_file; force=true)
         end
         if rundir_is_temp && rundir !== nothing
-            rm(dirname(rundir); recursive=true, force=true)
+            # rundir sits beside the example; removing its parent would delete examples/ itself
+            rm(rundir; recursive=true, force=true)
         end
         if own_worktree && worktree_path !== nothing
             remove_worktree(worktree_path, repo_root)

src/Analysis/CoilForcing.jl

@@ -238,11 +238,7 @@ function plot_mode_spectrum(forcing_modes; mlow=nothing, mhigh=nothing, save_pat
 end
 
 # Helicity = sign(Bt)·sign(Ip); the SFL→machine toroidal-angle handedness.
-function _equil_helicity(equil)
-    bt = isnothing(equil.params.bt_sign) ? 1 : Int(sign(equil.params.bt_sign))
-    ip = isnothing(equil.params.crnt) ? 1 : Int(sign(equil.params.crnt))
-    return bt * ip
-end
+_equil_helicity(equil) = equil.params.bt_sign * equil.params.ip_sign
 
 # Unpack `modes` (a Vector of ForcingMode-like objects filtered to toroidal mode `n`,
 # or a `(m_vals, amplitudes)` tuple) into parallel m and complex-amplitude vectors.