From 09cb856d3eb07464e01f24b505ef4297bc85738b Mon Sep 17 00:00:00 2001 From: logan-nc Date: Thu, 20 Aug 2026 11:47:01 -0400 Subject: [PATCH 1/4] KF - NEW FEATURE - Suppress kinetic terms where the zero-orbit-width ordering fails near the axis Physics ruling (issue #376 DIII-D kinetic pathology): the drift-kinetic model loses validity where thermal ion orbit widths reach the local minor radius. psi_c = outermost crossing of by max(potato width (q^2 rho^2 R0)^(1/3), banana width q rho/sqrt(eps), poloidal gyroradius q rho/eps), computed from the equilibrium and kinetic profiles at runtime -- no user tuning parameters (the Fortran ktanh_flag precedent needed four). A C2 quintic envelope zeroes the calculated kinetic increments below psi_c (kernel evaluation skipped) and rises to 1 at 2 psi_c; the same boundary and envelope apply to the NTV torque psi quadrature (one source of truth). One Bool (axis_validity_suppression, default true) to disable for debugging. Co-Authored-By: Claude Fable 5 Claude-Session: https://claude.ai/code/session_01LzbLFQKyuRE5DYZmLokKmk --- examples/DIIID-like_ideal_example/gpec.toml | 1 + .../Solovev_kinetic_NTV_example/gpec.toml | 1 + examples/a10_kinetic_example/gpec.toml | 1 + .../CalculatedKineticMatrices.jl | 68 +++++++++----- src/KineticForces/Compute.jl | 61 ++++++++----- src/KineticForces/KineticForcesStructs.jl | 88 ++++++++++--------- src/KineticForces/Utils.jl | 56 +++++++++++- 7 files changed, 188 insertions(+), 88 deletions(-) diff --git a/examples/DIIID-like_ideal_example/gpec.toml b/examples/DIIID-like_ideal_example/gpec.toml index b84e12872..1edc25231 100644 --- a/examples/DIIID-like_ideal_example/gpec.toml +++ b/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 diff --git a/examples/Solovev_kinetic_NTV_example/gpec.toml b/examples/Solovev_kinetic_NTV_example/gpec.toml index e52b0d537..7a54a0397 100644 --- a/examples/Solovev_kinetic_NTV_example/gpec.toml +++ b/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) diff --git a/examples/a10_kinetic_example/gpec.toml b/examples/a10_kinetic_example/gpec.toml index fd8881122..ae602799f 100644 --- a/examples/a10_kinetic_example/gpec.toml +++ b/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) diff --git a/src/KineticForces/CalculatedKineticMatrices.jl b/src/KineticForces/CalculatedKineticMatrices.jl index d36fead99..9d8a19def 100644 --- a/src/KineticForces/CalculatedKineticMatrices.jl +++ b/src/KineticForces/CalculatedKineticMatrices.jl @@ -31,22 +31,25 @@ 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, @@ -54,10 +57,13 @@ function compute_calculated_kinetic_matrices( 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)); $(count(iszero, env)) of $mpsi surfaces skipped" + 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 diff --git a/src/KineticForces/Compute.jl b/src/KineticForces/Compute.jl index 4175a62b4..3d9c3d8cf 100644 --- a/src/KineticForces/Compute.jl +++ b/src/KineticForces/Compute.jl @@ -41,7 +41,7 @@ Warn when the ψ torque quadrature terminated without satisfying the requested t silent-garbage scenario for weak applied fields, since NTV scales as δB². """ function check_psi_quadrature_convergence(total::ComplexF64, quad_err::Float64, - ctrl::KineticForcesControl, method::String) + ctrl::KineticForcesControl, method::String) if quad_err > max(ctrl.atol_psi, ctrl.rtol_psi * abs(total)) @warn "ψ torque quadrature ($method) did not converge within maxevals_psi=$(ctrl.maxevals_psi): " * "error estimate $quad_err vs |T|=$(abs(total)) N·m. Raise maxevals_psi or loosen rtol_psi." @@ -63,14 +63,16 @@ diagnostic T(ψ) profile at no extra cost (the values are computed anyway — we just keep them). # Returns + NamedTuple with: -- `total::ComplexF64`: Total integrated torque -- `torque_profile`: NamedTuple of (psi, dtdpsi, t_cumulative) from evaluation points -- `matrix_integrated`: Trapezoidal-integrated mpert×mpert×6 matrix (if matrix method) -- `psi_nsteps::Int`: Number of integrand evaluations -- `psi_quad_error::Float64`: Quadrature error estimate for the total torque -- `panel_psis::Vector{Float64}`: Quadrature panel boundaries actually used (bounds + interior resonant surfaces) -- `resonance_psis::Vector{Float64}`: Located kinetic-resonance ψ surfaces (Ω_ℓ(x=1)=0), for diagnostics/plotting + + - `total::ComplexF64`: Total integrated torque + - `torque_profile`: NamedTuple of (psi, dtdpsi, t_cumulative) from evaluation points + - `matrix_integrated`: Trapezoidal-integrated mpert×mpert×6 matrix (if matrix method) + - `psi_nsteps::Int`: Number of integrand evaluations + - `psi_quad_error::Float64`: Quadrature error estimate for the total torque + - `panel_psis::Vector{Float64}`: Quadrature panel boundaries actually used (bounds + interior resonant surfaces) + - `resonance_psis::Vector{Float64}`: Located kinetic-resonance ψ surfaces (Ω_ℓ(x=1)=0), for diagnostics/plotting The integral is paneled at the rational-surface ψ locations (`intr.sing_psis`) and capped at `ctrl.maxevals_psi` evaluations; a warning is emitted if the quadrature fails to reach @@ -95,9 +97,21 @@ function integrate_psi_quadgk( if x0 >= xout return (total=ComplexF64(0.0), torque_profile=nothing, matrix_integrated=nothing, psi_nsteps=0, psi_quad_error=0.0, - panel_psis=Float64[], resonance_psis=Float64[]) + panel_psis=Float64[], resonance_psis=Float64[]) end + # Near-axis validity suppression: same boundary and envelope as the EL kinetic + # matrices (one source of truth), applied to the torque density; the quadrature + # domain starts at the boundary since the integrand is identically zero below it. + psi_c = ctrl.axis_validity_suppression ? + kinetic_axis_validity_psi(kinetic_profiles, equil; zi=zi, mi=mi, electron=electron) : 0.0 + x0 = max(x0, psi_c) + if x0 >= xout + return (total=ComplexF64(0.0), torque_profile=nothing, matrix_integrated=nothing, psi_nsteps=0, psi_quad_error=0.0, + panel_psis=Float64[], resonance_psis=Float64[]) + end + psi_c > 0 && @info "Kinetic axis-validity suppression in ψ torque quadrature: domain starts at psi_c=$(round(psi_c; sigdigits=3))" + # The outer ψ-integral (the QuadGK batch / ψ-node loop) stays serial: QuadGK's refine # loop invokes the callback with small batches (~15 nodes), so threading it is # fork-join-bound. Instead thread the inner bounce-harmonic loop (2·nl+1 harmonics), @@ -133,9 +147,9 @@ function integrate_psi_quadgk( w = is_matrix_method ? thread_wtw[tid] : nothing is_matrix_method && fill!(w, 0) tpsi!(thread_tpsi[tid], psi, n, l, zi, mi, wdfac, divxfac, - electron, method, equil, thread_intrs[tid], kinetic_profiles; - op_wmats=w, - atol_xlmda=ctrl.atol_xlmda, rtol_xlmda=ctrl.rtol_xlmda) + electron, method, equil, thread_intrs[tid], kinetic_profiles; + op_wmats=w, + atol_xlmda=ctrl.atol_xlmda, rtol_xlmda=ctrl.rtol_xlmda) harm_vals[ell_idx] = thread_tpsi[tid][] is_matrix_method && (harm_elems[ell_idx] .= w) end @@ -145,6 +159,7 @@ function integrate_psi_quadgk( for ell_idx in 1:nharm total += harm_vals[ell_idx] end + total *= kinetic_axis_validity_envelope(psi, psi_c) y[k] = total push!(logged_psi, psi) @@ -154,6 +169,7 @@ function integrate_psi_quadgk( for ell_idx in 1:nharm elems_accum .+= harm_elems[ell_idx] end + elems_accum .*= kinetic_axis_validity_envelope(psi, psi_c) push!(logged_elems, elems_accum) end end @@ -162,7 +178,7 @@ function integrate_psi_quadgk( # Panels at the rational surfaces the run resolved plus the kinetic-resonance # surfaces (thermal-energy Ω_ℓ = 0 for ℓ ∈ -nl:nl) — both are torque-density peaks. resonance_psis = kinetic_resonance_psi_nodes(kinetic_profiles, equil; n, nl, zi, mi, electron, wdfac) - pts = psi_panel_points(vcat(intr.sing_psis, resonance_psis), x0, xout) + pts = psi_panel_points(vcat(intr.sing_psis, resonance_psis, [2 * psi_c]), x0, xout) bi = QuadGK.BatchIntegrand(psi_batch!, ComplexF64[], Float64[]) total, quad_err = quadgk(bi, pts...; atol=ctrl.atol_psi, rtol=ctrl.rtol_psi, maxevals=ctrl.maxevals_psi) @@ -205,7 +221,7 @@ function integrate_psi_quadgk( "$(length(pts) - 1) panels ($n_rational rational + $n_resonance kinetic resonance surfaces)" return (total=total, torque_profile=torque_profile, matrix_integrated=matrix_integrated, psi_nsteps=npts, psi_quad_error=quad_err, - panel_psis=pts, resonance_psis=sort(resonance_psis)) + panel_psis=pts, resonance_psis=sort(resonance_psis)) end @@ -224,15 +240,16 @@ For multi-n calculations, loops over toroidal mode numbers and assembles block-diagonal kinetic matrices. # Arguments -- `state::KineticForcesState`: Accumulates results for all methods -- `intr::KineticForcesInternal`: Internal state with equilibrium data -- `ctrl::KineticForcesControl`: Control parameters specifying which methods to run -- `equil`: PlasmaEquilibrium with 2D interpolants -- `kinetic_profiles::Equilibrium.KineticProfileSplines`: Named kinetic-profile splines + + - `state::KineticForcesState`: Accumulates results for all methods + - `intr::KineticForcesInternal`: Internal state with equilibrium data + - `ctrl::KineticForcesControl`: Control parameters specifying which methods to run + - `equil`: PlasmaEquilibrium with 2D interpolants + - `kinetic_profiles::Equilibrium.KineticProfileSplines`: Named kinetic-profile splines """ function compute_torque_all_methods!(state::KineticForcesState, intr::KineticForcesInternal, - ctrl::KineticForcesControl, equil, - kinetic_profiles::Equilibrium.KineticProfileSplines) + ctrl::KineticForcesControl, equil, + kinetic_profiles::Equilibrium.KineticProfileSplines) for entry in METHOD_REGISTRY getfield(ctrl, entry.flag) || continue @@ -311,7 +328,7 @@ function compute_torque_all_methods!(state::KineticForcesState, intr::KineticFor t_cumulative=t_cum_out, psi_nsteps=psi_nsteps_total, panel_psis=panel_psis_out, - resonance_psis=resonance_psis_out, + resonance_psis=resonance_psis_out ) state.method_results[method] = result_entry diff --git a/src/KineticForces/KineticForcesStructs.jl b/src/KineticForces/KineticForcesStructs.jl index 0d1b46464..f3e1fbf9b 100644 --- a/src/KineticForces/KineticForcesStructs.jl +++ b/src/KineticForces/KineticForcesStructs.jl @@ -3,15 +3,16 @@ Single source of truth for the NTV calculation methods. Each entry is a NamedTuple `(name, flag, kind, doc)`: -- `name` — short method identifier used as the HDF5 group key and in `intr.method` -- `flag` — the `KineticForcesControl` field symbol that enables the method -- `kind` — dispatch routing tag consumed by `method_kind` / `Torque.jl` - (`:gar` for the GAR/matrix family, `:fcgl`/`:rlar`/`:clar` for the - three special-cased methods) -- `doc` — one-line description printed in verbose output + + - `name` — short method identifier used as the HDF5 group key and in `intr.method` + - `flag` — the `KineticForcesControl` field symbol that enables the method + - `kind` — dispatch routing tag consumed by `method_kind` / `Torque.jl` + (`:gar` for the GAR/matrix family, `:fcgl`/`:rlar`/`:clar` for the + three special-cased methods) + - `doc` — one-line description printed in verbose output The method names/docs and the `Compute.jl` enable list are all derived from this -tuple, and `Torque.jl` routes on `kind`, so the methods are enumerated in one place. +tuple, and `Torque.jl` routes on `kind`, so the methods are enumerated in one place. To add a method: append an entry here and add the matching `*_flag` field to `KineticForcesControl`. """ @@ -56,6 +57,7 @@ User-facing control parameters from the TOML `[KineticForces]` section. Configures which NTV methods to run, species parameters, tolerances, and output options. Constructed via keyword arguments or from a TOML dict: + ```julia ctrl = KineticForcesControl(; (Symbol(k) => v for (k, v) in inputs["KineticForces"])...) ``` @@ -123,6 +125,7 @@ ctrl = KineticForcesControl(; (Symbol(k) => v for (k, v) in inputs["KineticForce nufac::Float64 = 1.0 # Collisionality scaling divxfac::Float64 = 1.0 # div(xi_perp) scaling + axis_validity_suppression::Bool = true # Suppress kinetic terms where the zero-orbit-width ordering fails near the axis (see kinetic_axis_validity_psi); envelope and boundary are profile-derived, no tuning parameters # Energy integration parameters nutype::String = "harmonic" # Collision operator: "zero", "small", "krook", "harmonic" @@ -176,12 +179,13 @@ Internal working state for KineticForces calculations. Holds equilibrium-derived quantities, profile interpolants, and integration results. Fields replacing former module-level globals: -- `ro`, `bo`, `chi1`: Equilibrium geometry parameters -- `mthsurf`, `mfac`: Poloidal grid info -- `dbob_m`, `divx_m`: Perturbation mode interpolants -- `sing_psis`: Rational-surface ψ locations (sorted, from the stability analysis), used as - panel boundaries for the outer ψ torque quadrature so the resonant peaks fall on - Gauss-Kronrod interval endpoints instead of driving deep adaptive bisection + + - `ro`, `bo`, `chi1`: Equilibrium geometry parameters + - `mthsurf`, `mfac`: Poloidal grid info + - `dbob_m`, `divx_m`: Perturbation mode interpolants + - `sing_psis`: Rational-surface ψ locations (sorted, from the stability analysis), used as + panel boundaries for the outer ψ torque quadrature so the resonant peaks fall on + Gauss-Kronrod interval endpoints instead of driving deep adaptive bisection Equilibrium and kinetic profile data are read directly from the `PlasmaEquilibrium` (`equil.profiles`, `equil.geometry`) and the @@ -271,17 +275,17 @@ function KineticForcesInternal(equil; verbose::Bool=false) # Axis toroidal field F(0)/ro that normalizes λ = μ·bo/E; F_spline stores 2πF. bo_axis = abs(equil.profiles.F_spline(0.0)) / (2π * equil.ro) KineticForcesInternal(; - ro = equil.ro, - bo = bo_axis, - chi1 = 2π * equil.psio, + ro=equil.ro, + bo=bo_axis, + chi1=2π * equil.psio, mthsurf, - tpsi_xs = collect(range(0.0, 1.0, length=nth)), - tpsi_B = Vector{Float64}(undef, nth), - tpsi_dBdpsi = Vector{Float64}(undef, nth), - tpsi_dBdtheta = Vector{Float64}(undef, nth), - tpsi_jac = Vector{Float64}(undef, nth), - tpsi_djdpsi = Vector{Float64}(undef, nth), - verbose, + tpsi_xs=collect(range(0.0, 1.0; length=nth)), + tpsi_B=Vector{Float64}(undef, nth), + tpsi_dBdpsi=Vector{Float64}(undef, nth), + tpsi_dBdtheta=Vector{Float64}(undef, nth), + tpsi_jac=Vector{Float64}(undef, nth), + tpsi_djdpsi=Vector{Float64}(undef, nth), + verbose ) end @@ -291,19 +295,21 @@ end Populate perturbation data from PerturbedEquilibriumState into KineticForcesInternal. Builds three interpolant sets from PE Clebsch displacements: -1. `xs_m` — [ξ^ψ, ∂ξ^ψ/∂ψ, ξ^α] CubicSeriesInterpolants over ψ -2. `dbob_m` — δB/B Fourier modes via JBB deweighting (Fortran set_peq) -3. `divx_m` — ∇·ξ⊥ Fourier modes via JBB deweighting + + 1. `xs_m` — [ξ^ψ, ∂ξ^ψ/∂ψ, ξ^α] CubicSeriesInterpolants over ψ + 2. `dbob_m` — δB/B Fourier modes via JBB deweighting (Fortran set_peq) + 3. `divx_m` — ∇·ξ⊥ Fourier modes via JBB deweighting The JBB deweighting algorithm (Fortran pentrc/inputs.f90:828-868): -1. Apply geometric matrices S,T,X,Y,Z in m-space -2. Inverse DFT to θ-space -3. Divide by J·B² at each θ -4. Forward DFT back to m-space + + 1. Apply geometric matrices S,T,X,Y,Z in m-space + 2. Inverse DFT to θ-space + 3. Divide by J·B² at each θ + 4. Forward DFT back to m-space """ function set_perturbation_data!(kf_intr::KineticForcesInternal, pe_state, ffs_intr, - equil::Equilibrium.PlasmaEquilibrium, - metric::ForceFreeStates.MetricData) + equil::Equilibrium.PlasmaEquilibrium, + metric::ForceFreeStates.MetricData) # Copy mode numbers from FFS kf_intr.mlow = ffs_intr.mlow kf_intr.mhigh = ffs_intr.mhigh @@ -372,9 +378,9 @@ function set_perturbation_data!(kf_intr::KineticForcesInternal, pe_state, ffs_in psi = psi_grid[ipsi] # Get Clebsch displacement vectors at this ψ - xsp = view(xi_modes.clebsch_psi, ipsi, :) # ξ^ψ [mpert] + xsp = view(xi_modes.clebsch_psi, ipsi, :) # ξ^ψ [mpert] xmp1 = view(xi_modes.clebsch_psi1, ipsi, :) # ∂ξ^ψ/∂ψ [mpert] - xms = view(clebsch_alpha_mat, ipsi, :) # ξ^α [mpert] + xms = view(clebsch_alpha_mat, ipsi, :) # ξ^α [mpert] # Evaluate geometric matrices at ψ → mpert² flat vectors, reshape to mpert×mpert geom_mats.smats(smat_flat, psi; hint=hint_s) @@ -394,8 +400,8 @@ function set_perturbation_data!(kf_intr::KineticForcesInternal, pe_state, ffs_in mul!(jbb_kapx, smat, xsp) mul!(jbb_kapx, tmat, xms, 1.0 + 0.0im, 1.0 + 0.0im) # += tmat * xms mul!(jbb_divx, xmat, xmp1) - mul!(jbb_divx, ymat, xsp, 1.0 + 0.0im, 1.0 + 0.0im) # += ymat * xsp - mul!(jbb_divx, zmat, xms, 1.0 + 0.0im, 1.0 + 0.0im) # += zmat * xms + mul!(jbb_divx, ymat, xsp, 1.0 + 0.0im, 1.0 + 0.0im) # += ymat * xsp + mul!(jbb_divx, zmat, xms, 1.0 + 0.0im, 1.0 + 0.0im) # += zmat * xms @. jbb_dbob = -(jbb_divx + jbb_kapx) # Inverse DFT to θ-space, divide by J·B², forward DFT back @@ -422,9 +428,9 @@ Matches Fortran set_peq lines 859-868: transforms JBB-weighted m-space data to θ-space, removes the J·B² weighting at each poloidal angle, and transforms back. """ function _jbb_deweight!(out::AbstractVector{ComplexF64}, jbb_modes::Vector{ComplexF64}, - ft::Utilities.FourierTransforms.FourierTransform, - psi::Float64, equil::Equilibrium.PlasmaEquilibrium, - mthsurf::Int, theta_buf::Vector{ComplexF64}) + ft::Utilities.FourierTransforms.FourierTransform, + psi::Float64, equil::Equilibrium.PlasmaEquilibrium, + mthsurf::Int, theta_buf::Vector{ComplexF64}) # Inverse DFT: m-space → θ-space theta_buf .= Utilities.FourierTransforms.inverse(ft, jbb_modes) @@ -494,8 +500,8 @@ Accumulated results from all KineticForces computations. Written to gpec.h5 under the "KineticForces" group. """ @kwdef mutable struct KineticForcesState - method_results::Dict{String, MethodResult} = Dict{String, MethodResult}() + method_results::Dict{String,MethodResult} = Dict{String,MethodResult}() # Block-diagonal kinetic matrices: key=method, value=(numpert_total, numpert_total, 6) - kinetic_matrices::Dict{String, Array{ComplexF64,3}} = Dict{String, Array{ComplexF64,3}}() + kinetic_matrices::Dict{String,Array{ComplexF64,3}} = Dict{String,Array{ComplexF64,3}}() completed::Bool = false end diff --git a/src/KineticForces/Utils.jl b/src/KineticForces/Utils.jl index 51c0cd3a8..97cfa4665 100644 --- a/src/KineticForces/Utils.jl +++ b/src/KineticForces/Utils.jl @@ -71,7 +71,7 @@ the estimate degenerates. Panel placement only needs ~peak-width accuracy, so th estimates (single spline evaluations) are sufficient and no bounce averaging is performed. """ function kinetic_resonance_psi_nodes(kinetic_profiles::Equilibrium.KineticProfileSplines, equil; - n::Int, nl::Int, zi::Int=1, mi::Int=2, electron::Bool=false, wdfac::Float64=1.0, xeval::Float64=2.5) + n::Int, nl::Int, zi::Int=1, mi::Int=2, electron::Bool=false, wdfac::Float64=1.0, xeval::Float64=2.5) chrg = electron ? -e : zi * e mass = electron ? me : mi * mp T_spline = electron ? kinetic_profiles.Te_spline : kinetic_profiles.Ti_spline @@ -89,3 +89,57 @@ function kinetic_resonance_psi_nodes(kinetic_profiles::Equilibrium.KineticProfil grid = filter(x -> x > 0, kinetic_profiles.xs) return _resonance_nodes_from_frequencies(wbhat_f, kinetic_profiles.omegaE_spline, wdhat_f, grid; n, nl, xeval) end + +""" + kinetic_axis_validity_psi(kinetic_profiles, equil; zi=1, mi=2, electron=false) → Float64 + +ψ_N below which the zero-orbit-width drift-kinetic ordering fails: the outermost ψ on the +kinetic-profile grid where any of the three thermal orbit-width scales reaches the local minor +radius ⟨r⟩ — potato width `(q²ρ²R₀)^(1/3)`, banana width `q·ρ/√ε`, and poloidal gyroradius +`q·ρ/ε` (ε = ⟨r⟩/⟨R⟩, clamped as in `kinetic_resonance_psi_nodes`; thermal gyroradius +ρ = √(2·m·T)·/(Z·e·B₀) of the computed species). Inside this boundary trapped bananas become +potato orbits with width comparable to r itself, so the bounce-averaged kinetic response is +evaluated outside its validity domain (measured consequence: diverging kinetic increments and a +pathological EL step count). Returns 0.0 when no criterion is met anywhere. The Fortran precedent +(`ktanh_flag`, dcon/fourfit.F) suppressed the same region with four hand-tuned knobs; here the +boundary is derived from the profiles with no user parameters. +""" +function kinetic_axis_validity_psi(kinetic_profiles::Equilibrium.KineticProfileSplines, equil; + zi::Int=1, mi::Int=2, electron::Bool=false) + chrg = electron ? e : zi * e + mass = electron ? me : mi * mp + T_spline = electron ? kinetic_profiles.Te_spline : kinetic_profiles.Ti_spline + q_spline = equil.profiles.q_spline + avg_r = equil.geometry.avg_r_spline + avg_R = equil.geometry.avg_R_spline + ro = abs(equil.ro) + bo = abs(equil.params.b0) + psi_c = 0.0 + for psi in kinetic_profiles.xs + psi <= 0 && continue + r = avg_r(psi) + r <= 0 && continue + eps = max(r / avg_R(psi), 1e-6) + rho = mass * sqrt(2 * T_spline(psi) / mass) / (abs(chrg) * bo) + q = abs(q_spline(psi)) + w_orbit = max(cbrt(q^2 * rho^2 * ro), q * rho / sqrt(eps), q * rho / eps) + w_orbit >= r && (psi_c = max(psi_c, psi)) + end + return psi_c +end + +""" + kinetic_axis_validity_envelope(psi, psi_c) → Float64 + +C² quintic smoothstep for the near-axis kinetic suppression: 0 for ψ ≤ ψ_c (drift-kinetic model +invalid; kernel evaluation may be skipped), rising over `[ψ_c, 2ψ_c]`, 1 above. The transition +width is tied to ψ_c itself, so there is no independent width parameter. `psi_c ≤ 0` returns 1 +(no suppression). +""" +function kinetic_axis_validity_envelope(psi::Float64, psi_c::Float64) + psi_c <= 0 && return 1.0 + t = (psi - psi_c) / psi_c + t <= 0 && return 0.0 + t >= 1 && return 1.0 + return t^3 * (10 + t * (6 * t - 15)) +end From f95230d38c1e5e0e7b8d61712f49fcbc3ae5b047 Mon Sep 17 00:00:00 2001 From: logan-nc Date: Thu, 20 Aug 2026 11:56:29 -0400 Subject: [PATCH 2/4] KF - NEW FEATURE - Output kinetic-model validity profiles under KineticForces/Validity Whenever kinetic profiles are used (self-consistent matrices or NTV post-processing), write the thermal orbit-width scales (rho_i, rho_banana, rho_theta, w_potato), the local geometry (r_minor, d_separatrix), the profile gradient lengths (L_p, L_q), the near-axis boundary psi_c with its applied envelope, and an is_valid array (orbit width < r, rho_banana < L_p and L_q, orbit width < distance to separatrix). Validity outside the near-axis envelope is flagged, never suppressed -- the far edge can dominate the physical NTV. Co-Authored-By: Claude Fable 5 Claude-Session: https://claude.ai/code/session_01LzbLFQKyuRE5DYZmLokKmk --- src/GeneralizedPerturbedEquilibrium.jl | 39 +++++++++-- src/HDF5Schema.jl | 38 +++++++++-- .../CalculatedKineticMatrices.jl | 2 +- src/KineticForces/Utils.jl | 67 +++++++++++++++++++ 4 files changed, 132 insertions(+), 14 deletions(-) diff --git a/src/GeneralizedPerturbedEquilibrium.jl b/src/GeneralizedPerturbedEquilibrium.jl index 7b17683cf..d1d85d22c 100755 --- a/src/GeneralizedPerturbedEquilibrium.jl +++ b/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,10 +419,11 @@ 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) make_kinetic_matrix(ctrl, equil, ffit, intr, metric; calculated_source=calculated_cb) @@ -495,6 +496,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 +708,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 +826,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. diff --git a/src/HDF5Schema.jl b/src/HDF5Schema.jl index a1a47ef44..06487b7ea 100644 --- a/src/HDF5Schema.jl +++ b/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⟩"), + "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 + ] ) """ diff --git a/src/KineticForces/CalculatedKineticMatrices.jl b/src/KineticForces/CalculatedKineticMatrices.jl index 9d8a19def..492a62403 100644 --- a/src/KineticForces/CalculatedKineticMatrices.jl +++ b/src/KineticForces/CalculatedKineticMatrices.jl @@ -108,7 +108,7 @@ function compute_calculated_kinetic_matrices( 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)); $(count(iszero, env)) of $mpsi surfaces skipped" + "psi=$(round(2 * psi_c; sigdigits=3)) (kernel evaluation skipped below psi_c)" maxlog = 1 end end diff --git a/src/KineticForces/Utils.jl b/src/KineticForces/Utils.jl index 97cfa4665..e2ab60c43 100644 --- a/src/KineticForces/Utils.jl +++ b/src/KineticForces/Utils.jl @@ -143,3 +143,70 @@ function kinetic_axis_validity_envelope(psi::Float64, psi_c::Float64) t >= 1 && return 1.0 return t^3 * (10 + t * (6 * t - 15)) end + +""" + kinetic_validity_profiles(kinetic_profiles, equil; zi=1, mi=2, electron=false) → NamedTuple + +Radial profiles of the drift-kinetic validity diagnostics, on the kinetic-profile ψ grid: +`psi`; the thermal orbit-width scales `rho_i` (gyroradius √(2mT)/(Z·e·B₀)), `rho_banana` +(q·ρ/√ε), `rho_theta` (poloidal gyroradius q·ρ/ε), `w_potato` ((q²ρ²R₀)^(1/3)); the local +geometry `r_minor` (⟨r⟩) and `d_separatrix` (⟨r⟩(1) − ⟨r⟩(ψ)); the profile gradient lengths +`L_p` and `L_q` (|X|·|dr/dψ|/|dX/dψ|, Inf where the profile is flat); the near-axis suppression +boundary `psi_c` (`kinetic_axis_validity_psi`) with its `envelope`; and `is_valid` — true where +every zero-orbit-width ordering holds at coefficient 1: max orbit width < r, ρ_banana < L_p and +L_q, and max(ρ_i, ρ_banana, ρ_θ) < d_separatrix. Validity is diagnostic only — nothing outside +the near-axis envelope is suppressed (the far edge and steep-gradient regions are flagged, not +zeroed, since they can dominate the physical NTV). +""" +function kinetic_validity_profiles(kinetic_profiles::Equilibrium.KineticProfileSplines, equil; + zi::Int=1, mi::Int=2, electron::Bool=false) + chrg = electron ? e : zi * e + mass = electron ? me : mi * mp + T_spline = electron ? kinetic_profiles.Te_spline : kinetic_profiles.Ti_spline + q_spline, q_deriv = equil.profiles.q_spline, equil.profiles.q_deriv + P_spline, P_deriv = equil.profiles.P_spline, equil.profiles.P_deriv + avg_r = equil.geometry.avg_r_spline + avg_R = equil.geometry.avg_R_spline + r_deriv = deriv1(avg_r) + ro = abs(equil.ro) + bo = abs(equil.params.b0) + r_sep = avg_r(1.0) + psi_c = kinetic_axis_validity_psi(kinetic_profiles, equil; zi=zi, mi=mi, electron=electron) + + psi = [x for x in kinetic_profiles.xs if 0 < x <= 1] + n = length(psi) + rho_i = zeros(n) + rho_banana = zeros(n) + rho_theta = zeros(n) + w_potato = zeros(n) + r_minor = zeros(n) + L_p = zeros(n) + L_q = zeros(n) + d_separatrix = zeros(n) + envelope = zeros(n) + is_valid = falses(n) + for (i, x) in pairs(psi) + r = avg_r(x) + eps = max(r / avg_R(x), 1e-6) + rho = mass * sqrt(2 * T_spline(x) / mass) / (abs(chrg) * bo) + q = abs(q_spline(x)) + drdpsi = abs(r_deriv(x)) + rho_i[i] = rho + rho_banana[i] = q * rho / sqrt(eps) + rho_theta[i] = q * rho / eps + w_potato[i] = cbrt(q^2 * rho^2 * ro) + r_minor[i] = r + d_separatrix[i] = max(r_sep - r, 0.0) + dP = abs(P_deriv(x)) + L_p[i] = dP > 0 ? abs(P_spline(x)) * drdpsi / dP : Inf + dq = abs(q_deriv(x)) + L_q[i] = dq > 0 ? q * drdpsi / dq : Inf + envelope[i] = kinetic_axis_validity_envelope(x, psi_c) + w_orbit = max(w_potato[i], rho_banana[i], rho_theta[i]) + is_valid[i] = w_orbit < r && rho_banana[i] < L_p[i] && rho_banana[i] < L_q[i] && + max(rho_i[i], rho_banana[i], rho_theta[i]) < d_separatrix[i] + end + return (psi=psi, rho_i=rho_i, rho_banana=rho_banana, rho_theta=rho_theta, w_potato=w_potato, + r_minor=r_minor, L_p=L_p, L_q=L_q, d_separatrix=d_separatrix, + psi_c=psi_c, envelope=envelope, is_valid=is_valid) +end From 3f1dbfef948d2b07511f9195e412bf0710767d95 Mon Sep 17 00:00:00 2001 From: logan-nc Date: Thu, 20 Aug 2026 12:45:05 -0400 Subject: [PATCH 3/4] KF - BUGFIX - Resolve the validity-envelope band on coarse grids; open-or-create the KineticForces group The envelope has structure on the psi_c scale; coarse kinetic decks (m16) cannot represent env*(increment) and the spline overshoot can land on a rational surface inside the transition band, corrupting the eigenvalues. Augment the kernel grid with knots across [psi_c, 2 psi_c] (band ends pinned -- the smoothstep is only C2 there) on the full-grid path and seed them on the certified path. Also open-or-create KineticForces in the NTV writer, which collided with the Validity group. Co-Authored-By: Claude Fable 5 Claude-Session: https://claude.ai/code/session_01LzbLFQKyuRE5DYZmLokKmk --- src/ForceFreeStates/Kinetic.jl | 27 ++++++++++++++++++++------ src/GeneralizedPerturbedEquilibrium.jl | 10 +++++++++- src/KineticForces/Output.jl | 25 +++++++++++++----------- 3 files changed, 44 insertions(+), 18 deletions(-) diff --git a/src/ForceFreeStates/Kinetic.jl b/src/ForceFreeStates/Kinetic.jl index da1202bae..fa9f954bb 100644 --- a/src/ForceFreeStates/Kinetic.jl +++ b/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 diff --git a/src/GeneralizedPerturbedEquilibrium.jl b/src/GeneralizedPerturbedEquilibrium.jl index d1d85d22c..52d255a07 100755 --- a/src/GeneralizedPerturbedEquilibrium.jl +++ b/src/GeneralizedPerturbedEquilibrium.jl @@ -424,8 +424,16 @@ function main_from_inputs( 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; diff --git a/src/KineticForces/Output.jl b/src/KineticForces/Output.jl index dc613c04f..975ed6f03 100644 --- a/src/KineticForces/Output.jl +++ b/src/KineticForces/Output.jl @@ -12,14 +12,15 @@ then write to gpec.h5 in a single pass. Write all KineticForces results to the "KineticForces" group in gpec.h5. # Arguments -- `h5file::HDF5.File`: Open HDF5 file handle -- `state::KineticForcesState`: Accumulated computation results -- `dVdpsi_spline`: Optional dV/dψ_N profile interpolant; when given, dV/dψ_N is - written at the quadrature points so the torque density dT/dV = (dT/dψ)/(dV/dψ) - is directly available + + - `h5file::HDF5.File`: Open HDF5 file handle + - `state::KineticForcesState`: Accumulated computation results + - `dVdpsi_spline`: Optional dV/dψ_N profile interpolant; when given, dV/dψ_N is + written at the quadrature points so the torque density dT/dV = (dT/dψ)/(dV/dψ) + is directly available """ function write_to_hdf5!(h5file::HDF5.File, state::KineticForcesState; dVdpsi_spline=nothing) - g = create_group(h5file, "KineticForces") + g = haskey(h5file, "KineticForces") ? h5file["KineticForces"] : create_group(h5file, "KineticForces") for (method_name, result) in state.method_results mg = create_group(g, method_name) @@ -113,8 +114,9 @@ Write variable-length integration trajectory records using offset-indexed concat This is the standard HDF5 ragged array pattern for storing variable-length data. # Arguments -- `mg::HDF5.Group`: HDF5 group for this method -- `records::Vector{EnergyIntegrationResult}`: Integration records to write + + - `mg::HDF5.Group`: HDF5 group for this method + - `records::Vector{EnergyIntegrationResult}`: Integration records to write """ function write_integration_records!(mg::HDF5.Group, records::Vector{EnergyIntegrationResult}) rg = create_group(mg, "EnergyIntegrals") @@ -146,13 +148,14 @@ end Print a summary of KineticForces results to stdout. # Arguments -- `state::KineticForcesState`: Accumulated computation results -- `verbose::Bool`: Print detailed per-surface results + + - `state::KineticForcesState`: Accumulated computation results + - `verbose::Bool`: Print detailed per-surface results """ function print_summary(state::KineticForcesState; verbose::Bool=false) for (method_name, result) in state.method_results @printf("%-8s T_phi = %11.3e 2n*dW_k = %11.3e\n", - method_name, real(result.total_torque), imag(result.total_torque)) + method_name, real(result.total_torque), imag(result.total_torque)) end if verbose for (method_name, _) in state.kinetic_matrices From 12fca467cc7c9e0aeb4022c42a7caedc90fda655 Mon Sep 17 00:00:00 2001 From: logan-nc Date: Thu, 20 Aug 2026 13:41:47 -0400 Subject: [PATCH 4/4] KF - MINOR - Add units to the Validity psi_c annotation (schema metadata contract) Co-Authored-By: Claude Fable 5 Claude-Session: https://claude.ai/code/session_01LzbLFQKyuRE5DYZmLokKmk --- src/HDF5Schema.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/HDF5Schema.jl b/src/HDF5Schema.jl index 06487b7ea..e1a73bb4e 100644 --- a/src/HDF5Schema.jl +++ b/src/HDF5Schema.jl @@ -151,7 +151,7 @@ const MAIN_H5_ANNOTATIONS = [ "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⟩"), + "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" => (;