KF - NEW FEATURE - Near-axis kinetic validity suppression with profile-derived boundary

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/ForceFreeStates/Kinetic.jl

@@ -26,11 +26,19 @@ function make_kinetic_matrix(
     ffit::FourFitVars,
     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=9) 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, mpsi, ctrl.kinetic_factor, intr.mlow, ffit, xs)
@@ -41,7 +49,14 @@ function make_kinetic_matrix(
             "calling make_kinetic_matrix directly, or pass " *
             "`calculated_source=KineticForces.compute_calculated_kinetic_matrices` explicitly."
         )
-        kw_flat, kt_flat = calculated_source(ctrl, equil, intr, metric, ffit)
+        band = band_knots(xs[1], xs[end])
+        if isempty(band)
+            kw_flat, kt_flat = calculated_source(ctrl, equil, intr, metric, ffit)
+        else
+            xs = sort!(unique!(vcat(collect(xs), band)))
+            mpsi = length(xs)
+            kw_flat, kt_flat = calculated_source(ctrl, equil, intr, metric, ffit; psis=xs)
+        end
         kw_flat .*= ctrl.kinetic_factor
         kt_flat .*= ctrl.kinetic_factor
     else
@@ -55,13 +70,13 @@ function make_kinetic_matrix(
     end
 
     # Pre-compute FKG derived matrices (corresponds to Fortran method=0)
-    _compute_fkg_matrices!(ffit, equil, intr, metric, kw_flat, kt_flat)
+    _compute_fkg_matrices!(ffit, equil, intr, metric, kw_flat, kt_flat; xs=xs)
 
     return nothing
 end
 
 """
-    _compute_fkg_matrices!(ffit, equil, intr, metric, kw_flat, kt_flat)
+    _compute_fkg_matrices!(ffit, equil, intr, metric, kw_flat, kt_flat; xs=xs)
 
 Pre-compute the derived F, K, G kinetic matrices at each ψ grid point and store as splines.
 This corresponds to `fourfit_kinetic_matrix` method=0 in the Fortran code (Fortran `fourfit.F` lines 1170-1260).
@@ -78,9 +93,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/GeneralizedPerturbedEquilibrium.jl

@@ -76,7 +76,7 @@ using .ForceFreeStates: eulerlagrange_integration, free_run, normalize_eigenfunc
 using .ForceFreeStates: galerkin_solve, write_galerkin!, GalerkinResult, gal_matched_odestate
 
 const _DEPRECATED_FFS_KEYS = ("mer_flag", "force_wv_symmetry", "ode_flag", "cyl_flag", "mat_flag", "reform_eq_with_psilim",
-                              "use_riccati", "use_parallel", "parallel_threads", "populate_dense_xi")
+    "use_riccati", "use_parallel", "parallel_threads", "populate_dense_xi")
 const _DEPRECATED_EQUIL_KEYS = ("power_bp", "power_b", "power_r", "power_rc")
 
 # Drop deprecated keys from a parsed gpec.toml section so legacy files keep parsing
@@ -419,12 +419,21 @@ function main_from_inputs(
         # Inject the KineticForces callback so the "calculated" source can
         # invoke compute_calculated_kinetic_matrices without ForceFreeStates
         # importing KineticForces (which would invert the load order).
-        calculated_cb = (c, e, i, m, f) ->
-            KineticForces.compute_calculated_kinetic_matrices(
-                c, e, i, m, f;
-                kf_ctrl=kf_ctrl, kinetic_profiles=kinetic_profiles)
+        calculated_cb =
+            (c, e, i, m, f; psis=nothing) ->
+                KineticForces.compute_calculated_kinetic_matrices(
+                    c, e, i, m, f;
+                    kf_ctrl=kf_ctrl, kinetic_profiles=kinetic_profiles, psis=psis)
+        # Near-axis validity boundary for the calculated kinetic matrices: the envelope band
+        # must be resolved by the kernel grid, so make_kinetic_matrix needs its location.
+        axis_psi_c = 0.0
+        if ctrl.kinetic_source == "calculated" && kf_ctrl.axis_validity_suppression && kinetic_profiles !== nothing
+            axis_psi_c = KineticForces.kinetic_axis_validity_psi(
+                kinetic_profiles, equil;
+                zi=kf_ctrl.zi, mi=kf_ctrl.mi, electron=kf_ctrl.electron)
+        end
         make_kinetic_matrix(ctrl, equil, ffit, intr, metric;
-            calculated_source=calculated_cb)
+            calculated_source=calculated_cb, axis_validity_psi_c=axis_psi_c)
 
         # Find kinetically-displaced singular surfaces (zeros of det(F̄)) for ODE crossings.
         # Matches Fortran ksing_find (sing.f:1486-1616). singfac_min > 0 gates crossings;
@@ -495,6 +504,8 @@ function main_from_inputs(
             inputs,
             forcing_modes_snapshot,
             gal_data;
+            kinetic_profiles=kinetic_profiles,
+            kf_ctrl=kf_ctrl,
             locstab=locstab,
             ballooning_boundary=ballooning_boundary
         )
@@ -705,7 +716,9 @@ function write_outputs_to_HDF5(
     forcing_modes::Union{Nothing,Vector{ForcingTerms.ForcingMode}}=nothing,
     gal_data::Union{GalerkinResult,Nothing}=nothing;
     locstab::Union{FastInterpolations.CubicSeriesInterpolant,Nothing}=nothing,
-    ballooning_boundary=(psi=Float64[], alpha=Float64[], alpha_critical=Float64[])
+    ballooning_boundary=(psi=Float64[], alpha=Float64[], alpha_critical=Float64[]),
+    kinetic_profiles=nothing,
+    kf_ctrl=nothing
 )
 
     # Idempotent: already done if a PerturbedEquilibrium stage ran. Leaves the stores empty
@@ -821,6 +834,28 @@ function write_outputs_to_HDF5(
         out_h5["$fwd/xi_s"] = odet.xi_s_store
         out_h5["$fwd/crit"] = odet.crit_store
 
+        # Kinetic-model validity diagnostics: orbit-width scales vs local geometry and profile
+        # gradient lengths, whenever kinetic profiles were used (self-consistent matrices or NTV
+        # post-processing). Diagnostic only — nothing outside the near-axis envelope is suppressed.
+        if kinetic_profiles !== nothing
+            kfc = kf_ctrl === nothing ? KineticForces.KineticForcesControl() : kf_ctrl
+            vp = KineticForces.kinetic_validity_profiles(kinetic_profiles, equil;
+                zi=kfc.zi, mi=kfc.mi, electron=kfc.electron)
+            vg = "KineticForces/Validity"
+            out_h5["$vg/psi"] = vp.psi
+            out_h5["$vg/rho_i"] = vp.rho_i
+            out_h5["$vg/rho_banana"] = vp.rho_banana
+            out_h5["$vg/rho_theta"] = vp.rho_theta
+            out_h5["$vg/w_potato"] = vp.w_potato
+            out_h5["$vg/r_minor"] = vp.r_minor
+            out_h5["$vg/L_p"] = vp.L_p
+            out_h5["$vg/L_q"] = vp.L_q
+            out_h5["$vg/d_separatrix"] = vp.d_separatrix
+            out_h5["$vg/psi_c"] = vp.psi_c
+            out_h5["$vg/envelope"] = kfc.axis_validity_suppression ? vp.envelope : ones(length(vp.psi))
+            out_h5["$vg/is_valid"] = Int8.(vp.is_valid)
+        end
+
         # Write edge stability scan data (only present when psiedge < psilim).
         # Generalized (W, N) pencil energies — power-normalized, Jacobian-invariant; these are
         # the values findmax_dW_edge! uses to choose the truncation point.

src/HDF5Schema.jl

@@ -141,6 +141,25 @@ const MAIN_H5_ANNOTATIONS = [
     "LocalStability/alpha_critical" =>
         (; long_name="critical normalized pressure gradient α for first ballooning stability", dims=("psi_ballooning",),
             attach=(1 => "LocalStability/ballooning_psi",)),
+    # --- KineticForces/Validity/ ---
+    "KineticForces/Validity/psi" => (; long_name="normalized poloidal flux ψ_N of the kinetic validity profiles", scale="psi"),
+    "KineticForces/Validity/rho_i" => (; long_name="thermal ion gyroradius √(2mT)/(Z·e·B₀)", units="m", attach=(1 => "KineticForces/Validity/psi",)),
+    "KineticForces/Validity/rho_banana" => (; long_name="thermal banana orbit width q·ρ_i/√ε", units="m", attach=(1 => "KineticForces/Validity/psi",)),
+    "KineticForces/Validity/rho_theta" => (; long_name="thermal poloidal gyroradius q·ρ_i/ε", units="m", attach=(1 => "KineticForces/Validity/psi",)),
+    "KineticForces/Validity/w_potato" => (; long_name="potato orbit width (q²ρ_i²R₀)^(1/3)", units="m", attach=(1 => "KineticForces/Validity/psi",)),
+    "KineticForces/Validity/r_minor" => (; long_name="surface-average minor radius ⟨r⟩", units="m", attach=(1 => "KineticForces/Validity/psi",)),
+    "KineticForces/Validity/L_p" => (; long_name="pressure gradient scale length |p|/|dp/dr|", units="m", attach=(1 => "KineticForces/Validity/psi",)),
+    "KineticForces/Validity/L_q" => (; long_name="safety-factor gradient scale length |q|/|dq/dr|", units="m", attach=(1 => "KineticForces/Validity/psi",)),
+    "KineticForces/Validity/d_separatrix" => (; long_name="distance to the separatrix ⟨r⟩(1) − ⟨r⟩(ψ)", units="m", attach=(1 => "KineticForces/Validity/psi",)),
+    "KineticForces/Validity/psi_c" => (; long_name="near-axis kinetic validity boundary: outermost ψ_N where a thermal orbit width reaches ⟨r⟩", units="1"),
+    "KineticForces/Validity/envelope" =>
+        (; long_name="near-axis suppression envelope applied to the calculated kinetic terms (1 = unsuppressed)", units="1", attach=(1 => "KineticForces/Validity/psi",)),
+    "KineticForces/Validity/is_valid" => (;
+        long_name="1 where every zero-orbit-width ordering holds: max orbit width < ⟨r⟩, ρ_banana < L_p and L_q, max orbit width < d_separatrix",
+        units="1",
+        attach=(1 => "KineticForces/Validity/psi",)
+    ),
+
     # --- ForceFreeStates/Solutions/ForwardIntegration/ ---
     "ForceFreeStates/Solutions/ForwardIntegration/nstep" => (; long_name="number of saved solution snapshots"),
     "ForceFreeStates/Solutions/ForwardIntegration/nstep_total" => (; long_name="total ODE solver steps taken"),
@@ -239,7 +258,7 @@ const MAIN_H5_ANNOTATIONS = [
     "SurfaceGeometries/Plasma/z" => (; long_name="Cartesian z of plasma-surface point cloud", units="m"),
     "SurfaceGeometries/Wall/x" => (; long_name="Cartesian x of wall point cloud", units="m"),
     "SurfaceGeometries/Wall/y" => (; long_name="Cartesian y of wall point cloud", units="m"),
-    "SurfaceGeometries/Wall/z" => (; long_name="Cartesian z of wall point cloud", units="m"),
+    "SurfaceGeometries/Wall/z" => (; long_name="Cartesian z of wall point cloud", units="m")
 ]
 
 # Euler-Lagrange operator matrices: same wording per letter, Ideal/ and Kinetic/ variants.
@@ -252,7 +271,7 @@ const _ELM_IDEAL_LETTERS = [
     ("H", "Euler-Lagrange primitive coefficient matrix H"),
     ("F", "Euler-Lagrange derived coefficient matrix F"),
     ("K", "Euler-Lagrange derived coefficient matrix K"),
-    ("G", "Euler-Lagrange derived coefficient matrix G"),
+    ("G", "Euler-Lagrange derived coefficient matrix G")
 ]
 # The kinetic branch overwrites only A, B, C, K, G and adds f0; D, E, H, F are shared
 # unchanged from the ideal set and are not re-emitted.
@@ -262,14 +281,19 @@ const _ELM_KINETIC_LETTERS = [
     ("C", "Euler-Lagrange primitive coefficient matrix C"),
     ("K", "Euler-Lagrange derived coefficient matrix K"),
     ("G", "Euler-Lagrange derived coefficient matrix G"),
-    ("f0", "raw kinetic component matrix f0"),
+    ("f0", "raw kinetic component matrix f0")
 ]
 const ELM_H5_ANNOTATIONS = vcat(
     ["ForceFreeStates/EulerLagrangeMatrices/psi" => (; long_name="normalized poloidal flux ψ_N grid of the operator matrices", scale="psi")],
-    ["ForceFreeStates/EulerLagrangeMatrices/Ideal/$l" =>
-        (; long_name="ideal " * d, dims=("psi", "mode_row", "mode_col"), attach=(1 => "ForceFreeStates/EulerLagrangeMatrices/psi",)) for (l, d) in _ELM_IDEAL_LETTERS],
-    ["ForceFreeStates/EulerLagrangeMatrices/Kinetic/$l" =>
-        (; long_name="kinetic-modified " * d, dims=("psi", "mode_row", "mode_col"), attach=(1 => "ForceFreeStates/EulerLagrangeMatrices/psi",)) for (l, d) in _ELM_KINETIC_LETTERS]
+    [
+        "ForceFreeStates/EulerLagrangeMatrices/Ideal/$l" =>
+            (; long_name="ideal " * d, dims=("psi", "mode_row", "mode_col"), attach=(1 => "ForceFreeStates/EulerLagrangeMatrices/psi",)) for (l, d) in _ELM_IDEAL_LETTERS
+    ],
+    [
+        "ForceFreeStates/EulerLagrangeMatrices/Kinetic/$l" =>
+            (; long_name="kinetic-modified " * d, dims=("psi", "mode_row", "mode_col"), attach=(1 => "ForceFreeStates/EulerLagrangeMatrices/psi",)) for
+        (l, d) in _ELM_KINETIC_LETTERS
+    ]
 )
 
 """

src/KineticForces/CalculatedKineticMatrices.jl

@@ -31,33 +31,39 @@ tracked as follow-up work blocked on PR #196 — see the plan's "Out of scope"
 section.
 
 # Arguments
-- `ffs_ctrl`: ForceFreeStatesControl (carries `kinetic_factor`, `kinetic_source`)
-- `equil`: PlasmaEquilibrium with 2D interpolants and named profile/geometry splines
-- `ffs_intr`: ForceFreeStatesInternal (mode indexing)
-- `metric`: MetricData (provides ψ grid via `metric.xs`)
-- `ffit`: FourFitVars (used only for `numpert_total` cross-check)
+
+  - `ffs_ctrl`: ForceFreeStatesControl (carries `kinetic_factor`, `kinetic_source`)
+  - `equil`: PlasmaEquilibrium with 2D interpolants and named profile/geometry splines
+  - `ffs_intr`: ForceFreeStatesInternal (mode indexing)
+  - `metric`: MetricData (provides ψ grid via `metric.xs`)
+  - `ffit`: FourFitVars (used only for `numpert_total` cross-check)
 
 # Keyword arguments
-- `kf_ctrl`: KineticForcesControl, defaults to `KineticForcesControl()`. Used to
-  carry NTV-specific knobs (nl, zi, mi, wdfac, divxfac, electron) that the
-  KineticForces kernel needs but ForceFreeStatesControl does not expose.
-- `kinetic_profiles::Equilibrium.KineticProfileSplines`: Required. Named kinetic-
-  profile splines loaded via `Equilibrium.load_kinetic_profiles`.
+
+  - `kf_ctrl`: KineticForcesControl, defaults to `KineticForcesControl()`. Used to
+    carry NTV-specific knobs (nl, zi, mi, wdfac, divxfac, electron) that the
+    KineticForces kernel needs but ForceFreeStatesControl does not expose.
+  - `kinetic_profiles::Equilibrium.KineticProfileSplines`: Required. Named kinetic-
+    profile splines loaded via `Equilibrium.load_kinetic_profiles`.
 
 # Returns
-- `kw_flat::Array{ComplexF64,3}`: Energy matrices, shape `(mpsi, np^2, 6)`
-- `kt_flat::Array{ComplexF64,3}`: Torque matrices, shape `(mpsi, np^2, 6)`
+
+  - `kw_flat::Array{ComplexF64,3}`: Energy matrices, shape `(mpsi, np^2, 6)`
+  - `kt_flat::Array{ComplexF64,3}`: Torque matrices, shape `(mpsi, np^2, 6)`
 """
 function compute_calculated_kinetic_matrices(
     _ffs_ctrl,
     equil,
     ffs_intr,
     metric,
     ffit;
-    kf_ctrl::KineticForcesControl = KineticForcesControl(),
+    kf_ctrl::KineticForcesControl=KineticForcesControl(),
     kinetic_profiles::Equilibrium.KineticProfileSplines,
+    psis::Union{Nothing,Vector{Float64}}=nothing
 )
-    xs = metric.xs
+    # The kernel is a pure function of psi (it evaluates equilibrium splines), so it can be
+    # driven over any knot list; default is the full equilibrium grid.
+    xs = psis === nothing ? metric.xs : psis
     mpsi = length(xs)
     mpert = ffs_intr.mpert
     npert = ffs_intr.npert
@@ -93,22 +99,36 @@ function compute_calculated_kinetic_matrices(
     # ipsi row of kw_flat/kt_flat. Per-thread copies of kf_intr provide isolated
     # tpsi_* θ-grid buffers and interpolant hint refs; geometric/profile splines
     # are read-only and safely shared through deepcopy semantics.
+    # Near-axis validity envelope: suppress the drift-kinetic increments where the
+    # zero-orbit-width ordering fails; kernel evaluation is skipped where it is identically 0.
+    env = ones(Float64, mpsi)
+    if kf_ctrl.axis_validity_suppression
+        psi_c = kinetic_axis_validity_psi(kinetic_profiles, equil;
+            zi=kf_ctrl.zi, mi=kf_ctrl.mi, electron=kf_ctrl.electron)
+        if psi_c > 0
+            env .= kinetic_axis_validity_envelope.(xs, psi_c)
+            @info "Kinetic axis-validity suppression: psi_c=$(round(psi_c; sigdigits=3)), envelope reaches 1 at " *
+                  "psi=$(round(2 * psi_c; sigdigits=3)) (kernel evaluation skipped below psi_c)" maxlog = 1
+        end
+    end
+
     nl = kf_ctrl.nl
     nthreads = Threads.maxthreadid()
-    thread_intrs   = [deepcopy(kf_intr) for _ in 1:nthreads]
-    thread_full_w  = [zeros(ComplexF64, mpert, mpert, 6) for _ in 1:nthreads]
-    thread_full_t  = [zeros(ComplexF64, mpert, mpert, 6) for _ in 1:nthreads]
+    thread_intrs = [deepcopy(kf_intr) for _ in 1:nthreads]
+    thread_full_w = [zeros(ComplexF64, mpert, mpert, 6) for _ in 1:nthreads]
+    thread_full_t = [zeros(ComplexF64, mpert, mpert, 6) for _ in 1:nthreads]
     thread_block_w = [zeros(ComplexF64, mpert, mpert, 6) for _ in 1:nthreads]
     thread_block_t = [zeros(ComplexF64, mpert, mpert, 6) for _ in 1:nthreads]
 
     Threads.@threads for ipsi in 1:mpsi
-        tid     = Threads.threadid()
-        intr_t  = thread_intrs[tid]
-        full_w  = thread_full_w[tid]
-        full_t  = thread_full_t[tid]
+        tid = Threads.threadid()
+        intr_t = thread_intrs[tid]
+        full_w = thread_full_w[tid]
+        full_t = thread_full_t[tid]
         block_w = thread_block_w[tid]
         block_t = thread_block_t[tid]
-        psi     = xs[ipsi]
+        psi = xs[ipsi]
+        env[ipsi] == 0.0 && continue
         for in_idx in 1:npert
             n = ffs_intr.nlow + in_idx - 1
             fill!(full_w, 0)
@@ -131,8 +151,8 @@ function compute_calculated_kinetic_matrices(
             row_offset = (in_idx - 1) * mpert
             for k in 1:6, j in 1:mpert, i in 1:mpert
                 idx = (row_offset + j - 1) * np + (row_offset + i)
-                kw_flat[ipsi, idx, k] = full_w[i, j, k]
-                kt_flat[ipsi, idx, k] = full_t[i, j, k]
+                kw_flat[ipsi, idx, k] = env[ipsi] * full_w[i, j, k]
+                kt_flat[ipsi, idx, k] = env[ipsi] * full_t[i, j, k]
             end
         end
     end