From 76370a6657e1b3383bb66914cb91ffb7f959dd33 Mon Sep 17 00:00:00 2001 From: logan-nc Date: Wed, 29 Jul 2026 10:54:30 -0400 Subject: [PATCH 01/17] KineticForces - BUG FIX - Coulomb logarithm uses natural log, not log10 The lnLambda formula 17.3 - 0.5*log(ne/1e20) + 1.5*log(Te/1keV) is the NRL Plasma Formulary electron-ion Coulomb logarithm (lambda_ei = 23 - ln(ne^1/2 * Te^-3/2)) re-normalized to (ne/1e20 m^-3, Te/1keV); its 17.3/-0.5/+1.5 coefficients are calibrated for natural log. The Julia port used log10, so lnLambda -- and the collision frequencies nue/nui that scale linearly with it -- ran ~13% low at a 20 keV core (19.2 vs the correct 21.7). Restore natural log, matching both Fortran PENTRC inputs.f90:238 and the formulary. Regression (diiid_n1): NTV FGAR and dW FGAR shift within tolerance (~0.06-0.2%; small only because that case sits near the ne=1e20/Te=1keV reference point where log10 ~= ln). Co-Authored-By: Claude Opus 4.8 (1M context) Claude-Session: https://claude.ai/code/session_01SsLUydP2gJbE1tGoaHDiS1 --- src/Equilibrium/KineticProfiles.jl | 5 ++++- 1 file changed, 4 insertions(+), 1 deletion(-) diff --git a/src/Equilibrium/KineticProfiles.jl b/src/Equilibrium/KineticProfiles.jl index 606e2449..65a105b2 100644 --- a/src/Equilibrium/KineticProfiles.jl +++ b/src/Equilibrium/KineticProfiles.jl @@ -267,7 +267,10 @@ function load_kinetic_profiles(kinetic_file::AbstractString; z = n_e > 0 ? zimp - (n_i / n_e) * zi * (zimp - zi) : Float64(zimp) zpitch = 1.0 + (1.0 + mimp) / (2.0 * mimp) * zimp * (z - 1.0) / (zimp - z) - ll = 17.3 - 0.5 * log10(n_e / 1.0e20) + 1.5 * log10(T_e / 1.602e-16) + # Coulomb logarithm — natural log, matching Fortran PENTRC inputs.f90:238 (the + # coefficients 17.3/-0.5/+1.5 are calibrated for ln, not log10; n_e in units of 1e20 m^-3, + # T_e in units of 1 keV = 1.602e-16 J). + ll = 17.3 - 0.5 * log(n_e / 1.0e20) + 1.5 * log(T_e / 1.602e-16) loglam[i] = ll nui[i] = T_i > 0 ? From ee5d211d43b9ebb76b73229db223a9be7af6e5f7 Mon Sep 17 00:00:00 2001 From: logan-nc Date: Wed, 29 Jul 2026 11:05:10 -0400 Subject: [PATCH 02/17] KineticForces - NEW FEATURE - ion_fraction for multi-main-ion NTV (fgar path) First increment toward correct multi-main-ion (e.g. D-T) NTV. Adds an ion_fraction control (default 1.0 = single main ion, bit-for-bit backward compatible) that scales only this species' resonant density out of the total main-ion n_i. The Coulomb collisionality and Zeff use the full n_i, so a multi-ion run supplies the TOTAL main-ion density in the kinetic file and runs each species (D, T) with its fraction; wdian/wdiat are unaffected (the fraction cancels in T*(dn/dpsi)/n), and nu_s stays the full-composition collisionality. Vetted against Logan-Park PoP 2013 (fortran-physics-reviewer). Wired through the fgar NTV torque path (Compute.integrate_psi_quadgk -> tpsi!) and auto-exposed via the [KineticForces] TOML splat. Not yet threaded (tracked in the feature plan): the matrix-method (_setup_surface_state) and self-consistent kinetic-DCON paths use the implicit 1.0; and an internal multi-species loop (D/T summation is still external for now, but with the correct Zeff). Co-Authored-By: Claude Opus 4.8 (1M context) Claude-Session: https://claude.ai/code/session_01SsLUydP2gJbE1tGoaHDiS1 --- src/KineticForces/Compute.jl | 2 +- src/KineticForces/KineticForcesStructs.jl | 6 ++++++ src/KineticForces/Torque.jl | 8 ++++++-- 3 files changed, 13 insertions(+), 3 deletions(-) diff --git a/src/KineticForces/Compute.jl b/src/KineticForces/Compute.jl index 4175a62b..f2b3f5b9 100644 --- a/src/KineticForces/Compute.jl +++ b/src/KineticForces/Compute.jl @@ -134,7 +134,7 @@ function integrate_psi_quadgk( 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, + op_wmats=w, ion_fraction=ctrl.ion_fraction, 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) diff --git a/src/KineticForces/KineticForcesStructs.jl b/src/KineticForces/KineticForcesStructs.jl index 7e533d8e..0183e9d5 100644 --- a/src/KineticForces/KineticForcesStructs.jl +++ b/src/KineticForces/KineticForcesStructs.jl @@ -90,6 +90,12 @@ ctrl = KineticForcesControl(; (Symbol(k) => v for (k, v) in inputs["KineticForce zimp::Int = 6 # Impurity charge mimp::Int = 12 # Impurity mass electron::Bool = false # Include electron contribution + # Fraction of the main-ion density carried by THIS ion species, for a multi-main-ion + # plasma (e.g. 50/50 D-T: run D and T each with ion_fraction=0.5). Scales only the + # resonant density prefactor (n_s = ion_fraction·n_i); the Coulomb collisionality and + # Zeff use the full n_i, so the kinetic file's n_i column must be the TOTAL main-ion + # density for a multi-ion run. Default 1.0 reproduces the single-main-ion behaviour. + ion_fraction::Float64 = 1.0 # Mode numbers nn::Int = 1 # Toroidal mode number diff --git a/src/KineticForces/Torque.jl b/src/KineticForces/Torque.jl index af22d165..82de0045 100644 --- a/src/KineticForces/Torque.jl +++ b/src/KineticForces/Torque.jl @@ -35,6 +35,7 @@ function tpsi!(tpsi_var::Ref{ComplexF64}, psi::Float64, n::Int, l::Int, op_wmats::Union{Nothing,Array{ComplexF64,3}}=nothing, rex_override::Union{Nothing,Float64}=nothing, imx_override::Union{Nothing,Float64}=nothing, + ion_fraction::Float64=1.0, atol_xlmda::Float64=1e-9, rtol_xlmda::Float64=1e-6) # Enforce bounds @@ -153,9 +154,12 @@ function tpsi!(tpsi_var::Ref{ComplexF64}, psi::Float64, n::Int, l::Int, dT_s_dpsi = kinetic_profiles.Te_deriv(psi) nu_s = kinetic_profiles.nue_spline(psi) else - n_s = kinetic_profiles.ni_spline(psi) + # ion_fraction scales this species' resonant density out of the total main-ion n_i + # (e.g. 0.5 for each of D/T). It cancels in wdian/wdiat (T·(dn/dψ)/n) and does NOT + # scale nu_s, which is the full-composition collisionality. + n_s = ion_fraction * kinetic_profiles.ni_spline(psi) T_s = kinetic_profiles.Ti_spline(psi) - dn_s_dpsi = kinetic_profiles.ni_deriv(psi) + dn_s_dpsi = ion_fraction * kinetic_profiles.ni_deriv(psi) dT_s_dpsi = kinetic_profiles.Ti_deriv(psi) nu_s = kinetic_profiles.nui_spline(psi) end From a3760fc4c6043ae8b92df8b4e2acfb4f2d496a63 Mon Sep 17 00:00:00 2001 From: logan-nc Date: Wed, 29 Jul 2026 11:53:18 -0400 Subject: [PATCH 03/17] KineticForces - NEW FEATURE - multi_ion_composition (Zeff/zpitch for arbitrary ion mix) Physics core for correct multi-main-ion NTV: full-composition Zeff and momentum-restoring pitch-angle enhancement from an arbitrary list of main-ion (z_s, n_s) plus one impurity that closes quasineutrality (n_imp = (ne - Sum z_s n_s)/zimp; Zeff = (Sum z_s^2 n_s + zimp^2 n_imp)/ne). Verified: reduces EXACTLY to the current single-ion formula Zeff = zimp - (ni/ne) zi (zimp-zi) (Zeff 1.44355 both ways); gives the correct 50/50 D-T Zeff=1.44 vs the wrong 3.72 that each split single-ion run currently sees (the ni double-duty bug). Per-species collisionality built on this shares zpitch/n_main/Zeff/lnLambda and varies only z_s^2, m_s, T_s (Logan-Park PoP 2013, fortran-physics-reviewer). Pure function, not yet wired into the loader/driver. Co-Authored-By: Claude Opus 4.8 (1M context) Claude-Session: https://claude.ai/code/session_01SsLUydP2gJbE1tGoaHDiS1 --- src/Equilibrium/KineticProfiles.jl | 32 ++++++++++++++++++++++++++++++ 1 file changed, 32 insertions(+) diff --git a/src/Equilibrium/KineticProfiles.jl b/src/Equilibrium/KineticProfiles.jl index 65a105b2..a10e6439 100644 --- a/src/Equilibrium/KineticProfiles.jl +++ b/src/Equilibrium/KineticProfiles.jl @@ -128,6 +128,38 @@ function write_kinetic_h5(path::AbstractString, data::KineticProfileData; return path end +""" + multi_ion_composition(zs, ns, ne; zimp, mimp) -> (zeff, zpitch, n_main, n_imp) + +Full-composition effective charge and pitch-angle enhancement for a plasma with an +arbitrary number of main-ion species. Inputs are per-point (scalars or broadcast arrays): + + - `zs`, `ns`: main-ion charges and number densities (one entry per species) + - `ne`: electron density + - `zimp`, `mimp`: single trailing impurity charge / mass; its density closes quasineutrality + +Composition (per point): + + n_main = Σ_s n_s total main-ion number density + n_imp = (n_e − Σ_s z_s n_s) / z_imp impurity density from quasineutrality + Zeff = (Σ_s z_s² n_s + z_imp² n_imp) / n_e + zpitch = 1 + (1+m_imp)/(2 m_imp)·z_imp·(Zeff−1)/(z_imp−Zeff) momentum-restoring closure + +Reduces **exactly** to the single-ion form `Zeff = z_imp − (n_i/n_e)·z_i·(z_imp−z_i)` for one +z-species (verified algebraically). The per-species pitch-angle collision frequency built from +this is `ν_s = (zpitch/3.5e17)·z_s²·n_main·lnΛ / (√m_s · (T_s)^{3/2})` — shared `zpitch`, `n_main`, +`Zeff`, `lnΛ`; the test species contributes only its own `z_s²`, `m_s`, `T_s`. +""" +function multi_ion_composition(zs::AbstractVector, ns::AbstractVector, ne::Real; + zimp::Real, mimp::Real) + n_main = sum(ns) + charge_density = sum(z * n for (z, n) in zip(zs, ns)) # Σ z_s n_s + n_imp = ne > 0 ? (ne - charge_density) / zimp : 0.0 + zeff = ne > 0 ? (sum(z^2 * n for (z, n) in zip(zs, ns)) + zimp^2 * n_imp) / ne : Float64(zimp) + zpitch = 1.0 + (1.0 + mimp) / (2.0 * mimp) * zimp * (zeff - 1.0) / (zimp - zeff) + return (zeff, zpitch, n_main, n_imp) +end + """ load_kinetic_profiles(kinetic_file::AbstractString; zi::Int=1, zimp::Int=6, mi::Int=2, mimp::Int=12, From 09f7d9c389cdace385cbddf9acb4605d1634a04e Mon Sep 17 00:00:00 2001 From: logan-nc Date: Wed, 29 Jul 2026 12:18:53 -0400 Subject: [PATCH 04/17] KineticForces - NEW FEATURE - IonSpecies list + per-species profile builder Replaces the scalar ion_fraction stopgap with the species-list design. Adds: - IonSpecies(z, m, fraction|density) + ion_species::Vector on KineticForcesControl, auto-converted from [[KineticForces.ion_species]] TOML tables (Vector{Dict}->Vector{IonSpecies}). Empty default => single ion from zi/mi (unchanged). - Equilibrium.build_species_profiles(file, ion_species; zimp, mimp) -> Vector{KineticProfileSplines}, one per species: ni_spline = that species' resonant density (fraction * total n_i), nui_spline = its full-composition collisionality (shared Zeff/zpitch/lnLambda via multi_ion_composition, per-species z_s^2, m_s, T_i); ne/Te/omegaE/zeff/nue shared. tpsi! consumes each view unchanged. Verified: single-species [z=1,m=2,fraction=1] reproduces load_kinetic_profiles exactly (ni/nui/zeff); 50/50 D-T gives n_D=n_T=0.5 n_i, shared Zeff=1.46 (not the buggy 3.7), nu_D/nu_T=sqrt(3/2). Reverts the ion_fraction wiring in tpsi!/Compute. Explicit per-species profiles (density=) validated but not yet wired (fraction path only) -- next increment. Co-Authored-By: Claude Opus 4.8 (1M context) Claude-Session: https://claude.ai/code/session_01SsLUydP2gJbE1tGoaHDiS1 --- src/Equilibrium/KineticProfiles.jl | 70 +++++++++++++++++++++++ src/KineticForces/Compute.jl | 2 +- src/KineticForces/KineticForcesStructs.jl | 37 ++++++++++-- src/KineticForces/Torque.jl | 10 ++-- 4 files changed, 106 insertions(+), 13 deletions(-) diff --git a/src/Equilibrium/KineticProfiles.jl b/src/Equilibrium/KineticProfiles.jl index a10e6439..129b8232 100644 --- a/src/Equilibrium/KineticProfiles.jl +++ b/src/Equilibrium/KineticProfiles.jl @@ -160,6 +160,76 @@ function multi_ion_composition(zs::AbstractVector, ns::AbstractVector, ne::Real; return (zeff, zpitch, n_main, n_imp) end +""" + build_species_profiles(kinetic_file, ion_species; zimp, mimp, ...) -> Vector{KineticProfileSplines} + +Build one [`KineticProfileSplines`](@ref) per main-ion species for a multi-ion NTV run. Each +returned view carries **this species' resonant density** `n_s` (as `ni_spline`) and its +**full-composition collision frequency** `ν_s` (as `nui_spline`), while `ne/Te/ωE/Zeff/lnΛ` +(and the electron `nue`) are shared across all species. Downstream `tpsi!` uses each view +unchanged; the driver loops over the returned vector and sums the torques. + +Species density is resolved from `IonSpecies`: `fraction` ⇒ `n_s = fraction · n_i` (the kinetic +file's `n_i` = total main-ion density); explicit per-species profiles (`density`) are not yet +wired here. Composition (Zeff, zpitch) comes from [`multi_ion_composition`]; the per-species +collisionality is `ν_s = (zpitch/3.5e17)·z_s²·n_main·lnΛ / (√m_s·(T_i)^{3/2})`. Reduces to the +single-ion `load_kinetic_profiles` result for one z=1 species with `fraction=1`. +""" +function build_species_profiles(kinetic_file::AbstractString, ion_species::AbstractVector; + zimp::Integer=6, mimp::Integer=12, + density_factor::Float64=1.0, temperature_factor::Float64=1.0, + ExB_rotation_factor::Float64=1.0, toroidal_rotation_factor::Float64=1.0) + + (density_factor == 1.0 && temperature_factor == 1.0 && ExB_rotation_factor == 1.0 && + toroidal_rotation_factor == 1.0) || + error("profile scaling factors are not yet supported together with a multi-ion ion_species list") + + data = read_kinetic_file(kinetic_file) + _need(f, n) = f === nothing ? error("kinetic file '$kinetic_file' missing required '$n'") : f + psi_in = data.psi + ne_in = _need(data.n_e, "n_e"); Ti_in = _need(data.T_i, "T_i") + Te_in = _need(data.T_e, "T_e"); omE_in = _need(data.omega_E, "omega_E") + ni_in = data.n_i === nothing ? copy(ne_in) : data.n_i # TOTAL main-ion density + + eV_to_J = 1.602e-19 + mp = 1.672_614e-27; me = 9.109_1e-31 + nkin = 100; psi_reg = collect(0:nkin) ./ nkin + ne = _cubic_resample(psi_in, ne_in, psi_reg) + ni_total = _cubic_resample(psi_in, ni_in, psi_reg) + Ti = _cubic_resample(psi_in, Ti_in, psi_reg) .* eV_to_J + Te = _cubic_resample(psi_in, Te_in, psi_reg) .* eV_to_J + omegaE = _cubic_resample(psi_in, omE_in, psi_reg) + + # Resolve each species' resonant density (fraction path). + zs = [Int(s.z) for s in ion_species] + ns = Vector{Vector{Float64}}(undef, length(ion_species)) + for (si, s) in enumerate(ion_species) + has_frac = !isnan(s.fraction) + has_dens = !isempty(s.density) + (has_frac ⊻ has_dens) || error("ion_species[$si]: specify exactly one of `fraction` or `density`") + has_dens && error("ion_species[$si]: explicit per-species density `$(s.density)` not yet wired; use `fraction`") + ns[si] = s.fraction .* ni_total + end + + # Shared composition + collisionality (natural-log Coulomb log), per grid point. + npts = length(psi_reg) + zeff = zeros(npts); zpitch = zeros(npts); loglam = zeros(npts); n_main = zeros(npts); nue = zeros(npts) + for i in 1:npts + z, zp, nm, _ = multi_ion_composition(zs, [ns[si][i] for si in eachindex(ion_species)], ne[i]; zimp=zimp, mimp=mimp) + zeff[i] = z; zpitch[i] = zp; n_main[i] = nm + loglam[i] = (ne[i] > 0 && Te[i] > 0) ? 17.3 - 0.5 * log(ne[i] / 1.0e20) + 1.5 * log(Te[i] / 1.602e-16) : 0.0 + nue[i] = Te[i] > 0 ? (zp / 3.5e17) * ne[i] * loglam[i] / (sqrt(me / mp) * (Te[i] / 1.602e-16)^1.5) : 0.0 + end + + # One KineticProfileSplines per species: ni = n_s, nui = ν_s, everything else shared. + profs = Vector{KineticProfileSplines}(undef, length(ion_species)) + for (si, s) in enumerate(ion_species) + nu_s = [Ti[i] > 0 ? (zpitch[i] / 3.5e17) * s.z^2 * n_main[i] * loglam[i] / (sqrt(Float64(s.m)) * (Ti[i] / 1.602e-16)^1.5) : 0.0 for i in 1:npts] + profs[si] = KineticProfileSplines(psi_reg, ns[si], ne, Ti, Te, omegaE, loglam, nu_s, nue, zeff) + end + return profs +end + """ load_kinetic_profiles(kinetic_file::AbstractString; zi::Int=1, zimp::Int=6, mi::Int=2, mimp::Int=12, diff --git a/src/KineticForces/Compute.jl b/src/KineticForces/Compute.jl index f2b3f5b9..4175a62b 100644 --- a/src/KineticForces/Compute.jl +++ b/src/KineticForces/Compute.jl @@ -134,7 +134,7 @@ function integrate_psi_quadgk( 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, ion_fraction=ctrl.ion_fraction, + 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) diff --git a/src/KineticForces/KineticForcesStructs.jl b/src/KineticForces/KineticForcesStructs.jl index 0183e9d5..d2c2e835 100644 --- a/src/KineticForces/KineticForcesStructs.jl +++ b/src/KineticForces/KineticForcesStructs.jl @@ -49,6 +49,31 @@ function method_kind(name::AbstractString) end +""" + IonSpecies(; z, m, fraction=NaN, density="") + +One main-ion species in a multi-ion NTV run. `z`/`m` are the charge (e) and mass (proton +masses). The density is given by **exactly one** of: + + - `fraction`: this species' share of the total main-ion density (the kinetic file's `n_i` + column), e.g. 0.5 for each of 50/50 D-T. Convenient when species share a profile shape. + - `density`: the name of an explicit per-species density profile in the (HDF5) kinetic file, + for measured profiles with distinct shapes. + +Both resolve to a per-species density profile `n_s(ψ)` at load time (one downstream path). +Constructed from a TOML `[[KineticForces.ion_species]]` table via [`convert`]. +""" +struct IonSpecies + z::Int + m::Int + fraction::Float64 + density::String +end +IonSpecies(; z::Integer, m::Integer, fraction=NaN, density="") = + IonSpecies(Int(z), Int(m), Float64(fraction), String(density)) +IonSpecies(d::AbstractDict) = IonSpecies(; (Symbol(k) => v for (k, v) in d)...) +Base.convert(::Type{Vector{IonSpecies}}, v::AbstractVector{<:AbstractDict}) = IonSpecies.(v) + """ KineticForcesControl @@ -90,12 +115,12 @@ ctrl = KineticForcesControl(; (Symbol(k) => v for (k, v) in inputs["KineticForce zimp::Int = 6 # Impurity charge mimp::Int = 12 # Impurity mass electron::Bool = false # Include electron contribution - # Fraction of the main-ion density carried by THIS ion species, for a multi-main-ion - # plasma (e.g. 50/50 D-T: run D and T each with ion_fraction=0.5). Scales only the - # resonant density prefactor (n_s = ion_fraction·n_i); the Coulomb collisionality and - # Zeff use the full n_i, so the kinetic file's n_i column must be the TOTAL main-ion - # density for a multi-ion run. Default 1.0 reproduces the single-main-ion behaviour. - ion_fraction::Float64 = 1.0 + # Multi-main-ion plasma: list of IonSpecies, each with its own z, m, and density + # (fraction of the total main-ion n_i, or an explicit per-species profile). When + # non-empty the NTV is computed per species with one shared full-composition Zeff and + # summed. Empty (default) ⇒ single main ion from zi/mi (unchanged behaviour). The + # kinetic file's n_i column is the TOTAL main-ion density for a multi-ion run. + ion_species::Vector{IonSpecies} = IonSpecies[] # Mode numbers nn::Int = 1 # Toroidal mode number diff --git a/src/KineticForces/Torque.jl b/src/KineticForces/Torque.jl index 82de0045..17c93133 100644 --- a/src/KineticForces/Torque.jl +++ b/src/KineticForces/Torque.jl @@ -35,7 +35,6 @@ function tpsi!(tpsi_var::Ref{ComplexF64}, psi::Float64, n::Int, l::Int, op_wmats::Union{Nothing,Array{ComplexF64,3}}=nothing, rex_override::Union{Nothing,Float64}=nothing, imx_override::Union{Nothing,Float64}=nothing, - ion_fraction::Float64=1.0, atol_xlmda::Float64=1e-9, rtol_xlmda::Float64=1e-6) # Enforce bounds @@ -154,12 +153,11 @@ function tpsi!(tpsi_var::Ref{ComplexF64}, psi::Float64, n::Int, l::Int, dT_s_dpsi = kinetic_profiles.Te_deriv(psi) nu_s = kinetic_profiles.nue_spline(psi) else - # ion_fraction scales this species' resonant density out of the total main-ion n_i - # (e.g. 0.5 for each of D/T). It cancels in wdian/wdiat (T·(dn/dψ)/n) and does NOT - # scale nu_s, which is the full-composition collisionality. - n_s = ion_fraction * kinetic_profiles.ni_spline(psi) + # ni_spline / nui_spline carry THIS species' resonant density and its + # full-composition collisionality (per-species views built by the loader). + n_s = kinetic_profiles.ni_spline(psi) T_s = kinetic_profiles.Ti_spline(psi) - dn_s_dpsi = ion_fraction * kinetic_profiles.ni_deriv(psi) + dn_s_dpsi = kinetic_profiles.ni_deriv(psi) dT_s_dpsi = kinetic_profiles.Ti_deriv(psi) nu_s = kinetic_profiles.nui_spline(psi) end From 23c82d6ebcea3fc867eeddcb7363a76a07f88d9b Mon Sep 17 00:00:00 2001 From: logan-nc Date: Wed, 29 Jul 2026 12:23:09 -0400 Subject: [PATCH 05/17] KineticForces - NEW FEATURE - internal multi-species NTV loop + summed output Wires the ion_species list end to end for the fgar NTV (post-PE diagnostic): when ctrl.ion_species is non-empty, build per-species profiles (build_species_profiles), run compute_torque_all_methods! once per species (+ electron if enabled) sharing the perturbed field and the full-composition Zeff, and sum via combine_species_states (exact scalar totals; dT/dpsi summed on the union grid via the per-species torque_profile interpolants; T(psi) re-integrated). Output: per-species groups kinetic_forces_ion_z_m / kinetic_forces_electron, plus the summed total in the standard kinetic_forces group (so existing analysis reads the total). Empty ion_species => unchanged single path. write_to_hdf5! gains a group_name kwarg. The self-consistent kinetic-DCON path stays single-ion. Co-Authored-By: Claude Opus 4.8 (1M context) Claude-Session: https://claude.ai/code/session_01SsLUydP2gJbE1tGoaHDiS1 --- src/GeneralizedPerturbedEquilibrium.jl | 46 ++++++++++++++++++++++---- src/KineticForces/Compute.jl | 34 +++++++++++++++++++ src/KineticForces/Output.jl | 4 +-- 3 files changed, 76 insertions(+), 8 deletions(-) diff --git a/src/GeneralizedPerturbedEquilibrium.jl b/src/GeneralizedPerturbedEquilibrium.jl index edee8303..cfb08057 100755 --- a/src/GeneralizedPerturbedEquilibrium.jl +++ b/src/GeneralizedPerturbedEquilibrium.jl @@ -605,12 +605,46 @@ function main_from_inputs( kf_intr = KineticForces.KineticForcesInternal(equil; verbose=kf_ctrl.verbose) KineticForces.set_perturbation_data!(kf_intr, pe_state, intr, equil, metric) - kf_state = KineticForces.KineticForcesState() - KineticForces.compute_torque_all_methods!(kf_state, kf_intr, kf_ctrl, equil, kinetic_profiles) - - if kf_ctrl.write_outputs_to_HDF5 - h5open(joinpath(intr.dir_path, kf_ctrl.HDF5_filename), "cw") do h5file - KineticForces.write_to_hdf5!(h5file, kf_state) + if isempty(kf_ctrl.ion_species) + # Single main ion (from zi/mi) — unchanged path. + kf_state = KineticForces.KineticForcesState() + KineticForces.compute_torque_all_methods!(kf_state, kf_intr, kf_ctrl, equil, kinetic_profiles) + if kf_ctrl.write_outputs_to_HDF5 + h5open(joinpath(intr.dir_path, kf_ctrl.HDF5_filename), "cw") do h5file + KineticForces.write_to_hdf5!(h5file, kf_state) + end + end + else + # Multi main-ion: NTV per species on one shared full-composition Zeff, then summed. + kfile = joinpath(intr.dir_path, kf_ctrl.kinetic_file) + species_profs = Equilibrium.build_species_profiles(kfile, kf_ctrl.ion_species; + zimp=kf_ctrl.zimp, mimp=kf_ctrl.mimp) + states = KineticForces.KineticForcesState[] + labels = String[] + for (si, sp) in enumerate(kf_ctrl.ion_species) + sctrl = deepcopy(kf_ctrl) + sctrl.zi = sp.z; sctrl.mi = sp.m; sctrl.electron = false + sctrl.ion_species = KineticForces.IonSpecies[] + st = KineticForces.KineticForcesState() + KineticForces.compute_torque_all_methods!(st, kf_intr, sctrl, equil, species_profs[si]) + push!(states, st); push!(labels, "ion$(si)_z$(sp.z)_m$(sp.m)") + end + if kf_ctrl.electron + ectrl = deepcopy(kf_ctrl) + ectrl.electron = true; ectrl.ion_species = KineticForces.IonSpecies[] + st = KineticForces.KineticForcesState() + KineticForces.compute_torque_all_methods!(st, kf_intr, ectrl, equil, species_profs[1]) + push!(states, st); push!(labels, "electron") + end + @info "Multi-ion NTV: summed $(length(kf_ctrl.ion_species)) ion species$(kf_ctrl.electron ? " + electrons" : "")" + kf_state = KineticForces.combine_species_states(states) + if kf_ctrl.write_outputs_to_HDF5 + h5open(joinpath(intr.dir_path, kf_ctrl.HDF5_filename), "cw") do h5file + for (st, lbl) in zip(states, labels) + KineticForces.write_to_hdf5!(h5file, st; group_name="kinetic_forces_$lbl") + end + KineticForces.write_to_hdf5!(h5file, kf_state) # "kinetic_forces" = total + end end end end diff --git a/src/KineticForces/Compute.jl b/src/KineticForces/Compute.jl index 4175a62b..c1e768fa 100644 --- a/src/KineticForces/Compute.jl +++ b/src/KineticForces/Compute.jl @@ -213,6 +213,40 @@ end # High-level orchestration # ============================================================================ +""" + combine_species_states(states) -> KineticForcesState + +Sum per-species [`KineticForcesState`](@ref) results into a single total (τ = Σ_s τ_s). Per +method: the scalar `total_torque`/`total_energy` are summed exactly; the dT/dψ profile is summed +by evaluating each species' `torque_profile` interpolant on the sorted union of the species ψ +grids, and the cumulative T(ψ) is re-integrated (trapezoid) from that summed profile. All species +share the same ψ-integration range (`ctrl.psilims`), so the grids differ only in adaptive nodes. +""" +function combine_species_states(states::AbstractVector{KineticForcesState}) + combined = KineticForcesState() + isempty(states) && return combined + method_names = collect(keys(first(states).method_results)) + for mname in method_names + results = [s.method_results[mname] for s in states if haskey(s.method_results, mname)] + isempty(results) && continue + grid = sort(unique(reduce(vcat, (r.psi_grid for r in results); init=Float64[]))) + _eval(r, ψ) = r.torque_profile === nothing ? zero(ComplexF64) : ComplexF64(r.torque_profile(ψ)) + dtdpsi = ComplexF64[sum(_eval(r, ψ) for r in results) for ψ in grid] + tcum = zeros(ComplexF64, length(grid)) + for j in 2:length(grid) + tcum[j] = tcum[j-1] + 0.5 * (dtdpsi[j] + dtdpsi[j-1]) * (grid[j] - grid[j-1]) + end + combined.method_results[mname] = MethodResult(; + method=mname, nn=first(results).nn, + total_torque=sum(r.total_torque for r in results), + total_energy=sum(r.total_energy for r in results), + psi_grid=grid, dtdpsi=dtdpsi, t_cumulative=tcum, + psi_nsteps=sum(r.psi_nsteps for r in results)) + end + combined.completed = true + return combined +end + """ compute_torque_all_methods!(state::KineticForcesState, intr::KineticForcesInternal, ctrl::KineticForcesControl, equil, kinetic_profiles) diff --git a/src/KineticForces/Output.jl b/src/KineticForces/Output.jl index d1fef4e0..6fadf73a 100644 --- a/src/KineticForces/Output.jl +++ b/src/KineticForces/Output.jl @@ -15,8 +15,8 @@ Write all KineticForces results to the "kinetic_forces" group in gpec.h5. - `h5file::HDF5.File`: Open HDF5 file handle - `state::KineticForcesState`: Accumulated computation results """ -function write_to_hdf5!(h5file::HDF5.File, state::KineticForcesState) - g = create_group(h5file, "kinetic_forces") +function write_to_hdf5!(h5file::HDF5.File, state::KineticForcesState; group_name::AbstractString="kinetic_forces") + g = create_group(h5file, group_name) for (method_name, result) in state.method_results mg = create_group(g, method_name) From de8e2661f1a02d8052bb7ef1c3b521d1fe694d14 Mon Sep 17 00:00:00 2001 From: logan-nc Date: Wed, 29 Jul 2026 12:28:04 -0400 Subject: [PATCH 06/17] KineticForces - BUG FIX - accept Vector{Any} for ion_species TOML conversion TOML.jl parses [[KineticForces.ion_species]] as Vector{Any} (not Vector{Dict}), so the AbstractVector{<:AbstractDict} convert did not match. Broaden to AbstractVector with per-element construction (pass IonSpecies through). Co-Authored-By: Claude Opus 4.8 (1M context) Claude-Session: https://claude.ai/code/session_01SsLUydP2gJbE1tGoaHDiS1 --- src/KineticForces/KineticForcesStructs.jl | 5 ++++- 1 file changed, 4 insertions(+), 1 deletion(-) diff --git a/src/KineticForces/KineticForcesStructs.jl b/src/KineticForces/KineticForcesStructs.jl index d2c2e835..7b990762 100644 --- a/src/KineticForces/KineticForcesStructs.jl +++ b/src/KineticForces/KineticForcesStructs.jl @@ -72,7 +72,10 @@ end IonSpecies(; z::Integer, m::Integer, fraction=NaN, density="") = IonSpecies(Int(z), Int(m), Float64(fraction), String(density)) IonSpecies(d::AbstractDict) = IonSpecies(; (Symbol(k) => v for (k, v) in d)...) -Base.convert(::Type{Vector{IonSpecies}}, v::AbstractVector{<:AbstractDict}) = IonSpecies.(v) +# TOML.jl parses [[KineticForces.ion_species]] as a Vector{Any} of Dicts, so accept any +# AbstractVector and build each element (pass IonSpecies through unchanged). +Base.convert(::Type{Vector{IonSpecies}}, v::AbstractVector) = + IonSpecies[x isa IonSpecies ? x : IonSpecies(x) for x in v] """ KineticForcesControl From 2b1be4c648fb11eb633b4718e7723786d4e6e892 Mon Sep 17 00:00:00 2001 From: logan-nc Date: Wed, 29 Jul 2026 12:41:01 -0400 Subject: [PATCH 07/17] KineticForces - NEW FEATURE - unified multi-species sum (KF quadrature + kinetic matrices) incl impurity Both NTV paths now loop over and sum the SAME full species set. resolve_ntv_species (was build_species_profiles) returns descriptors (z, m, profiles, electron, label) for: the main ions, the neutrality-closing impurity (zimp/mimp, density = the quasineutrality n_imp) as its own resonant species, and (optionally) electrons -- all on one shared full-composition Zeff/zpitch/lnLambda. Resolved once in main() and passed to: - the KineticForces quadrature path (per-species compute + combine_species_states), and - compute_calculated_kinetic_matrices (self-consistent DCON Mode A): the kinetic W/torque matrices are additive over species, so the threaded psi-loop is wrapped in a species loop that accumulates (+=) each species' block-diagonal contribution (species=nothing => single-species, unchanged). Empty ion_species => single-species from zi/mi (bit-for-bit unchanged). Verified: D-T-e config resolves to 4 summed species (D, T, C impurity, e). Co-Authored-By: Claude Opus 4.8 (1M context) Claude-Session: https://claude.ai/code/session_01SsLUydP2gJbE1tGoaHDiS1 --- src/Equilibrium/KineticProfiles.jl | 61 +++++++++------ src/GeneralizedPerturbedEquilibrium.jl | 35 ++++----- .../CalculatedKineticMatrices.jl | 74 +++++++++++-------- 3 files changed, 98 insertions(+), 72 deletions(-) diff --git a/src/Equilibrium/KineticProfiles.jl b/src/Equilibrium/KineticProfiles.jl index 129b8232..a0bdf6c5 100644 --- a/src/Equilibrium/KineticProfiles.jl +++ b/src/Equilibrium/KineticProfiles.jl @@ -161,22 +161,24 @@ function multi_ion_composition(zs::AbstractVector, ns::AbstractVector, ne::Real; end """ - build_species_profiles(kinetic_file, ion_species; zimp, mimp, ...) -> Vector{KineticProfileSplines} - -Build one [`KineticProfileSplines`](@ref) per main-ion species for a multi-ion NTV run. Each -returned view carries **this species' resonant density** `n_s` (as `ni_spline`) and its -**full-composition collision frequency** `ν_s` (as `nui_spline`), while `ne/Te/ωE/Zeff/lnΛ` -(and the electron `nue`) are shared across all species. Downstream `tpsi!` uses each view -unchanged; the driver loops over the returned vector and sums the torques. - -Species density is resolved from `IonSpecies`: `fraction` ⇒ `n_s = fraction · n_i` (the kinetic -file's `n_i` = total main-ion density); explicit per-species profiles (`density`) are not yet -wired here. Composition (Zeff, zpitch) comes from [`multi_ion_composition`]; the per-species -collisionality is `ν_s = (zpitch/3.5e17)·z_s²·n_main·lnΛ / (√m_s·(T_i)^{3/2})`. Reduces to the -single-ion `load_kinetic_profiles` result for one z=1 species with `fraction=1`. + resolve_ntv_species(kinetic_file, ion_species; electron, zimp, mimp, ...) -> Vector{NamedTuple} + +Resolve the **full NTV species set** for a multi-ion run: the main ions, the neutrality-closing +impurity, and (if `electron`) the electrons. Returns a vector of descriptors +`(z, m, profiles::KineticProfileSplines, electron::Bool, label)` — the single set that BOTH NTV +paths (the KineticForces quadrature and the self-consistent kinetic-matrix build) loop over and sum. + +Each descriptor's `profiles` view carries **that species' resonant density** `n_s` (as `ni_spline`) +and its **full-composition collision frequency** `ν_s` (as `nui_spline`); `ne/Te/ωE/Zeff/lnΛ` (and +the electron `nue`) are shared. Main-ion density is resolved from `IonSpecies`: `fraction` ⇒ +`n_s = fraction · n_i` (the kinetic file's `n_i` = total main-ion density); explicit per-species +profiles (`density`) are not yet wired. Composition (Zeff, zpitch, `n_imp`) comes from +[`multi_ion_composition`]; per-species `ν_s = (zpitch/3.5e17)·z_s²·n_main·lnΛ / (√m_s·(T_i)^{3/2})`. +The impurity (`zimp`,`mimp`, density = the quasineutrality `n_imp`) is included as its own resonant +species. Reduces to `load_kinetic_profiles` for one z=1 main ion with `fraction=1`. """ -function build_species_profiles(kinetic_file::AbstractString, ion_species::AbstractVector; - zimp::Integer=6, mimp::Integer=12, +function resolve_ntv_species(kinetic_file::AbstractString, ion_species::AbstractVector; + electron::Bool=false, zimp::Integer=6, mimp::Integer=12, density_factor::Float64=1.0, temperature_factor::Float64=1.0, ExB_rotation_factor::Float64=1.0, toroidal_rotation_factor::Float64=1.0) @@ -213,21 +215,34 @@ function build_species_profiles(kinetic_file::AbstractString, ion_species::Abstr # Shared composition + collisionality (natural-log Coulomb log), per grid point. npts = length(psi_reg) - zeff = zeros(npts); zpitch = zeros(npts); loglam = zeros(npts); n_main = zeros(npts); nue = zeros(npts) + zeff = zeros(npts); zpitch = zeros(npts); loglam = zeros(npts) + n_main = zeros(npts); n_imp = zeros(npts); nue = zeros(npts) for i in 1:npts - z, zp, nm, _ = multi_ion_composition(zs, [ns[si][i] for si in eachindex(ion_species)], ne[i]; zimp=zimp, mimp=mimp) - zeff[i] = z; zpitch[i] = zp; n_main[i] = nm + z, zp, nm, nimp = multi_ion_composition(zs, [ns[si][i] for si in eachindex(ion_species)], ne[i]; zimp=zimp, mimp=mimp) + zeff[i] = z; zpitch[i] = zp; n_main[i] = nm; n_imp[i] = max(nimp, 0.0) loglam[i] = (ne[i] > 0 && Te[i] > 0) ? 17.3 - 0.5 * log(ne[i] / 1.0e20) + 1.5 * log(Te[i] / 1.602e-16) : 0.0 nue[i] = Te[i] > 0 ? (zp / 3.5e17) * ne[i] * loglam[i] / (sqrt(me / mp) * (Te[i] / 1.602e-16)^1.5) : 0.0 end - # One KineticProfileSplines per species: ni = n_s, nui = ν_s, everything else shared. - profs = Vector{KineticProfileSplines}(undef, length(ion_species)) + # Per-species collision frequency ν_s (shared zpitch/n_main/Zeff/lnΛ; species z², m, T_i) and a + # KineticProfileSplines view carrying that species' density as `ni_spline`. + _nu(zsp, msp) = [Ti[i] > 0 ? (zpitch[i] / 3.5e17) * zsp^2 * n_main[i] * loglam[i] / (sqrt(Float64(msp)) * (Ti[i] / 1.602e-16)^1.5) : 0.0 for i in 1:npts] + _view(dens, nu) = KineticProfileSplines(psi_reg, dens, ne, Ti, Te, omegaE, loglam, nu, nue, zeff) + + # Full NTV species set: main ions, the neutrality-closing impurity, and (optionally) electrons. + species = NamedTuple[] for (si, s) in enumerate(ion_species) - nu_s = [Ti[i] > 0 ? (zpitch[i] / 3.5e17) * s.z^2 * n_main[i] * loglam[i] / (sqrt(Float64(s.m)) * (Ti[i] / 1.602e-16)^1.5) : 0.0 for i in 1:npts] - profs[si] = KineticProfileSplines(psi_reg, ns[si], ne, Ti, Te, omegaE, loglam, nu_s, nue, zeff) + push!(species, (z=Int(s.z), m=Int(s.m), profiles=_view(ns[si], _nu(s.z, s.m)), + electron=false, label="ion$(si)_z$(s.z)_m$(s.m)")) + end + if any(>(0), n_imp) + push!(species, (z=Int(zimp), m=Int(mimp), profiles=_view(n_imp, _nu(zimp, mimp)), + electron=false, label="impurity_z$(zimp)_m$(mimp)")) + end + if electron + push!(species, (z=-1, m=1, profiles=_view(ns[1], nue), electron=true, label="electron")) end - return profs + return species end """ diff --git a/src/GeneralizedPerturbedEquilibrium.jl b/src/GeneralizedPerturbedEquilibrium.jl index cfb08057..104a00e0 100755 --- a/src/GeneralizedPerturbedEquilibrium.jl +++ b/src/GeneralizedPerturbedEquilibrium.jl @@ -218,6 +218,10 @@ function main_from_inputs( KineticForces.KineticForcesControl() kinetic_profiles = nothing + # Multi-species NTV set (main ions + neutrality impurity + electrons), resolved once and + # summed by BOTH the self-consistent kinetic-matrix build and the KineticForces quadrature. + # `nothing` ⇒ single-species from kf_ctrl.zi/mi (unchanged behaviour). + kf_species = nothing needs_kinetic_profiles = haskey(inputs, "KineticForces") || (ctrl.kinetic_factor > 0 && ctrl.kinetic_source == "calculated") if needs_kinetic_profiles @@ -229,6 +233,12 @@ function main_from_inputs( density_factor=kf_ctrl.density_factor, temperature_factor=kf_ctrl.temperature_factor, ExB_rotation_factor=kf_ctrl.ExB_rotation_factor, toroidal_rotation_factor=kf_ctrl.toroidal_rotation_factor, chi1=2π * equil.psio) + if !isempty(kf_ctrl.ion_species) + kf_species = Equilibrium.resolve_ntv_species(kinetic_file, kf_ctrl.ion_species; + electron=kf_ctrl.electron, zimp=kf_ctrl.zimp, mimp=kf_ctrl.mimp, + density_factor=kf_ctrl.density_factor, temperature_factor=kf_ctrl.temperature_factor, + ExB_rotation_factor=kf_ctrl.ExB_rotation_factor, toroidal_rotation_factor=kf_ctrl.toroidal_rotation_factor) + end end # Two-pass auto grid: measure the pass-1 equilibrium's curvature (profiles, geometry, @@ -418,7 +428,7 @@ function main_from_inputs( 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) + kf_ctrl=kf_ctrl, kinetic_profiles=kinetic_profiles, species=kf_species) make_kinetic_matrix(ctrl, equil, ffit, intr, metric; calculated_source=calculated_cb) @@ -615,28 +625,19 @@ function main_from_inputs( end end else - # Multi main-ion: NTV per species on one shared full-composition Zeff, then summed. - kfile = joinpath(intr.dir_path, kf_ctrl.kinetic_file) - species_profs = Equilibrium.build_species_profiles(kfile, kf_ctrl.ion_species; - zimp=kf_ctrl.zimp, mimp=kf_ctrl.mimp) + # Multi-species: summed over main ions + neutrality impurity + electrons + # (kf_species, resolved once above and shared with the kinetic-matrix build). states = KineticForces.KineticForcesState[] labels = String[] - for (si, sp) in enumerate(kf_ctrl.ion_species) + for sp in kf_species sctrl = deepcopy(kf_ctrl) - sctrl.zi = sp.z; sctrl.mi = sp.m; sctrl.electron = false + sctrl.zi = sp.z; sctrl.mi = sp.m; sctrl.electron = sp.electron sctrl.ion_species = KineticForces.IonSpecies[] st = KineticForces.KineticForcesState() - KineticForces.compute_torque_all_methods!(st, kf_intr, sctrl, equil, species_profs[si]) - push!(states, st); push!(labels, "ion$(si)_z$(sp.z)_m$(sp.m)") - end - if kf_ctrl.electron - ectrl = deepcopy(kf_ctrl) - ectrl.electron = true; ectrl.ion_species = KineticForces.IonSpecies[] - st = KineticForces.KineticForcesState() - KineticForces.compute_torque_all_methods!(st, kf_intr, ectrl, equil, species_profs[1]) - push!(states, st); push!(labels, "electron") + KineticForces.compute_torque_all_methods!(st, kf_intr, sctrl, equil, sp.profiles) + push!(states, st); push!(labels, sp.label) end - @info "Multi-ion NTV: summed $(length(kf_ctrl.ion_species)) ion species$(kf_ctrl.electron ? " + electrons" : "")" + @info "Multi-species NTV: summed $(length(kf_species)) species ($(join((s.label for s in kf_species), ", ")))" kf_state = KineticForces.combine_species_states(states) if kf_ctrl.write_outputs_to_HDF5 h5open(joinpath(intr.dir_path, kf_ctrl.HDF5_filename), "cw") do h5file diff --git a/src/KineticForces/CalculatedKineticMatrices.jl b/src/KineticForces/CalculatedKineticMatrices.jl index d36fead9..76b66008 100644 --- a/src/KineticForces/CalculatedKineticMatrices.jl +++ b/src/KineticForces/CalculatedKineticMatrices.jl @@ -56,6 +56,7 @@ function compute_calculated_kinetic_matrices( ffit; kf_ctrl::KineticForcesControl = KineticForcesControl(), kinetic_profiles::Equilibrium.KineticProfileSplines, + species::Union{Nothing,AbstractVector} = nothing, ) xs = metric.xs mpsi = length(xs) @@ -101,38 +102,47 @@ function compute_calculated_kinetic_matrices( 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] - block_w = thread_block_w[tid] - block_t = thread_block_t[tid] - psi = xs[ipsi] - for in_idx in 1:npert - n = ffs_intr.nlow + in_idx - 1 - fill!(full_w, 0) - fill!(full_t, 0) - for ell in -nl:nl - fill!(block_w, 0) - fill!(block_t, 0) - compute_kinetic_matrices_at_psi!( - block_w, block_t, psi, n, ell, - kf_ctrl.zi, kf_ctrl.mi, kf_ctrl.wdfac, kf_ctrl.divxfac, - kf_ctrl.electron, equil, intr_t, kinetic_profiles; - nutype=kf_ctrl.nutype, f0type=kf_ctrl.f0type, nufac=kf_ctrl.nufac, - atol_xlmda=kf_ctrl.atol_xlmda, rtol_xlmda=kf_ctrl.rtol_xlmda - ) - full_w .+= block_w - full_t .+= block_t - end - - # Place the n-block on the diagonal of the full np×np matrix and flatten. - 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] + # Species set summed into the kinetic matrices: the resolved multi-species set + # (main ions + neutrality impurity + electrons) when provided, else a single species + # from kf_ctrl. The kinetic W/torque matrices are additive over species, so we + # accumulate (+=) each species' block-diagonal contribution. + splist = species === nothing ? + [(z=kf_ctrl.zi, m=kf_ctrl.mi, electron=kf_ctrl.electron, profiles=kinetic_profiles)] : species + + for sp in splist + 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] + block_w = thread_block_w[tid] + block_t = thread_block_t[tid] + psi = xs[ipsi] + for in_idx in 1:npert + n = ffs_intr.nlow + in_idx - 1 + fill!(full_w, 0) + fill!(full_t, 0) + for ell in -nl:nl + fill!(block_w, 0) + fill!(block_t, 0) + compute_kinetic_matrices_at_psi!( + block_w, block_t, psi, n, ell, + sp.z, sp.m, kf_ctrl.wdfac, kf_ctrl.divxfac, + sp.electron, equil, intr_t, sp.profiles; + nutype=kf_ctrl.nutype, f0type=kf_ctrl.f0type, nufac=kf_ctrl.nufac, + atol_xlmda=kf_ctrl.atol_xlmda, rtol_xlmda=kf_ctrl.rtol_xlmda + ) + full_w .+= block_w + full_t .+= block_t + end + + # Place the n-block on the diagonal of the full np×np matrix and accumulate. + 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] + end end end end From 2fe72e495ea16b24723edcdf89b0f532ccbb9a40 Mon Sep 17 00:00:00 2001 From: logan-nc Date: Wed, 29 Jul 2026 00:22:42 -0400 Subject: [PATCH 08/17] ForceFreeStates - BUG FIX - Pass direction arg to IntegrationChunk in kinetic EL crossing cross_kinetic_singular_surf! built its placeholder IntegrationChunk with a stale 4-arg signature (psi_start, psi_end, needs_crossing, ising), but the struct's @kwdef positional constructor requires the direction field added later for bidirectional parallel FM. The 4-arg call throws MethodError. This path is only reached on the serial Euler-Lagrange shooting solver (use_parallel=false) with a kinetic singular surface (kinetic_factor>0), which previously errored earlier in the Riccati branch, so the defect stayed latent. Now passes direction=1, matching the sibling placeholder call sites in cross_ideal_singular_surf! and the Riccati crossing. Co-Authored-By: Claude Opus 4.8 (1M context) Claude-Session: https://claude.ai/code/session_01SsLUydP2gJbE1tGoaHDiS1 --- src/ForceFreeStates/EulerLagrange.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/ForceFreeStates/EulerLagrange.jl b/src/ForceFreeStates/EulerLagrange.jl index 9f93e60a..0eab1935 100644 --- a/src/ForceFreeStates/EulerLagrange.jl +++ b/src/ForceFreeStates/EulerLagrange.jl @@ -676,7 +676,7 @@ function cross_kinetic_singular_surf!( ksurf = intr.kinsing[ising] dpsi = ksurf.psifac - odet.psifac - params = (ctrl, equil, ffit, intr, odet, IntegrationChunk(0.0, 0.0, false, 0)) + params = (ctrl, equil, ffit, intr, odet, IntegrationChunk(0.0, 0.0, false, 0, 1)) du1 = zeros(ComplexF64, intr.numpert_total, intr.numpert_total, 2) du2 = zeros(ComplexF64, intr.numpert_total, intr.numpert_total, 2) From ee20627fa7e7bacf370b704302c3678a46b89782 Mon Sep 17 00:00:00 2001 From: logan-nc Date: Wed, 29 Jul 2026 13:31:10 -0400 Subject: [PATCH 09/17] KineticForces - NEW FEATURE - explicit per-species density profiles (HDF5) Completes the multi-ion IO: KineticProfileData gains a species_densities dict populated by _read_kinetic_h5 from any non-standard HDF5 datasets (e.g. "n_D", "n_T"). resolve_ntv_species' `density=` branch resolves a species from its named profile (resampled), the same downstream path as the `fraction` shorthand -- so measured, differently-shaped per-species profiles are supported via the HDF5 container while ASCII+fraction stays the simple case. One-of-{fraction,density} per species is validated; a missing named profile errors with the available names. Verified: n_D=0.6*Ni, n_T=0.4*Ni explicit datasets resolve to the correct per-species densities, impurity closes quasineutrality from them. Co-Authored-By: Claude Opus 4.8 (1M context) Claude-Session: https://claude.ai/code/session_01SsLUydP2gJbE1tGoaHDiS1 --- src/Equilibrium/KineticProfiles.jl | 33 ++++++++++++++++++++++++------ 1 file changed, 27 insertions(+), 6 deletions(-) diff --git a/src/Equilibrium/KineticProfiles.jl b/src/Equilibrium/KineticProfiles.jl index a0bdf6c5..eccd7989 100644 --- a/src/Equilibrium/KineticProfiles.jl +++ b/src/Equilibrium/KineticProfiles.jl @@ -41,14 +41,18 @@ struct KineticProfileData omega_tor::Union{Nothing,Vector{Float64}} chi_e::Union{Nothing,Vector{Float64}} chi_phi::Union{Nothing,Vector{Float64}} + # Named per-species density profiles (e.g. "n_D", "n_T") for explicit multi-ion input, + # read from non-standard datasets in the HDF5 kinetic file. `nothing` for ASCII / no extras. + species_densities::Union{Nothing,Dict{String,Vector{Float64}}} provenance::String end function KineticProfileData(; psi, n_i=nothing, n_e=nothing, T_i=nothing, T_e=nothing, - omega_E=nothing, omega_tor=nothing, chi_e=nothing, chi_phi=nothing, provenance="") + omega_E=nothing, omega_tor=nothing, chi_e=nothing, chi_phi=nothing, + species_densities=nothing, provenance="") _v(x) = x === nothing ? nothing : Float64.(collect(x)) return KineticProfileData(Float64.(collect(psi)), _v(n_i), _v(n_e), _v(T_i), _v(T_e), - _v(omega_E), _v(omega_tor), _v(chi_e), _v(chi_phi), String(provenance)) + _v(omega_E), _v(omega_tor), _v(chi_e), _v(chi_phi), species_densities, String(provenance)) end const _KINETIC_H5_EXTS = (".h5", ".hdf5", ".he5") @@ -89,10 +93,20 @@ function _read_kinetic_h5(path::AbstractString; group::AbstractString="/") prov = "" ats = attributes(g) haskey(ats, "provenance") && (prov = string(read(ats["provenance"]))) + # Any dataset outside the standard schema is treated as a named per-species density + # profile (e.g. "n_D", "n_T") for explicit multi-ion input. + standard = ("psi", "n_i", "n_e", "T_i", "T_e", "omega_E", "omega_tor", "chi_e", "chi_phi") + extras = Dict{String,Vector{Float64}}() + for k in keys(g) + k in standard && continue + v = read(g[k]) + v isa AbstractArray && (extras[k] = Float64.(vec(v))) + end return KineticProfileData(; psi=Float64.(vec(read(g["psi"]))), n_i=rd("n_i"), n_e=rd("n_e"), T_i=rd("T_i"), T_e=rd("T_e"), omega_E=rd("omega_E"), omega_tor=rd("omega_tor"), - chi_e=rd("chi_e"), chi_phi=rd("chi_phi"), provenance=prov) + chi_e=rd("chi_e"), chi_phi=rd("chi_phi"), + species_densities=(isempty(extras) ? nothing : extras), provenance=prov) end end @@ -202,15 +216,22 @@ function resolve_ntv_species(kinetic_file::AbstractString, ion_species::Abstract Te = _cubic_resample(psi_in, Te_in, psi_reg) .* eV_to_J omegaE = _cubic_resample(psi_in, omE_in, psi_reg) - # Resolve each species' resonant density (fraction path). + # Resolve each species' resonant density: `fraction` ⇒ share of the total n_i; `density` ⇒ + # an explicit named profile from the HDF5 kinetic file (resampled to the working grid). zs = [Int(s.z) for s in ion_species] ns = Vector{Vector{Float64}}(undef, length(ion_species)) for (si, s) in enumerate(ion_species) has_frac = !isnan(s.fraction) has_dens = !isempty(s.density) (has_frac ⊻ has_dens) || error("ion_species[$si]: specify exactly one of `fraction` or `density`") - has_dens && error("ion_species[$si]: explicit per-species density `$(s.density)` not yet wired; use `fraction`") - ns[si] = s.fraction .* ni_total + if has_dens + (data.species_densities !== nothing && haskey(data.species_densities, s.density)) || + error("ion_species[$si]: density profile `$(s.density)` not found in kinetic file " * + "(available: $(data.species_densities === nothing ? "none — ASCII files carry only n_i" : join(keys(data.species_densities), ", ")))") + ns[si] = _cubic_resample(psi_in, data.species_densities[s.density], psi_reg) + else + ns[si] = s.fraction .* ni_total + end end # Shared composition + collisionality (natural-log Coulomb log), per grid point. From 976a362bf90351105e0477373c336ab21f7b0944 Mon Sep 17 00:00:00 2001 From: logan-nc Date: Wed, 29 Jul 2026 13:34:04 -0400 Subject: [PATCH 10/17] KineticForces - TEST - multi-ion NTV unit tests (25) Covers: multi_ion_composition (exact single-ion reduction for z=1,2; 50/50 D-T Zeff; quasineutrality closes the impurity), IonSpecies + TOML Vector{Any}->Vector{IonSpecies} conversion, resolve_ntv_species (single-species reduces to load_kinetic_profiles; D-T-e set = {D,T,impurity,electron}; nu ~ 1/sqrt(m); one-of-{fraction,density} validation), and explicit per-species HDF5 density profiles. Co-Authored-By: Claude Opus 4.8 (1M context) Claude-Session: https://claude.ai/code/session_01SsLUydP2gJbE1tGoaHDiS1 --- test/runtests.jl | 1 + test/runtests_multiion.jl | 87 +++++++++++++++++++++++++++++++++++++++ 2 files changed, 88 insertions(+) create mode 100644 test/runtests_multiion.jl diff --git a/test/runtests.jl b/test/runtests.jl index 9bd7b0ed..221ab8de 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -33,6 +33,7 @@ else include("./runtests_innerlayer.jl") include("./runtests_tj_analytic.jl") include("./runtests_kinetic.jl") + include("./runtests_multiion.jl") include("./runtests_fullruns.jl") include("./runtests_coils.jl") include("./runtests_imas.jl") diff --git a/test/runtests_multiion.jl b/test/runtests_multiion.jl new file mode 100644 index 00000000..0a052198 --- /dev/null +++ b/test/runtests_multiion.jl @@ -0,0 +1,87 @@ +@testset "Multi-ion NTV" begin + EQ = GeneralizedPerturbedEquilibrium.Equilibrium + KF = GeneralizedPerturbedEquilibrium.KineticForces + + @testset "multi_ion_composition" begin + ne = 1.24e20; ni = 1.13e20; zimp = 6; mimp = 12 + # Reduces exactly to the single-ion formula Zeff = zimp - (ni/ne)*zi*(zimp-zi). + for zi in (1, 2) + zeff_old = zimp - (ni / ne) * zi * (zimp - zi) + zeff, zpitch, nmain, nimp = EQ.multi_ion_composition([zi], [ni], ne; zimp=zimp, mimp=mimp) + @test zeff ≈ zeff_old rtol = 1e-12 + @test nmain ≈ ni + end + # 50/50 D-T with total ni == single ion at total ni; quasineutrality closes the impurity. + zdt, _, nmdt, nimpdt = EQ.multi_ion_composition([1, 1], [ni / 2, ni / 2], ne; zimp=zimp, mimp=mimp) + zsingle, _, _, _ = EQ.multi_ion_composition([1], [ni], ne; zimp=zimp, mimp=mimp) + @test zdt ≈ zsingle rtol = 1e-12 + @test nmdt ≈ ni + @test ne ≈ ni + zimp * nimpdt rtol = 1e-12 # ne = Σ z_s n_s + z_imp n_imp + end + + @testset "IonSpecies + TOML conversion" begin + # TOML.jl parses [[KineticForces.ion_species]] as a Vector{Any} of Dicts. + v = Any[Dict("z" => 1, "m" => 2, "fraction" => 0.5), Dict("z" => 1, "m" => 3, "fraction" => 0.5)] + c = KF.KineticForcesControl(; ion_species=v) + @test length(c.ion_species) == 2 + @test (c.ion_species[1].z, c.ion_species[1].m, c.ion_species[1].fraction) == (1, 2, 0.5) + @test isempty(KF.KineticForcesControl().ion_species) # default single-ion + end + + @testset "resolve_ntv_species" begin + # Synthetic ASCII kinetic profile: psi, n_i(total), n_e, T_i, T_e, omega_E. + psi = collect(0.0:0.05:1.0) + Ni = @. 1.0e20 * (1 - 0.5psi); ne = @. 1.1e20 * (1 - 0.5psi) + Ti = @. 2000.0 * (1 - 0.5psi); Te = @. 2500.0 * (1 - 0.5psi); wE = @. 1.0e4 * (1 - psi) + f = tempname() * ".txt" + open(f, "w") do io + println(io, "# psi ni ne Ti Te wexb") + for i in eachindex(psi) + println(io, join((psi[i], Ni[i], ne[i], Ti[i], Te[i], wE[i]), " ")) + end + end + + # Single species [z=1,m=2,fraction=1] reproduces load_kinetic_profiles. + single = EQ.load_kinetic_profiles(f; zi=1, zimp=6, mi=2, mimp=12) + sp1 = EQ.resolve_ntv_species(f, [KF.IonSpecies(z=1, m=2, fraction=1.0)]; electron=false, zimp=6, mimp=12) + main = sp1[1] + for ψ in (0.3, 0.6) + @test main.profiles.ni_spline(ψ) ≈ single.ni_spline(ψ) rtol = 1e-8 + @test main.profiles.nui_spline(ψ) ≈ single.nui_spline(ψ) rtol = 1e-8 + @test main.profiles.zeff_spline(ψ) ≈ single.zeff_spline(ψ) rtol = 1e-8 + end + + # D-T-e set = {D, T, impurity, electron}; ν scales as 1/√m at fixed z, T. + sp = EQ.resolve_ntv_species(f, [KF.IonSpecies(z=1, m=2, fraction=0.5), KF.IonSpecies(z=1, m=3, fraction=0.5)]; + electron=true, zimp=6, mimp=12) + @test length(sp) == 4 + @test count(s -> !s.electron && s.z == 1, sp) == 2 # D, T + @test any(s -> s.z == 6 && !s.electron, sp) # impurity + @test any(s -> s.electron, sp) # electron + iD = findfirst(s -> s.m == 2, sp); iT = findfirst(s -> s.m == 3, sp) + @test sp[iD].profiles.nui_spline(0.5) / sp[iT].profiles.nui_spline(0.5) ≈ sqrt(3 / 2) rtol = 1e-6 + + # Exactly one of {fraction, density} required. + @test_throws ErrorException EQ.resolve_ntv_species(f, [KF.IonSpecies(z=1, m=2)]; zimp=6, mimp=12) + rm(f) + end + + @testset "explicit per-species profiles (HDF5)" begin + psi = collect(0.0:0.05:1.0) + Ni = @. 1.0e20 * (1 - 0.5psi); ne = @. 1.1e20 * (1 - 0.5psi) + Ti = @. 2000.0 * (1 - 0.5psi); Te = @. 2500.0 * (1 - 0.5psi); wE = @. 1.0e4 * (1 - psi) + h5 = tempname() * ".h5" + HDF5.h5open(h5, "w") do file + file["psi"] = psi; file["n_e"] = ne; file["T_i"] = Ti; file["T_e"] = Te; file["omega_E"] = wE + file["n_i"] = Ni; file["n_D"] = 0.6 .* Ni; file["n_T"] = 0.4 .* Ni # unequal explicit shapes + end + sp = EQ.resolve_ntv_species(h5, [KF.IonSpecies(z=1, m=2, density="n_D"), KF.IonSpecies(z=1, m=3, density="n_T")]; + electron=false, zimp=6, mimp=12) + iD = findfirst(s -> s.m == 2, sp); iT = findfirst(s -> s.m == 3, sp) + @test sp[iD].profiles.ni_spline(0.5) ≈ 0.6 * 1.0e20 * (1 - 0.25) rtol = 1e-6 + @test sp[iT].profiles.ni_spline(0.5) ≈ 0.4 * 1.0e20 * (1 - 0.25) rtol = 1e-6 + # A missing named profile errors clearly. + @test_throws ErrorException EQ.resolve_ntv_species(h5, [KF.IonSpecies(z=1, m=2, density="n_He")]; zimp=6, mimp=12) + rm(h5) + end +end From 34d4895be25e02cf04fd51533b47ce4c43ce109f Mon Sep 17 00:00:00 2001 From: logan-nc Date: Wed, 29 Jul 2026 13:47:02 -0400 Subject: [PATCH 11/17] KineticForces - TEST - multi-ion NTV regression case (solovev_kinetic_multiion) New example examples/Solovev_kinetic_multiion_example (Solovev main ion modelled as 50/50 D-T) and regression case solovev_kinetic_multiion tracking the summed total plus each per-species NTV contribution (D, T, electron; pure fixture so no impurity), pinning the multi-species resolve_ntv_species + combine_species_states loop end to end. Verified: total 1.318e-4 = D+T+e sum. Co-Authored-By: Claude Opus 4.8 (1M context) Claude-Session: https://claude.ai/code/session_01SsLUydP2gJbE1tGoaHDiS1 --- .../forcing.dat | 5 + .../gpec.toml | 109 ++++++++++++++++++ .../kinetic.dat | 23 ++++ .../run_example.jl | 4 + .../cases/solovev_kinetic_multiion.toml | 67 +++++++++++ 5 files changed, 208 insertions(+) create mode 100644 examples/Solovev_kinetic_multiion_example/forcing.dat create mode 100644 examples/Solovev_kinetic_multiion_example/gpec.toml create mode 100644 examples/Solovev_kinetic_multiion_example/kinetic.dat create mode 100644 examples/Solovev_kinetic_multiion_example/run_example.jl create mode 100644 regression-harness/cases/solovev_kinetic_multiion.toml diff --git a/examples/Solovev_kinetic_multiion_example/forcing.dat b/examples/Solovev_kinetic_multiion_example/forcing.dat new file mode 100644 index 00000000..42eb9470 --- /dev/null +++ b/examples/Solovev_kinetic_multiion_example/forcing.dat @@ -0,0 +1,5 @@ +# Forcing data for perturbed equilibrium calculations +# normalization: normal_field_T +# Format: n m amplitude_real amplitude_imag +# Single mode test case: n=1, m=2, amplitude=1e-4 T +1 2 1e-4 0.0 diff --git a/examples/Solovev_kinetic_multiion_example/gpec.toml b/examples/Solovev_kinetic_multiion_example/gpec.toml new file mode 100644 index 00000000..b869e04e --- /dev/null +++ b/examples/Solovev_kinetic_multiion_example/gpec.toml @@ -0,0 +1,109 @@ +# Solovev analytical equilibrium — n=1 ideal stability + perturbed equilibrium + NTV torque. +# Extends the Solovev ideal example with a [KineticForces] section, so after the plasma +# response is computed the neoclassical toroidal viscosity torque is integrated over ψ as a +# post-processing diagnostic (adaptive quadrature paneled at the rational surfaces). +# The kinetic profiles (kinetic.dat) are constructed to be consistent with this equilibrium's +# pressure; the equilibrium is generated analytically from the embedded [SOL_INPUT] section. + +[Equilibrium] +eq_type = "sol" # Type of the input 2D equilibrium file +jac_type = "pest" # Coordinate system (hamada, pest, boozer, equal_arc, park, other) +grid_type = "ldp" # Radial grid packing type +psilow = 1e-4 # Lower limit of normalized poloidal flux +psihigh = 0.9995 # Upper limit of normalized poloidal flux +mpsi = 128 # Number of radial grid points (0 = auto-compute from psi_accuracy) +mtheta = 256 # Number of poloidal grid points +newq0 = 0 # Override for on-axis safety factor (0 = use input value) +etol = 1e-7 # Error tolerance for equilibrium solver +force_termination = false # Terminate after equilibrium setup (skip stability calculations) + +[Wall] +# Close conformal wall is required to stabilize this Solovev fixture's n=1 external kink: +# with nowall, et[1] = -6.8 (strongly unstable); with this wall, et[1] = +0.24 (barely stable). +# The plasma is near marginal stability, so the BVP Δ' matrix values are pathological +# (dpm magnitudes ~ 10¹¹, |Im/Re| ≫ 1). This fixture's role is integration-pipeline +# smoke testing + et[1] regression, NOT BVP Δ' regression — DIIID-like is the canonical +# Δ'-matrix fixture (stable et[1] = +1.6, clean BVP Δ'). +shape = "conformal" # Wall shape (nowall, conformal, elliptical, dee, mod_dee, filepath) +a = 0.2415 # Distance from plasma (conformal) or shape parameter +aw = 0.05 # Half-thickness parameter for Dee-shaped walls +bw = 1.5 # Elongation parameter for wall shapes +cw = 0 # Offset of wall center from major radius +dw = 0.5 # Triangularity parameter for wall shapes +tw = 0.05 # Sharpness of wall corners (try 0.05 as initial value) +equal_arc_wall = true # Equal arc length distribution of nodes on wall + +[ForcingTerms] +forcing_data_file = "forcing.dat" # Manual mode table (n, m, amplitude_real, amplitude_imag) +forcing_data_format = "ascii" # Format: "ascii", "hdf5", or "coil" (Biot-Savart from 3D wires) + +[PerturbedEquilibrium] +fixed_boundary = false # Use fixed boundary conditions +output_eigenmodes = true # Output eigenmode fields as b-fields +compute_response = true # Compute plasma response to forcing +compute_singular_coupling = true # Compute singular layer coupling metrics +verbose = true # Enable verbose logging +write_outputs_to_HDF5 = true # Write perturbed equilibrium outputs to HDF5 +reg_spot = 0.05 # Regularization width for singular surfaces (0 = disabled) + +[ForceFreeStates] +local_stability_flag = true # Perform local stability analysis (Mercier and ballooning) across the ψ profile +mat_flag = true # Construct coefficient matrices for diagnostic purposes +ode_flag = true # Integrate ODEs for stability of the internal long-wavelength mode (must be true for GPEC) +vac_flag = true # Compute plasma, vacuum, and total energies for free-boundary modes + +psiedge = 0.99 # Edge dW(ψ) diagnostic scan band [psiedge, psilim]; set ≥ psilim to disable +qlow = 1.02 # Integration initiated at q determined by min(q0, qlow) +qhigh = 1e3 # Integration terminated at q limit determined by min(qa, qhigh) +sing_start = 0 # Start integration at the sing_start'th rational from the axis (psilow) + +nn_low = 1 # Smallest toroidal mode number to include +nn_high = 1 # Largest toroidal mode number to include +delta_mlow = 8 # Expands lower bound of Fourier harmonics +delta_mhigh = 8 # Expands upper bound of Fourier harmonics +mthvac = 960 # Number of points used in splines over poloidal angle at the plasma-vacuum interface + +kinetic_source = "fixed" # Kinetic matrix source: "fixed" test matrices, or "calculated" from the kinetic NTV model +kinetic_factor = 0.0 # Scaling of kinetic matrices (0 = ideal path; >0 enables kinetic mode) +eulerlagrange_tolerance = 1e-7 # Relative tolerance for ODE integration of Euler-Lagrange equations +singfac_min = 1e-4 # Fractional distance from rational q at which ideal jump enforced +ucrit = 1e3 # Column-norm threshold that triggers solution renormalization +force_wv_symmetry = true # Enforce symmetry of the vacuum energy matrix +save_interval = 3 # Save every Nth ODE step (1=all). Always saves near rational surfaces. + +# Δ' BVP + parallel integration (see ForceFreeStatesControl docstring for details) +use_parallel = true # Run parallel FM-propagator BVP path (unlocks singular/delta_prime_matrix) +parallel_threads = 2 # BVP thread cap (1 = serial/bit-deterministic; 2 ≈ +20% speedup; ≥3 saturates) +populate_dense_xi = true # Append serial-EL pass so dense ξ is stored — REQUIRED with a [PerturbedEquilibrium] section +set_psilim_via_dmlim = false # FALSE for limited/analytical equilibria — rationals sparse, dmlim would chop too much edge +dmlim = 0.2 # Truncate integration at (last_rational_q + dmlim)/n (used when set_psilim_via_dmlim = true) +# Solovev analytic equilibrium parameters (eq_type = "sol"); see SolovevConfig in src/Equilibrium. +[SOL_INPUT] +mr = 128 # Number of radial grid zones +mz = 128 # Number of axial grid zones +ma = 128 # Number of flux grid zones +e = 1.6 # Elongation +a = 0.33 # Minor radius +r0 = 1.0 # Major radius +q0 = 1.9 # Safety factor at the magnetic axis (O-point) +p0fac = 1 # Scale on-axis pressure (β changes; Φ, q fixed) +b0fac = 1 # Scale toroidal field at constant β (Bt changes; shape, β fixed) +f0fac = 1 # Scale toroidal field at constant pressure (β, q change; Φ, p, Bp fixed) + +[KineticForces] +# Multi-ion NTV: the analytic Solovev main ion is modelled as a 50/50 D-T mix (n_i = total), +# so the summed NTV runs over D + T + the neutrality-closing impurity + electrons. +kinetic_file = "kinetic.dat" +electron = true +zimp = 6 +mimp = 12 + +[[KineticForces.ion_species]] +z = 1 +m = 2 +fraction = 0.5 + +[[KineticForces.ion_species]] +z = 1 +m = 3 +fraction = 0.5 diff --git a/examples/Solovev_kinetic_multiion_example/kinetic.dat b/examples/Solovev_kinetic_multiion_example/kinetic.dat new file mode 100644 index 00000000..96fd9fa3 --- /dev/null +++ b/examples/Solovev_kinetic_multiion_example/kinetic.dat @@ -0,0 +1,23 @@ +# Kinetic profiles consistent with this Solovev equilibrium's pressure. +# Constraint: P(psi_n) = n_e*T_e + n_i*T_i must equal the equilibrium pressure +# P(psi_n) = P0*(1 - psi_n), P0 = 4.273e4 Pa (linear, zero at the edge). +# Construction (quasineutral, n_e=n_i, T_e=T_i): +# T(psi_n) = 100 + 2600*(1 - psi_n) eV (core 2700 eV, edge 100 eV) +# n(psi_n) = P0*(1 - psi_n) / (2*e*T) (core ~4.9e19; product n0*T0 = P0/2e fixed) +# Density carries the edge falloff to a small floor at psi_n=1 so T stays finite +# (avoids T->0 in the collisionality / drift-frequency formulas); the residual edge +# pressure from the floor is ~1e-6 of P0. omega_E is a free (non-pressure) profile. +psi_n n_i_m3 n_e_m3 T_i_eV T_e_eV omega_E_rads +0.00 4.939e19 4.939e19 2700.0 2700.0 1.0e4 +0.10 4.919e19 4.919e19 2440.0 2440.0 9.0e3 +0.20 4.894e19 4.894e19 2180.0 2180.0 8.0e3 +0.30 4.862e19 4.862e19 1920.0 1920.0 7.0e3 +0.40 4.820e19 4.820e19 1660.0 1660.0 6.0e3 +0.50 4.763e19 4.763e19 1400.0 1400.0 5.0e3 +0.60 4.679e19 4.679e19 1140.0 1140.0 4.0e3 +0.70 4.546e19 4.546e19 880.0 880.0 3.0e3 +0.80 4.302e19 4.302e19 620.0 620.0 2.0e3 +0.90 3.705e19 3.705e19 360.0 360.0 1.0e3 +0.95 2.899e19 2.899e19 230.0 230.0 5.0e2 +0.99 1.058e19 1.058e19 126.0 126.0 1.0e2 +1.00 2.000e18 2.000e18 100.0 100.0 0.0 diff --git a/examples/Solovev_kinetic_multiion_example/run_example.jl b/examples/Solovev_kinetic_multiion_example/run_example.jl new file mode 100644 index 00000000..fa0a9e53 --- /dev/null +++ b/examples/Solovev_kinetic_multiion_example/run_example.jl @@ -0,0 +1,4 @@ +using Pkg; +Pkg.activate(joinpath(@__DIR__, "../..")) +using GeneralizedPerturbedEquilibrium +GeneralizedPerturbedEquilibrium.main([dirname(@__FILE__)]) diff --git a/regression-harness/cases/solovev_kinetic_multiion.toml b/regression-harness/cases/solovev_kinetic_multiion.toml new file mode 100644 index 00000000..df97e674 --- /dev/null +++ b/regression-harness/cases/solovev_kinetic_multiion.toml @@ -0,0 +1,67 @@ +# Regression case: Solovev n=1 + PE + MULTI-ION NTV torque. The Solovev main ion is modelled as a +# 50/50 D-T mix, so the KineticForces quadrature sums over D + T + electrons (this fixture is a pure +# plasma, n_e = n_i, so no impurity species is present). Tracks the summed total AND each per-species +# contribution, pinning the multi-species loop (resolve_ntv_species + combine_species_states) and the +# full-composition Zeff. The impurity-species path is covered by the unit tests + the TC-24 case. +[case] +name = "solovev_kinetic_multiion" +description = "Solovev n=1, ideal + PE + multi-ion NTV (D + T + electron summed)" +example_dir = "examples/Solovev_kinetic_multiion_example" + +# Summed total NTV over all species — the headline multi-ion quantity. +[quantities.ntv_total] +h5path = "kinetic_forces/fgar/total_torque" +type = "real_scalar" +extract = "value" +label = "NTV total (D+T+imp+e) [N·m]" +noise_threshold = 1e-8 +order = 10 + +[quantities.ntv_D] +h5path = "kinetic_forces_ion1_z1_m2/fgar/total_torque" +type = "real_scalar" +extract = "value" +label = "NTV Deuterium [N·m]" +noise_threshold = 1e-8 +order = 11 + +[quantities.ntv_T] +h5path = "kinetic_forces_ion2_z1_m3/fgar/total_torque" +type = "real_scalar" +extract = "value" +label = "NTV Tritium [N·m]" +noise_threshold = 1e-8 +order = 12 + +[quantities.ntv_electron] +h5path = "kinetic_forces_electron/fgar/total_torque" +type = "real_scalar" +extract = "value" +label = "NTV electron [N·m]" +noise_threshold = 1e-8 +order = 14 + +# Stability anchors — confirm the upstream FFS/PE stages feeding the NTV diagnostic. +[quantities.et_real] +h5path = "FreeBoundaryStability/eigenmode_energies" +type = "complex_vector" +extract = "real_first" +label = "total energy Re(et[1])" +noise_threshold = 1e-10 +order = 20 + +[quantities.msing] +h5path = "singular/msing" +type = "int_scalar" +extract = "value" +label = "# singular surfaces" +noise_threshold = 0 +order = 21 + +[quantities.runtime] +h5path = "" +type = "runtime" +extract = "value" +label = "Runtime (s)" +noise_threshold = 0.0 +order = 90 From bc1c3956c964278214fb10b92691487e7f02dc0c Mon Sep 17 00:00:00 2001 From: logan-nc Date: Wed, 29 Jul 2026 13:55:10 -0400 Subject: [PATCH 12/17] KineticForces - CLEANUP - annotate zpitch domain + electron view density (per physics review) Per fortran-physics-reviewer sign-off of the multi-species NTV: (1) flag that zpitch is a main-ion momentum-restoring closure applied approximately to the impurity/electron test species; (2) the electron descriptor now carries ne (not ns[1]) in the unused ni_spline slot for clarity. No numerical change. Reviewer verdict: collision-freq z^2/n_main form CORRECT (n_main is the faithful Fortran choice, not ne*Zeff); the self-consistent dW +0.066->-0.10 shift is additive-channel physics, not a summation bug (ideal F/K/G added once downstream, kw/kt purely kinetic). Co-Authored-By: Claude Opus 4.8 (1M context) Claude-Session: https://claude.ai/code/session_01SsLUydP2gJbE1tGoaHDiS1 --- src/Equilibrium/KineticProfiles.jl | 8 ++++++-- 1 file changed, 6 insertions(+), 2 deletions(-) diff --git a/src/Equilibrium/KineticProfiles.jl b/src/Equilibrium/KineticProfiles.jl index eccd7989..c40abde8 100644 --- a/src/Equilibrium/KineticProfiles.jl +++ b/src/Equilibrium/KineticProfiles.jl @@ -246,7 +246,10 @@ function resolve_ntv_species(kinetic_file::AbstractString, ion_species::Abstract end # Per-species collision frequency ν_s (shared zpitch/n_main/Zeff/lnΛ; species z², m, T_i) and a - # KineticProfileSplines view carrying that species' density as `ni_spline`. + # KineticProfileSplines view carrying that species' density as `ni_spline`. NB: `zpitch` is + # strictly a main-ion momentum-restoring closure; applying it to the impurity/electron test + # species is an approximation beyond the single-ion theory (acceptable since the impurity NTV + # scales with its small resonant density n_imp). _nu(zsp, msp) = [Ti[i] > 0 ? (zpitch[i] / 3.5e17) * zsp^2 * n_main[i] * loglam[i] / (sqrt(Float64(msp)) * (Ti[i] / 1.602e-16)^1.5) : 0.0 for i in 1:npts] _view(dens, nu) = KineticProfileSplines(psi_reg, dens, ne, Ti, Te, omegaE, loglam, nu, nue, zeff) @@ -261,7 +264,8 @@ function resolve_ntv_species(kinetic_file::AbstractString, ion_species::Abstract electron=false, label="impurity_z$(zimp)_m$(mimp)")) end if electron - push!(species, (z=-1, m=1, profiles=_view(ns[1], nue), electron=true, label="electron")) + # The electron view's `ni_spline` is unused (the electron path reads `ne_spline`); carry ne. + push!(species, (z=-1, m=1, profiles=_view(ne, nue), electron=true, label="electron")) end return species end From 851b4d257db30fd35c4eaf3a88a3db427080ae49 Mon Sep 17 00:00:00 2001 From: logan-nc Date: Wed, 29 Jul 2026 15:13:36 -0400 Subject: [PATCH 13/17] KineticForces - BUG FIX - combine_species_states summed profile was all zeros MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit The multi-species combined dT/dψ profile (written to the "kinetic_forces" total group) was identically zero: combine summed via each MethodResult.torque_profile interpolant, but that field is not populated by the per-species compute (only the dtdpsi array is). total_torque was correct (summed scalar), so the bug was invisible until the combined mid-radius profile was integrated. Fix: linear-interpolate each species' (psi_grid, dtdpsi) arrays onto the union grid and sum (zero outside a species' range). Unit test added (overlap sum + nonzero + total). Per-species output groups were always correct. Co-Authored-By: Claude Opus 4.8 (1M context) Claude-Session: https://claude.ai/code/session_01SsLUydP2gJbE1tGoaHDiS1 --- src/KineticForces/Compute.jl | 20 ++++++++++++++++++-- test/runtests_multiion.jl | 17 +++++++++++++++++ 2 files changed, 35 insertions(+), 2 deletions(-) diff --git a/src/KineticForces/Compute.jl b/src/KineticForces/Compute.jl index c1e768fa..869f16a8 100644 --- a/src/KineticForces/Compute.jl +++ b/src/KineticForces/Compute.jl @@ -230,8 +230,24 @@ function combine_species_states(states::AbstractVector{KineticForcesState}) results = [s.method_results[mname] for s in states if haskey(s.method_results, mname)] isempty(results) && continue grid = sort(unique(reduce(vcat, (r.psi_grid for r in results); init=Float64[]))) - _eval(r, ψ) = r.torque_profile === nothing ? zero(ComplexF64) : ComplexF64(r.torque_profile(ψ)) - dtdpsi = ComplexF64[sum(_eval(r, ψ) for r in results) for ψ in grid] + # Sum each species' dT/dψ onto the union grid by linear interpolation of its own + # (psi_grid, dtdpsi) arrays (torque_profile interpolants are not populated here); zero + # outside a species' ψ range. + dtdpsi = zeros(ComplexF64, length(grid)) + for r in results + length(r.psi_grid) >= 2 || continue + o = sortperm(r.psi_grid); xs = r.psi_grid[o]; ys = r.dtdpsi[o] + for (k, ψ) in enumerate(grid) + (ψ < xs[1] || ψ > xs[end]) && continue + j = searchsortedlast(xs, ψ) + if j == length(xs) + dtdpsi[k] += ys[end] + else + t = (ψ - xs[j]) / (xs[j+1] - xs[j]) + dtdpsi[k] += (1 - t) * ys[j] + t * ys[j+1] + end + end + end tcum = zeros(ComplexF64, length(grid)) for j in 2:length(grid) tcum[j] = tcum[j-1] + 0.5 * (dtdpsi[j] + dtdpsi[j-1]) * (grid[j] - grid[j-1]) diff --git a/test/runtests_multiion.jl b/test/runtests_multiion.jl index 0a052198..5adb97f1 100644 --- a/test/runtests_multiion.jl +++ b/test/runtests_multiion.jl @@ -84,4 +84,21 @@ @test_throws ErrorException EQ.resolve_ntv_species(h5, [KF.IonSpecies(z=1, m=2, density="n_He")]; zimp=6, mimp=12) rm(h5) end + + @testset "combine_species_states" begin + # Combined profile = Σ species dT/dψ interpolated onto the union grid; total = Σ total. + mk(psi, dt) = (st = KF.KineticForcesState(); + st.method_results["fgar"] = KF.MethodResult(; method="fgar", nn=1, + total_torque=ComplexF64(sum(dt)), psi_grid=collect(Float64, psi), + dtdpsi=ComplexF64.(dt), t_cumulative=zeros(ComplexF64, length(dt))); + st) + s1 = mk(0.0:0.1:1.0, fill(1.0, 11)) + s2 = mk(0.05:0.1:0.95, fill(2.0, 10)) # different (offset) grid + c = KF.combine_species_states([s1, s2]) + r = c.method_results["fgar"] + @test real(r.total_torque) ≈ 31.0 # 11·1 + 10·2 + @test maximum(abs.(r.dtdpsi)) > 0 # not the all-zeros regression + overlap = [real(d) for (p, d) in zip(r.psi_grid, r.dtdpsi) if 0.1 <= p <= 0.9] + @test all(x -> isapprox(x, 3.0; atol=1e-9), overlap) # 1 + 2 on the overlap + end end From 2edbb120b689494c1fd1884a7e789834a09d7e8f Mon Sep 17 00:00:00 2001 From: logan-nc Date: Wed, 29 Jul 2026 18:17:12 -0400 Subject: [PATCH 14/17] KineticForces - CLEANUP - concise Coulomb-log annotation (CLAUDE.md comment rule) Collapse the 3-line lnLambda comment to one self-contained line: drop the 'not log10' bug-history note and the stale Fortran file:line, per the keep-comments-concise rule. Co-Authored-By: Claude Opus 4.8 (1M context) Claude-Session: https://claude.ai/code/session_01SsLUydP2gJbE1tGoaHDiS1 --- src/Equilibrium/KineticProfiles.jl | 4 +--- 1 file changed, 1 insertion(+), 3 deletions(-) diff --git a/src/Equilibrium/KineticProfiles.jl b/src/Equilibrium/KineticProfiles.jl index c40abde8..67006ed7 100644 --- a/src/Equilibrium/KineticProfiles.jl +++ b/src/Equilibrium/KineticProfiles.jl @@ -409,9 +409,7 @@ function load_kinetic_profiles(kinetic_file::AbstractString; z = n_e > 0 ? zimp - (n_i / n_e) * zi * (zimp - zi) : Float64(zimp) zpitch = 1.0 + (1.0 + mimp) / (2.0 * mimp) * zimp * (z - 1.0) / (zimp - z) - # Coulomb logarithm — natural log, matching Fortran PENTRC inputs.f90:238 (the - # coefficients 17.3/-0.5/+1.5 are calibrated for ln, not log10; n_e in units of 1e20 m^-3, - # T_e in units of 1 keV = 1.602e-16 J). + # Coulomb logarithm (NRL formulary form), natural log; n_e in 1e20 m^-3, T_e in keV. ll = 17.3 - 0.5 * log(n_e / 1.0e20) + 1.5 * log(T_e / 1.602e-16) loglam[i] = ll From 66cc9deefdbe3152e32a34a2538a1ea619e3ddf2 Mon Sep 17 00:00:00 2001 From: logan-nc Date: Thu, 30 Jul 2026 00:12:44 -0400 Subject: [PATCH 15/17] ForceFreeStates - CLEANUP - pass real ising to placeholder IntegrationChunk (PR #336 review) Match the sibling call sites (cross_ideal_singular_surf!, Riccati crossing): pass the in-scope ising instead of 0. Functionally inert (sing_der! discards the chunk); keeps the cherry-picked IntegrationChunk fix in sync with the amended bugfix branch. Co-Authored-By: Claude Opus 4.8 (1M context) Claude-Session: https://claude.ai/code/session_01SsLUydP2gJbE1tGoaHDiS1 --- src/ForceFreeStates/EulerLagrange.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/ForceFreeStates/EulerLagrange.jl b/src/ForceFreeStates/EulerLagrange.jl index 0eab1935..1f8fe3bc 100644 --- a/src/ForceFreeStates/EulerLagrange.jl +++ b/src/ForceFreeStates/EulerLagrange.jl @@ -676,7 +676,7 @@ function cross_kinetic_singular_surf!( ksurf = intr.kinsing[ising] dpsi = ksurf.psifac - odet.psifac - params = (ctrl, equil, ffit, intr, odet, IntegrationChunk(0.0, 0.0, false, 0, 1)) + params = (ctrl, equil, ffit, intr, odet, IntegrationChunk(0.0, 0.0, false, ising, 1)) du1 = zeros(ComplexF64, intr.numpert_total, intr.numpert_total, 2) du2 = zeros(ComplexF64, intr.numpert_total, intr.numpert_total, 2) From 5f9f40eeac6c0623ec914a608ff04dceeb0c0a56 Mon Sep 17 00:00:00 2001 From: logan-nc Date: Fri, 31 Jul 2026 12:05:19 -0400 Subject: [PATCH 16/17] =?UTF-8?q?KineticForces=20-=20IMPROVEMENT=20-=20Add?= =?UTF-8?q?ress=20multi-ion=20NTV=20review:=20test-particle=20z=C2=B2,=20d?= =?UTF-8?q?ensity=20hard-stop,=20typed=20species?= MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit - ν_s carries the test-particle z² (Krook deflection, Logan & Park 2013 Eq. 6); the single-ion path now equals the multi-ion limit (identical at z=1). - Hard error on any negative density (bad input or cubic-resample overshoot) instead of silently clamping the neutrality-closing impurity. - resolve_ntv_species returns concrete-eltype ResolvedNTVSpecies structs, so the self-consistent kinetic-matrix species loop is type-stable. - Share physical constants and the Coulomb-log helper across load_kinetic_profiles and resolve_ntv_species; drop duplicate local literals. - Fix combine_species_states docstring (linear interp of dtdpsi, dropped fields); move the multi-ion contract into the KineticForcesControl / KineticProfileData docstrings; trim inline field comments. - Make runtests_multiion self-contained (using HDF5). Co-Authored-By: Claude Opus 4.8 (1M context) Claude-Session: https://claude.ai/code/session_01SsLUydP2gJbE1tGoaHDiS1 --- src/Equilibrium/KineticProfiles.jl | 180 +++++++++++------- .../CalculatedKineticMatrices.jl | 49 ++--- src/KineticForces/Compute.jl | 58 +++--- src/KineticForces/KineticForcesStructs.jl | 102 +++++----- test/runtests_multiion.jl | 53 ++++-- 5 files changed, 264 insertions(+), 178 deletions(-) diff --git a/src/Equilibrium/KineticProfiles.jl b/src/Equilibrium/KineticProfiles.jl index 67006ed7..359d26c0 100644 --- a/src/Equilibrium/KineticProfiles.jl +++ b/src/Equilibrium/KineticProfiles.jl @@ -12,6 +12,42 @@ shimming. using DelimitedFiles using HDF5 +# Physical constants and collisionality normalizations shared by `load_kinetic_profiles` and +# `resolve_ntv_species` (Krook collision frequency, Logan & Park 2013 Eq. 6). +const _EV_J = 1.602e-19 # eV → J +const _KEV_J = 1.602e-16 # 1 keV in J (temperature normalization in the collision frequency) +const _NU_PREFAC = 3.5e17 # Krook collision-frequency prefactor +const _MP = 1.672_614e-27 # proton mass [kg] +const _ME = 9.109_1e-31 # electron mass [kg] + +# NRL-formulary Coulomb logarithm (natural log); n_e in m⁻³, T_e in J. +_coulomb_log(ne, Te) = 17.3 - 0.5 * log(ne / 1.0e20) + 1.5 * log(Te / _KEV_J) + +# Hard stop on an unphysical negative density (bad input or a cubic-resample overshoot). +function _assert_nonneg_density(kinetic_file, label, arr, psi_reg) + i = findfirst(<(0), arr) + i === nothing || + error("kinetic file '$kinetic_file': negative $label = $(arr[i]) at ψ_n = $(psi_reg[i]) — unphysical") +end + +""" + ResolvedNTVSpecies{P} + +One fully-resolved species in the NTV sum — a main ion, the neutrality-closing impurity, or the +electrons. Fields: charge `z`, mass `m` (proton masses), `electron` flag, `label`, and `profiles` +(a `KineticProfileSplines` view carrying this species' density as `ni_spline` and its +full-composition collision frequency as `nui_spline`). Produced by `resolve_ntv_species`; both NTV +paths — the KineticForces ψ-quadrature and the self-consistent kinetic-matrix build — iterate this +and sum τ = Σ_s τ_s. +""" +struct ResolvedNTVSpecies{P} + z::Int + m::Int + electron::Bool + label::String + profiles::P +end + """ KineticProfileData @@ -21,15 +57,17 @@ present; every profile field is `nothing` when the source omits it, so each consumer validates only the fields it needs. Units: densities m⁻³, temperatures eV, frequencies rad/s, diffusivities m²/s. -| field | meaning | -|-------------|--------------------------------------| -| `n_i` | main-ion density | -| `n_e` | electron density | -| `T_i`/`T_e` | ion / electron temperature | -| `omega_E` | ExB rotation | -| `omega_tor` | toroidal rotation (optional) | -| `chi_e` | perpendicular heat diffusivity χ⊥ | -| `chi_phi` | toroidal momentum diffusivity χ_φ | +| field | meaning | +|:------------------- |:----------------------------------------------------------------------------------------------------------------- | +| `n_i` | main-ion density | +| `n_e` | electron density | +| `T_i`/`T_e` | ion / electron temperature | +| `omega_E` | ExB rotation | +| `omega_tor` | toroidal rotation (optional) | +| `chi_e` | perpendicular heat diffusivity χ⊥ | +| `chi_phi` | toroidal momentum diffusivity χ_φ | +| `species_densities` | named per-species densities (e.g. `"n_D"`, `"n_T"`) for explicit multi-ion input; `nothing` for ASCII / no extras | +| `provenance` | short string recording the source file/format | """ struct KineticProfileData psi::Vector{Float64} @@ -41,8 +79,6 @@ struct KineticProfileData omega_tor::Union{Nothing,Vector{Float64}} chi_e::Union{Nothing,Vector{Float64}} chi_phi::Union{Nothing,Vector{Float64}} - # Named per-species density profiles (e.g. "n_D", "n_T") for explicit multi-ion input, - # read from non-standard datasets in the HDF5 kinetic file. `nothing` for ASCII / no extras. species_densities::Union{Nothing,Dict{String,Vector{Float64}}} provenance::String end @@ -76,15 +112,19 @@ function read_kinetic_file(path::AbstractString; group::AbstractString="/") return ext in _KINETIC_H5_EXTS ? _read_kinetic_h5(path; group=group) : _read_kinetic_ascii(path) end -"""Read a legacy 6-column ASCII kinetic table into a `KineticProfileData`.""" +""" +Read a legacy 6-column ASCII kinetic table into a `KineticProfileData`. +""" function _read_kinetic_ascii(path::AbstractString) psi, ni, ne, Ti, Te, omegaE = _read_kinetic_table(path) return KineticProfileData(; psi=psi, n_i=ni, n_e=ne, T_i=Ti, T_e=Te, omega_E=omegaE, provenance="ASCII 6-column: $(basename(path))") end -"""Read the GPEC HDF5 kinetic schema. Only `psi` is required; other datasets are -optional and surfaced as-is.""" +""" +Read the GPEC HDF5 kinetic schema. Only `psi` is required; other datasets are +optional and surfaced as-is. +""" function _read_kinetic_h5(path::AbstractString; group::AbstractString="/") h5open(path, "r") do f g = group == "/" ? f : f[group] @@ -165,7 +205,7 @@ this is `ν_s = (zpitch/3.5e17)·z_s²·n_main·lnΛ / (√m_s · (T_s)^{3/2})` `Zeff`, `lnΛ`; the test species contributes only its own `z_s²`, `m_s`, `T_s`. """ function multi_ion_composition(zs::AbstractVector, ns::AbstractVector, ne::Real; - zimp::Real, mimp::Real) + zimp::Real, mimp::Real) n_main = sum(ns) charge_density = sum(z * n for (z, n) in zip(zs, ns)) # Σ z_s n_s n_imp = ne > 0 ? (ne - charge_density) / zimp : 0.0 @@ -175,19 +215,19 @@ function multi_ion_composition(zs::AbstractVector, ns::AbstractVector, ne::Real; end """ - resolve_ntv_species(kinetic_file, ion_species; electron, zimp, mimp, ...) -> Vector{NamedTuple} + resolve_ntv_species(kinetic_file, ion_species; electron, zimp, mimp, ...) -> Vector{ResolvedNTVSpecies} Resolve the **full NTV species set** for a multi-ion run: the main ions, the neutrality-closing -impurity, and (if `electron`) the electrons. Returns a vector of descriptors -`(z, m, profiles::KineticProfileSplines, electron::Bool, label)` — the single set that BOTH NTV -paths (the KineticForces quadrature and the self-consistent kinetic-matrix build) loop over and sum. +impurity, and (if `electron`) the electrons. Returns `ResolvedNTVSpecies` descriptors — the +single set that BOTH NTV paths (the KineticForces quadrature and the self-consistent kinetic-matrix +build) loop over and sum. Errors on any negative density (unphysical input / resample overshoot). Each descriptor's `profiles` view carries **that species' resonant density** `n_s` (as `ni_spline`) and its **full-composition collision frequency** `ν_s` (as `nui_spline`); `ne/Te/ωE/Zeff/lnΛ` (and the electron `nue`) are shared. Main-ion density is resolved from `IonSpecies`: `fraction` ⇒ `n_s = fraction · n_i` (the kinetic file's `n_i` = total main-ion density); explicit per-species profiles (`density`) are not yet wired. Composition (Zeff, zpitch, `n_imp`) comes from -[`multi_ion_composition`]; per-species `ν_s = (zpitch/3.5e17)·z_s²·n_main·lnΛ / (√m_s·(T_i)^{3/2})`. +`multi_ion_composition`; per-species `ν_s = (zpitch/3.5e17)·z_s²·n_main·lnΛ / (√m_s·(T_i)^{3/2})`. The impurity (`zimp`,`mimp`, density = the quasineutrality `n_imp`) is included as its own resonant species. Reduces to `load_kinetic_profiles` for one z=1 main ion with `fraction=1`. """ @@ -203,21 +243,24 @@ function resolve_ntv_species(kinetic_file::AbstractString, ion_species::Abstract data = read_kinetic_file(kinetic_file) _need(f, n) = f === nothing ? error("kinetic file '$kinetic_file' missing required '$n'") : f psi_in = data.psi - ne_in = _need(data.n_e, "n_e"); Ti_in = _need(data.T_i, "T_i") - Te_in = _need(data.T_e, "T_e"); omE_in = _need(data.omega_E, "omega_E") + ne_in = _need(data.n_e, "n_e") + Ti_in = _need(data.T_i, "T_i") + Te_in = _need(data.T_e, "T_e") + omE_in = _need(data.omega_E, "omega_E") ni_in = data.n_i === nothing ? copy(ne_in) : data.n_i # TOTAL main-ion density - eV_to_J = 1.602e-19 - mp = 1.672_614e-27; me = 9.109_1e-31 - nkin = 100; psi_reg = collect(0:nkin) ./ nkin + nkin = 100 + psi_reg = collect(0:nkin) ./ nkin ne = _cubic_resample(psi_in, ne_in, psi_reg) ni_total = _cubic_resample(psi_in, ni_in, psi_reg) - Ti = _cubic_resample(psi_in, Ti_in, psi_reg) .* eV_to_J - Te = _cubic_resample(psi_in, Te_in, psi_reg) .* eV_to_J + Ti = _cubic_resample(psi_in, Ti_in, psi_reg) .* _EV_J + Te = _cubic_resample(psi_in, Te_in, psi_reg) .* _EV_J omegaE = _cubic_resample(psi_in, omE_in, psi_reg) + _assert_nonneg_density(kinetic_file, "n_e", ne, psi_reg) + _assert_nonneg_density(kinetic_file, "n_i (total)", ni_total, psi_reg) # Resolve each species' resonant density: `fraction` ⇒ share of the total n_i; `density` ⇒ - # an explicit named profile from the HDF5 kinetic file (resampled to the working grid). + # an explicit named profile from the (HDF5) kinetic file (resampled to the working grid). zs = [Int(s.z) for s in ion_species] ns = Vector{Vector{Float64}}(undef, length(ion_species)) for (si, s) in enumerate(ion_species) @@ -226,48 +269,56 @@ function resolve_ntv_species(kinetic_file::AbstractString, ion_species::Abstract (has_frac ⊻ has_dens) || error("ion_species[$si]: specify exactly one of `fraction` or `density`") if has_dens (data.species_densities !== nothing && haskey(data.species_densities, s.density)) || - error("ion_species[$si]: density profile `$(s.density)` not found in kinetic file " * - "(available: $(data.species_densities === nothing ? "none — ASCII files carry only n_i" : join(keys(data.species_densities), ", ")))") + error( + "ion_species[$si]: density profile `$(s.density)` not found in kinetic file " * + "(available: $(data.species_densities === nothing ? "none — ASCII files carry only n_i" : join(keys(data.species_densities), ", ")))" + ) ns[si] = _cubic_resample(psi_in, data.species_densities[s.density], psi_reg) else ns[si] = s.fraction .* ni_total end + _assert_nonneg_density(kinetic_file, "ion_species[$si] density", ns[si], psi_reg) end # Shared composition + collisionality (natural-log Coulomb log), per grid point. npts = length(psi_reg) - zeff = zeros(npts); zpitch = zeros(npts); loglam = zeros(npts) - n_main = zeros(npts); n_imp = zeros(npts); nue = zeros(npts) + zeff = zeros(npts) + zpitch = zeros(npts) + loglam = zeros(npts) + n_main = zeros(npts) + n_imp = zeros(npts) + nue = zeros(npts) for i in 1:npts z, zp, nm, nimp = multi_ion_composition(zs, [ns[si][i] for si in eachindex(ion_species)], ne[i]; zimp=zimp, mimp=mimp) - zeff[i] = z; zpitch[i] = zp; n_main[i] = nm; n_imp[i] = max(nimp, 0.0) - loglam[i] = (ne[i] > 0 && Te[i] > 0) ? 17.3 - 0.5 * log(ne[i] / 1.0e20) + 1.5 * log(Te[i] / 1.602e-16) : 0.0 - nue[i] = Te[i] > 0 ? (zp / 3.5e17) * ne[i] * loglam[i] / (sqrt(me / mp) * (Te[i] / 1.602e-16)^1.5) : 0.0 + zeff[i] = z + zpitch[i] = zp + n_main[i] = nm + n_imp[i] = nimp + loglam[i] = (ne[i] > 0 && Te[i] > 0) ? _coulomb_log(ne[i], Te[i]) : 0.0 + nue[i] = Te[i] > 0 ? (zp / _NU_PREFAC) * ne[i] * loglam[i] / (sqrt(_ME / _MP) * (Te[i] / _KEV_J)^1.5) : 0.0 end + # A negative impurity density means the main ions over-neutralize (Σ z_s n_s > n_e) — unphysical. + _assert_nonneg_density(kinetic_file, "impurity density (n_e < Σ z_s n_s)", n_imp, psi_reg) - # Per-species collision frequency ν_s (shared zpitch/n_main/Zeff/lnΛ; species z², m, T_i) and a - # KineticProfileSplines view carrying that species' density as `ni_spline`. NB: `zpitch` is - # strictly a main-ion momentum-restoring closure; applying it to the impurity/electron test - # species is an approximation beyond the single-ion theory (acceptable since the impurity NTV - # scales with its small resonant density n_imp). - _nu(zsp, msp) = [Ti[i] > 0 ? (zpitch[i] / 3.5e17) * zsp^2 * n_main[i] * loglam[i] / (sqrt(Float64(msp)) * (Ti[i] / 1.602e-16)^1.5) : 0.0 for i in 1:npts] + # Per-species ν_s: shared zpitch/n_main/Zeff/lnΛ, species-specific z²/m/T_i (Krook deflection + # frequency, test-particle z², Logan & Park 2013 Eq. 6). `zpitch` is a main-ion closure, applied + # approximately to the impurity/electron test species. + _nu(zsp, msp) = [Ti[i] > 0 ? (zpitch[i] / _NU_PREFAC) * zsp^2 * n_main[i] * loglam[i] / (sqrt(Float64(msp)) * (Ti[i] / _KEV_J)^1.5) : 0.0 for i in 1:npts] _view(dens, nu) = KineticProfileSplines(psi_reg, dens, ne, Ti, Te, omegaE, loglam, nu, nue, zeff) # Full NTV species set: main ions, the neutrality-closing impurity, and (optionally) electrons. - species = NamedTuple[] + species = ResolvedNTVSpecies[] for (si, s) in enumerate(ion_species) - push!(species, (z=Int(s.z), m=Int(s.m), profiles=_view(ns[si], _nu(s.z, s.m)), - electron=false, label="ion$(si)_z$(s.z)_m$(s.m)")) + push!(species, ResolvedNTVSpecies(Int(s.z), Int(s.m), false, "ion$(si)_z$(s.z)_m$(s.m)", _view(ns[si], _nu(s.z, s.m)))) end if any(>(0), n_imp) - push!(species, (z=Int(zimp), m=Int(mimp), profiles=_view(n_imp, _nu(zimp, mimp)), - electron=false, label="impurity_z$(zimp)_m$(mimp)")) + push!(species, ResolvedNTVSpecies(Int(zimp), Int(mimp), false, "impurity_z$(zimp)_m$(mimp)", _view(n_imp, _nu(zimp, mimp)))) end if electron # The electron view's `ni_spline` is unused (the electron path reads `ne_spline`); carry ne. - push!(species, (z=-1, m=1, profiles=_view(ne, nue), electron=true, label="electron")) + push!(species, ResolvedNTVSpecies(-1, 1, true, "electron", _view(ne, nue))) end - return species + return identity.(species) # narrow to a concrete-eltype Vector{ResolvedNTVSpecies{...}} end """ @@ -304,11 +355,12 @@ six-column whitespace-separated ASCII table (header rows are filtered out): # Scaling sequence When any of `density_factor`, `temperature_factor`, `toroidal_rotation_factor` differ from 1.0: -1. Build first-pass cubic splines from unscaled profiles (for derivatives) -2. Compute diamagnetic frequencies `wdian`, `wdiat` and total toroidal rotation `wphi` -3. Scale: `wdian_new = temperature_factor * wdian`, `wdiat_new = temperature_factor * wdiat` (`density_factor` cancels in `T*(dn/dψ)/n`) -4. Reform: `omegaE = toroidal_rotation_factor * wphi - wdian_new - wdiat_new` -5. Scale density/temperature arrays: `ni *= density_factor`, `Ti *= temperature_factor`, etc. + + 1. Build first-pass cubic splines from unscaled profiles (for derivatives) + 2. Compute diamagnetic frequencies `wdian`, `wdiat` and total toroidal rotation `wphi` + 3. Scale: `wdian_new = temperature_factor * wdian`, `wdiat_new = temperature_factor * wdiat` (`density_factor` cancels in `T*(dn/dψ)/n`) + 4. Reform: `omegaE = toroidal_rotation_factor * wphi - wdian_new - wdiat_new` + 5. Scale density/temperature arrays: `ni *= density_factor`, `Ti *= temperature_factor`, etc. Then `ExB_rotation_factor` is applied independently: `omegaE *= ExB_rotation_factor`. @@ -334,10 +386,6 @@ function load_kinetic_profiles(kinetic_file::AbstractString; omegaE_input = _need(data.omega_E, "omega_E") ni_input = data.n_i === nothing ? copy(ne_input) : data.n_i - eV_to_J = 1.602e-19 - mp = 1.672_614e-27 - me = 9.109_1e-31 - nkin = 100 psi_reg = collect(0:nkin) ./ nkin @@ -349,8 +397,8 @@ function load_kinetic_profiles(kinetic_file::AbstractString; # denominator). See feedback_kf_kin_profile_linear_interp.md. ni = _cubic_resample(psi_input, ni_input, psi_reg) ne = _cubic_resample(psi_input, ne_input, psi_reg) - Ti = _cubic_resample(psi_input, Ti_input_eV, psi_reg) .* eV_to_J - Te = _cubic_resample(psi_input, Te_input_eV, psi_reg) .* eV_to_J + Ti = _cubic_resample(psi_input, Ti_input_eV, psi_reg) .* _EV_J + Te = _cubic_resample(psi_input, Te_input_eV, psi_reg) .* _EV_J omegaE = _cubic_resample(psi_input, omegaE_input, psi_reg) needs_rotation_reform = density_factor != 1.0 || temperature_factor != 1.0 || toroidal_rotation_factor != 1.0 @@ -370,8 +418,7 @@ function load_kinetic_profiles(kinetic_file::AbstractString; dTi_dpsi = deriv1(Ti_spl) # Compute original wdian, wdiat, wphi and reform omegaE at each grid point - echarge = 1.602e-19 - chrg_ion = zi * echarge + chrg_ion = zi * _EV_J for i in eachindex(omegaE) ψ = psi_reg[i] wdian_i = ni[i] > 0 ? -2π * Ti[i] * dni_dpsi(ψ) / (chrg_ion * chi1 * ni[i]) : 0.0 @@ -400,7 +447,7 @@ function load_kinetic_profiles(kinetic_file::AbstractString; nue = zeros(Float64, nkin + 1) zeff = zeros(Float64, nkin + 1) - for i in 1:(nkin + 1) + for i in 1:(nkin+1) n_i = ni[i] n_e = ne[i] T_i = Ti[i] @@ -409,14 +456,15 @@ function load_kinetic_profiles(kinetic_file::AbstractString; z = n_e > 0 ? zimp - (n_i / n_e) * zi * (zimp - zi) : Float64(zimp) zpitch = 1.0 + (1.0 + mimp) / (2.0 * mimp) * zimp * (z - 1.0) / (zimp - z) - # Coulomb logarithm (NRL formulary form), natural log; n_e in 1e20 m^-3, T_e in keV. - ll = 17.3 - 0.5 * log(n_e / 1.0e20) + 1.5 * log(T_e / 1.602e-16) + ll = _coulomb_log(n_e, T_e) loglam[i] = ll + # Test-particle z² pitch-angle (Krook deflection) scaling — Logan & Park 2013 Eq. 6; + # generalizes the implicit zi=1 form; identical to the single-species limit of resolve_ntv_species. nui[i] = T_i > 0 ? - (zpitch / 3.5e17) * n_i * ll / (sqrt(1.0 * mi) * (T_i / 1.602e-16)^1.5) : 0.0 + (zpitch / _NU_PREFAC) * zi^2 * n_i * ll / (sqrt(1.0 * mi) * (T_i / _KEV_J)^1.5) : 0.0 nue[i] = T_e > 0 ? - (zpitch / 3.5e17) * n_e * ll / (sqrt(me / mp) * (T_e / 1.602e-16)^1.5) : 0.0 + (zpitch / _NU_PREFAC) * n_e * ll / (sqrt(_ME / _MP) * (T_e / _KEV_J)^1.5) : 0.0 zeff[i] = z end diff --git a/src/KineticForces/CalculatedKineticMatrices.jl b/src/KineticForces/CalculatedKineticMatrices.jl index 76b66008..71f58f0b 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,9 +57,9 @@ function compute_calculated_kinetic_matrices( ffs_intr, metric, ffit; - kf_ctrl::KineticForcesControl = KineticForcesControl(), + kf_ctrl::KineticForcesControl=KineticForcesControl(), kinetic_profiles::Equilibrium.KineticProfileSplines, - species::Union{Nothing,AbstractVector} = nothing, + species::Union{Nothing,AbstractVector}=nothing ) xs = metric.xs mpsi = length(xs) @@ -96,9 +99,9 @@ function compute_calculated_kinetic_matrices( # are read-only and safely shared through deepcopy semantics. 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] @@ -107,17 +110,17 @@ function compute_calculated_kinetic_matrices( # from kf_ctrl. The kinetic W/torque matrices are additive over species, so we # accumulate (+=) each species' block-diagonal contribution. splist = species === nothing ? - [(z=kf_ctrl.zi, m=kf_ctrl.mi, electron=kf_ctrl.electron, profiles=kinetic_profiles)] : species + [Equilibrium.ResolvedNTVSpecies(kf_ctrl.zi, kf_ctrl.mi, kf_ctrl.electron, "single", kinetic_profiles)] : species for sp in splist 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] for in_idx in 1:npert n = ffs_intr.nlow + in_idx - 1 fill!(full_w, 0) diff --git a/src/KineticForces/Compute.jl b/src/KineticForces/Compute.jl index 869f16a8..d2ce59ec 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,7 +97,7 @@ 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 # The outer ψ-integral (the QuadGK batch / ψ-node loop) stays serial: QuadGK's refine @@ -133,9 +135,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 @@ -205,7 +207,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 @@ -218,9 +220,12 @@ end Sum per-species [`KineticForcesState`](@ref) results into a single total (τ = Σ_s τ_s). Per method: the scalar `total_torque`/`total_energy` are summed exactly; the dT/dψ profile is summed -by evaluating each species' `torque_profile` interpolant on the sorted union of the species ψ -grids, and the cumulative T(ψ) is re-integrated (trapezoid) from that summed profile. All species -share the same ψ-integration range (`ctrl.psilims`), so the grids differ only in adaptive nodes. +by linearly interpolating each species' own `(psi_grid, dtdpsi)` arrays onto the sorted union of +the species ψ grids (zero outside a species' range), and the cumulative T(ψ) is re-integrated +(trapezoid) from that summed profile. All species share the same ψ-integration range +(`ctrl.psilims`), so the grids differ only in adaptive nodes. The combined `MethodResult` carries +only the summed scalars and profile; per-species diagnostics (`torque_profile`, `records`, +`panel_psis`, `resonance_psis`) are not aggregated and are left at their defaults. """ function combine_species_states(states::AbstractVector{KineticForcesState}) combined = KineticForcesState() @@ -236,7 +241,9 @@ function combine_species_states(states::AbstractVector{KineticForcesState}) dtdpsi = zeros(ComplexF64, length(grid)) for r in results length(r.psi_grid) >= 2 || continue - o = sortperm(r.psi_grid); xs = r.psi_grid[o]; ys = r.dtdpsi[o] + o = sortperm(r.psi_grid) + xs = r.psi_grid[o] + ys = r.dtdpsi[o] for (k, ψ) in enumerate(grid) (ψ < xs[1] || ψ > xs[end]) && continue j = searchsortedlast(xs, ψ) @@ -274,15 +281,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 @@ -361,7 +369,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 7b990762..9dcc6a58 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`. """ @@ -61,7 +62,7 @@ masses). The density is given by **exactly one** of: for measured profiles with distinct shapes. Both resolve to a per-species density profile `n_s(ψ)` at load time (one downstream path). -Constructed from a TOML `[[KineticForces.ion_species]]` table via [`convert`]. +Constructed from a TOML `[[KineticForces.ion_species]]` table via `convert`. """ struct IonSpecies z::Int @@ -83,7 +84,14 @@ Base.convert(::Type{Vector{IonSpecies}}, v::AbstractVector) = User-facing control parameters from the TOML `[KineticForces]` section. Configures which NTV methods to run, species parameters, tolerances, and output options. +`ion_species` selects the plasma model: empty (default) ⇒ a single main ion from `zi`/`mi` +(unchanged behaviour); non-empty ⇒ a multi-main-ion run where the NTV is computed per +`IonSpecies` under one shared full-composition Zeff and summed. For a multi-ion run the kinetic +file's `n_i` column is the **total** main-ion density, split across species by their `fraction` +or explicit `density` profile. + Constructed via keyword arguments or from a TOML dict: + ```julia ctrl = KineticForcesControl(; (Symbol(k) => v for (k, v) in inputs["KineticForces"])...) ``` @@ -118,12 +126,7 @@ ctrl = KineticForcesControl(; (Symbol(k) => v for (k, v) in inputs["KineticForce zimp::Int = 6 # Impurity charge mimp::Int = 12 # Impurity mass electron::Bool = false # Include electron contribution - # Multi-main-ion plasma: list of IonSpecies, each with its own z, m, and density - # (fraction of the total main-ion n_i, or an explicit per-species profile). When - # non-empty the NTV is computed per species with one shared full-composition Zeff and - # summed. Empty (default) ⇒ single main ion from zi/mi (unchanged behaviour). The - # kinetic file's n_i column is the TOTAL main-ion density for a multi-ion run. - ion_species::Vector{IonSpecies} = IonSpecies[] + ion_species::Vector{IonSpecies} = IonSpecies[] # multi-main-ion set (see docstring); empty ⇒ single ion # Mode numbers nn::Int = 1 # Toroidal mode number @@ -210,12 +213,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 @@ -305,17 +309,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 @@ -325,19 +329,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 @@ -406,9 +412,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) @@ -428,8 +434,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 @@ -456,9 +462,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) @@ -528,8 +534,8 @@ Accumulated results from all KineticForces computations. Written to gpec.h5 under the "kinetic_forces" 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/test/runtests_multiion.jl b/test/runtests_multiion.jl index 5adb97f1..d02bb7bf 100644 --- a/test/runtests_multiion.jl +++ b/test/runtests_multiion.jl @@ -1,9 +1,14 @@ +using HDF5 + @testset "Multi-ion NTV" begin EQ = GeneralizedPerturbedEquilibrium.Equilibrium KF = GeneralizedPerturbedEquilibrium.KineticForces @testset "multi_ion_composition" begin - ne = 1.24e20; ni = 1.13e20; zimp = 6; mimp = 12 + ne = 1.24e20 + ni = 1.13e20 + zimp = 6 + mimp = 12 # Reduces exactly to the single-ion formula Zeff = zimp - (ni/ne)*zi*(zimp-zi). for zi in (1, 2) zeff_old = zimp - (ni / ne) * zi * (zimp - zi) @@ -31,8 +36,11 @@ @testset "resolve_ntv_species" begin # Synthetic ASCII kinetic profile: psi, n_i(total), n_e, T_i, T_e, omega_E. psi = collect(0.0:0.05:1.0) - Ni = @. 1.0e20 * (1 - 0.5psi); ne = @. 1.1e20 * (1 - 0.5psi) - Ti = @. 2000.0 * (1 - 0.5psi); Te = @. 2500.0 * (1 - 0.5psi); wE = @. 1.0e4 * (1 - psi) + Ni = @. 1.0e20 * (1 - 0.5psi) + ne = @. 1.1e20 * (1 - 0.5psi) + Ti = @. 2000.0 * (1 - 0.5psi) + Te = @. 2500.0 * (1 - 0.5psi) + wE = @. 1.0e4 * (1 - psi) f = tempname() * ".txt" open(f, "w") do io println(io, "# psi ni ne Ti Te wexb") @@ -43,7 +51,7 @@ # Single species [z=1,m=2,fraction=1] reproduces load_kinetic_profiles. single = EQ.load_kinetic_profiles(f; zi=1, zimp=6, mi=2, mimp=12) - sp1 = EQ.resolve_ntv_species(f, [KF.IonSpecies(z=1, m=2, fraction=1.0)]; electron=false, zimp=6, mimp=12) + sp1 = EQ.resolve_ntv_species(f, [KF.IonSpecies(; z=1, m=2, fraction=1.0)]; electron=false, zimp=6, mimp=12) main = sp1[1] for ψ in (0.3, 0.6) @test main.profiles.ni_spline(ψ) ≈ single.ni_spline(ψ) rtol = 1e-8 @@ -52,46 +60,59 @@ end # D-T-e set = {D, T, impurity, electron}; ν scales as 1/√m at fixed z, T. - sp = EQ.resolve_ntv_species(f, [KF.IonSpecies(z=1, m=2, fraction=0.5), KF.IonSpecies(z=1, m=3, fraction=0.5)]; + sp = EQ.resolve_ntv_species(f, [KF.IonSpecies(; z=1, m=2, fraction=0.5), KF.IonSpecies(; z=1, m=3, fraction=0.5)]; electron=true, zimp=6, mimp=12) @test length(sp) == 4 @test count(s -> !s.electron && s.z == 1, sp) == 2 # D, T @test any(s -> s.z == 6 && !s.electron, sp) # impurity @test any(s -> s.electron, sp) # electron - iD = findfirst(s -> s.m == 2, sp); iT = findfirst(s -> s.m == 3, sp) + iD = findfirst(s -> s.m == 2, sp) + iT = findfirst(s -> s.m == 3, sp) @test sp[iD].profiles.nui_spline(0.5) / sp[iT].profiles.nui_spline(0.5) ≈ sqrt(3 / 2) rtol = 1e-6 # Exactly one of {fraction, density} required. - @test_throws ErrorException EQ.resolve_ntv_species(f, [KF.IonSpecies(z=1, m=2)]; zimp=6, mimp=12) + @test_throws ErrorException EQ.resolve_ntv_species(f, [KF.IonSpecies(; z=1, m=2)]; zimp=6, mimp=12) rm(f) end @testset "explicit per-species profiles (HDF5)" begin psi = collect(0.0:0.05:1.0) - Ni = @. 1.0e20 * (1 - 0.5psi); ne = @. 1.1e20 * (1 - 0.5psi) - Ti = @. 2000.0 * (1 - 0.5psi); Te = @. 2500.0 * (1 - 0.5psi); wE = @. 1.0e4 * (1 - psi) + Ni = @. 1.0e20 * (1 - 0.5psi) + ne = @. 1.1e20 * (1 - 0.5psi) + Ti = @. 2000.0 * (1 - 0.5psi) + Te = @. 2500.0 * (1 - 0.5psi) + wE = @. 1.0e4 * (1 - psi) h5 = tempname() * ".h5" HDF5.h5open(h5, "w") do file - file["psi"] = psi; file["n_e"] = ne; file["T_i"] = Ti; file["T_e"] = Te; file["omega_E"] = wE - file["n_i"] = Ni; file["n_D"] = 0.6 .* Ni; file["n_T"] = 0.4 .* Ni # unequal explicit shapes + file["psi"] = psi + file["n_e"] = ne + file["T_i"] = Ti + file["T_e"] = Te + file["omega_E"] = wE + file["n_i"] = Ni + file["n_D"] = 0.6 .* Ni + file["n_T"] = 0.4 .* Ni # unequal explicit shapes end - sp = EQ.resolve_ntv_species(h5, [KF.IonSpecies(z=1, m=2, density="n_D"), KF.IonSpecies(z=1, m=3, density="n_T")]; + sp = EQ.resolve_ntv_species(h5, [KF.IonSpecies(; z=1, m=2, density="n_D"), KF.IonSpecies(; z=1, m=3, density="n_T")]; electron=false, zimp=6, mimp=12) - iD = findfirst(s -> s.m == 2, sp); iT = findfirst(s -> s.m == 3, sp) + iD = findfirst(s -> s.m == 2, sp) + iT = findfirst(s -> s.m == 3, sp) @test sp[iD].profiles.ni_spline(0.5) ≈ 0.6 * 1.0e20 * (1 - 0.25) rtol = 1e-6 @test sp[iT].profiles.ni_spline(0.5) ≈ 0.4 * 1.0e20 * (1 - 0.25) rtol = 1e-6 # A missing named profile errors clearly. - @test_throws ErrorException EQ.resolve_ntv_species(h5, [KF.IonSpecies(z=1, m=2, density="n_He")]; zimp=6, mimp=12) + @test_throws ErrorException EQ.resolve_ntv_species(h5, [KF.IonSpecies(; z=1, m=2, density="n_He")]; zimp=6, mimp=12) rm(h5) end @testset "combine_species_states" begin # Combined profile = Σ species dT/dψ interpolated onto the union grid; total = Σ total. - mk(psi, dt) = (st = KF.KineticForcesState(); + mk(psi, dt) = ( + st = KF.KineticForcesState(); st.method_results["fgar"] = KF.MethodResult(; method="fgar", nn=1, total_torque=ComplexF64(sum(dt)), psi_grid=collect(Float64, psi), dtdpsi=ComplexF64.(dt), t_cumulative=zeros(ComplexF64, length(dt))); - st) + st + ) s1 = mk(0.0:0.1:1.0, fill(1.0, 11)) s2 = mk(0.05:0.1:0.95, fill(2.0, 10)) # different (offset) grid c = KF.combine_species_states([s1, s2]) From b864d38d350f78864bab25541fcfb3755c4c654b Mon Sep 17 00:00:00 2001 From: logan-nc Date: Fri, 31 Jul 2026 14:59:52 -0400 Subject: [PATCH 17/17] AGENTS - MEMORY - Record physics-reviewer findings from multi-ion NTV work MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit - kinetic_ntv_map.md: multi-ion NTV audit — test-particle z² justification, field-density (zpitch·n_main vs n_e·Zeff) reconciliation, and the electron- driven self-consistent δW sign-flip analysis (physics, not a double-count). - reg_spot_regularization.md: verdict that the reg_spot field-reconstruction smoothing port is correct and distinct from the singfac_min ODE crossing gate. - MEMORY.md: index pointer to the new reg_spot audit. Co-Authored-By: Claude Opus 4.8 (1M context) Claude-Session: https://claude.ai/code/session_01SsLUydP2gJbE1tGoaHDiS1 --- .../fortran-physics-reviewer/MEMORY.md | 1 + .../kinetic_ntv_map.md | 56 +++++++++++++++++++ .../reg_spot_regularization.md | 24 ++++++++ 3 files changed, 81 insertions(+) create mode 100644 .claude/agent-memory/fortran-physics-reviewer/reg_spot_regularization.md diff --git a/.claude/agent-memory/fortran-physics-reviewer/MEMORY.md b/.claude/agent-memory/fortran-physics-reviewer/MEMORY.md index 1a3bc1fc..9d17a0bd 100644 --- a/.claude/agent-memory/fortran-physics-reviewer/MEMORY.md +++ b/.claude/agent-memory/fortran-physics-reviewer/MEMORY.md @@ -5,3 +5,4 @@ - [KineticForces (NTV) Audit Checklist](kinetic_ntv_map.md) — pentrc->KineticForces map, Logan 2015 matrices, what to verify - [InnerLayer (Resistive) Audit Checklist](resistive_layer_map.md) — rmatch->InnerLayer map, GGJ Wasow basis / Δ′, what to verify - [Galerkin Δ′ Assembly Map](galerkin_assembly_map.md) — gal.f<->GalerkinAssembly.jl; resonant sign chain verified; PASS +- [reg_spot vs singfac_min audit](reg_spot_regularization.md) — reg_spot field-recon smoothing CORRECT; distinct from singfac_min ODE gate diff --git a/.claude/agent-memory/fortran-physics-reviewer/kinetic_ntv_map.md b/.claude/agent-memory/fortran-physics-reviewer/kinetic_ntv_map.md index 414be7a9..a3b73369 100644 --- a/.claude/agent-memory/fortran-physics-reviewer/kinetic_ntv_map.md +++ b/.claude/agent-memory/fortran-physics-reviewer/kinetic_ntv_map.md @@ -31,3 +31,59 @@ kinetic energy/matrices from trapped-particle nonambipolar transport. Fortran co - **Normalization**: torque Im(T)=2n·δW_k (Eq. 19) and diamagnetic-frequency sign/factor conventions match Fortran `torque.F90`. - **Mode indexing**: m, n ranges and block packing over n stay consistent with ForceFreeStates. - **Method variants**: FGAR/TGAR/PGAR/RLAR/CLAR/FCGL/TMM/WMM each present, not stubbed to a single fallback. + +## Multi-ion (D-T) NTV — composition/collisionality (2026-07) +Both Fortran PENTRC (`inputs.f90:236-243`, `read_kin`) and Julia (`KineticProfiles.jl:261-278`) +support ONE main ion (zi,mi) + ONE impurity (zimp,mimp) per run. Correct multi-main-ion NTV is +an EXTENSION beyond Fortran, but consistent with Logan-Park 2013's pitch-angle (Lorentz) model. +- Zeff = Σ_s Z_s² n_s / n_e; quasineutrality n_e = Σ_s Z_s n_s (all ions incl. impurity). +- Bug in current D-T split (run D, then T, each with ni=Ni/2): the `ni` column sets BOTH species + density AND Zeff via z = zimp-(n_i/n_e)zi(zimp-zi). With ni=Ni/2 it treats the missing half as + high-Z impurity → Zeff≈3.7 instead of true ≈1–1.4. Corrupts zpitch → nue,nui. +- Correct collisionality: ν_a ∝ Z_a² lnΛ · Σ_b n_b Z_b² /(√m_a T_a^{3/2}) = ∝ Z_a² lnΛ·(n_e·Zeff_true). + The zpitch(Zeff) polynomial is a main-ion+impurity closure (momentum-restoring correction); the + MINIMAL fix is to feed the TRUE (full-composition) Zeff into the existing zpitch/ν formulas. +- Additivity: τ = Σ_s τ_s. Lorentz operator is additive over field species; species couple only + through shared δB, shared ω_E, shared Zeff/ν. Additive at this theory's order. +- What changes single→multi: Zeff, zpitch, nue, nui (→nueff). UNCHANGED for equal-shape D-T: + wdian/wdiat (log-derivative, density factor cancels), wtran/wbhat/wdhat/wgyro (already per-species), + and the resonant-density prefactor (Ni/2 per species is correct). +- ASIDE (separate fidelity bug, same code block): Fortran `inputs.f90:238` uses natural log for lnΛ; + Julia `KineticProfiles.jl:270` uses log10 — diverges away from the n=1e20,T=1keV reference point. + +## Multi-ion NTV — full-composition species set (2026-07 audit, PASS-with-caveats) +Reviewed `resolve_ntv_species` (KineticProfiles.jl ~214-266) + `compute_calculated_kinetic_matrices` +(CalculatedKineticMatrices.jl ~96-151). Verdict: physics is sound. +- ν_s field-density RECONCILIATION (supersedes the earlier "should be n_e·Zeff" note above): + code uses `ν_s = (zpitch/3.5e17)·z_s²·n_main·lnΛ/(√m_s·T_i^1.5)`, i.e. field density = zpitch·n_main + (n_main = Σ MAIN-ion densities, no impurity, unweighted). This is CORRECT and MORE faithful to + single-ion PENTRC than n_e·Zeff: PENTRC's design is zpitch·n_i, NOT Σ_b n_b Z_b². Numerically + zpitch·n_main ≈ n_e·Zeff (e.g. Zeff=1.5,C6 → 1.55 vs 1.5); the gap IS the intended momentum- + restoring correction that zpitch(Zeff) carries. So n_main vs n_e·Zeff is immaterial at D-T (z=1). +- z_s² test-particle factor: CORRECT and correctly placed (deflection freq ∝ test charge²). Single-ion + had no z² only because zi=1. Reduces EXACTLY to single-ion nui for one z=1,fraction=1 ion (verified). +- Impurity as its own test species: field density zpitch·n_main is NOT undercounting — the impurity's + z_imp² contribution is already folded into Zeff inside zpitch. CAVEAT: zpitch is strictly a main-ion + momentum-restoring closure; reusing it for the impurity/electron ν is an approximation beyond the + single-ion theory. Acceptable (impurity δW ∝ n_imp is small); worth a one-line annotation. +- Electron descriptor passes `ns[1]` (first ion's density) as its `ni_spline` — HARMLESS: the kernel + `_setup_surface_state` (Torque.jl:659-664) reads `ne_spline` (full shared n_e) when electron=true, + never ni_spline. Electron n_s = full n_e (correct). Cosmetic smell only. +- Self-consistent δW summation (CalculatedKineticMatrices.jl:112-148): kw_flat/kt_flat are PURELY the + kinetic matrices (Logan 7.30-7.35), every term ∝ species phase-space density n_s·f0_s — NO species- + independent baseline. Fluid F,K,G added ONCE downstream in _compute_fkg_matrices!, outside the + species loop. So `+=` over species is clean additivity; NOTHING in kw is wrongly ×species-count. + → TC-24 n=3 δW +0.066(D) → −0.10(D+T+C+e) sign flip is PLAUSIBLE physics near marginal stability, + NOT a double-count. Dominant driver: the newly-ON electron channel (full n_e, opposite precession). + DECISIVE cheap diagnostic to confirm: run D(½)+T(½) with electron=OFF, impurity absent — should + return ≈ +0.066 (single-ion D). If D+T alone ≈ +0.066, the whole shift is electron+impurity = physical. + +## Single-ion nui vs multi-species z_s² (#339, 2026-07 decision) +Fortran PENTRC `inputs.f90:240-241` single-ion nui = (zpitch/3.5e17)·n_i·lnΛ/(√mi·T_i^1.5) has +NO explicit zi² (implicitly assumes zi=1, main ion hydrogenic). Julia `load_kinetic_profiles` +(~L416) is a faithful exact port. Multi-species `_nu` (~L253) adds explicit test-particle z_s² +(the physically correct pitch-angle/Lorentz form; zpitch is the field-side momentum-restoring +factor of Zeff, independent of test charge — no double count). RECOMMENDATION: option (a) ADD +zi² to single-ion nui. It is numerically identical for the default zi=1 (regression byte-identical), +makes single-ion the true 1-species special case of `_nu`, and only "deviates" from Fortran in the +exotic zi≠1 case where Fortran is physically wrong anyway. Annotate as a documented improvement. diff --git a/.claude/agent-memory/fortran-physics-reviewer/reg_spot_regularization.md b/.claude/agent-memory/fortran-physics-reviewer/reg_spot_regularization.md new file mode 100644 index 00000000..137301e2 --- /dev/null +++ b/.claude/agent-memory/fortran-physics-reviewer/reg_spot_regularization.md @@ -0,0 +1,24 @@ +--- +name: reg_spot vs singfac_min regularization audit +description: Verdict on GPEC-Julia reg_spot (field-reconstruction smoothing) port and its distinction from singfac_min (ODE crossing gate) +metadata: + type: project +--- + +## Verdict (audited FieldReconstruction.jl, no local Fortran repo available) +`reg_spot` port judged CORRECT-WITH-CAVEATS (caveat = Fortran verified from GPEC `xm*` convention + physics + known default, not from on-disk source). + +## reg_spot — field-reconstruction smoothing (PerturbedEquilibrium) +- Factor form: `reg_factor = singfac²/(singfac²+reg_spot²)`, singfac = m - n·q. Correct GPEC form. → 0 at rational surface, → 1 away. Applied in singfac-space, NOT ψ-space. +- Default 5e-2 = 0.05 (PerturbedEquilibriumStructs.jl:50). Matches GPEC gpec_input namelist default. +- Quantities regularized (matches GPEC `xm*` modified-quantity convention): + - xmp1 = ξ^ψ' (FieldReconstruction.jl:393-395) + - xms (clebsch_alpha) computed FROM regularized xmp1 via xms=-A⁻¹(B·xmp1_reg+C·xsp) (line 398-415, gpeq_sol) + - xmt/xmz = regularized ξ^θ,ξ^ζ (line 584-594, gpeq_contra); xwt/xwz kept unregularized + - b^θ_reg,b^ζ_reg inherit reg via xmp1/xms (compute_modified_field_modes) +- LEFT UNREGULARIZED: ξ^ψ primitive (xsp / clebsch_psi = copy, line 356; xwp uses raw xsp). CORRECT — ξ^ψ (normal disp) is finite at rational surface; only 1/singfac-divergent tangential/derivative quantities are softened. + +## singfac_min — DISTINCT mechanism (ForceFreeStates ODE crossing gate) +- Default 1e-4 (ForceFreeStatesStructs.jl:267), "Matches Fortran STRIDE". Gates singular-surface crossings during EL ODE integration (EulerLagrange.jl:466-519, GeneralizedPerturbedEquilibrium.jl:426). +- SingularCoupling.jl:204 uses `spot = 5e-4` as an offset (5e-4/|n·q1|) to evaluate bwp1 at lpsi/rpsi — labeled "matches Fortran default singfac_min". This is the ψ-offset for one-sided derivative eval, a THIRD distinct use, not reg_spot. +- reg_spot (0.05, singfac-space smoothing of reconstructed fields) and singfac_min (1e-4, ODE crossing gate) are separate; both ported.