KineticForces - NEW FEATURE - Multi-main-ion NTV (summed ions + impurity + electrons, both paths)

[logan-nc]

Summary

Adds correct multi-main-ion NTV to KineticForces. A plasma can now declare an arbitrary list of main-ion species (each with its own z, m, and density), and the NTV is computed per species on one shared full-composition Zeff and summed — τ = Σ_s τ_s — over the main ions, the neutrality-closing impurity, and the electrons. This applies identically to both NTV paths:

• the KineticForces ψ-quadrature (post-PE diagnostic), and
• the self-consistent DCON kinetic matrices (Mode A) — the kinetic W/torque matrices are additive over species, accumulated in a species loop.

Previously only a single main ion was supported (D-T had to be run per-species externally and summed, which double-used the n_i column as both the species density and the Zeff source — giving each ion run the wrong Zeff ≈ 3.7 instead of the composition Zeff ≈ 1.4). This is an extension beyond Fortran PENTRC (also single main ion), consistent with Logan & Park PoP 2013.

Design (physics-reviewed)

• One canonical path: each species resolves to a per-species density profile n_s(ψ); fraction (share of the total n_i) and explicit HDF5 density arrays both feed it (normalize-at-the-boundary).
• Two distinct densities kept separate: resonant n_s (per-species torque prefactor) vs the collision field density (total n_main + shared composition Zeff + species z_s², m_s).
• Composition multi_ion_composition: Zeff = (Σ z_s² n_s + z_imp² n_imp)/n_e, impurity closes quasineutrality; reduces exactly to the single-ion formula for one z=1 ion.

API

[KineticForces]
kinetic_file = "kinetic.dat"   # n_i column = TOTAL main-ion density
electron = true
zimp = 6; mimp = 12            # impurity (closes quasineutrality; summed as its own species)

[[KineticForces.ion_species]]
z = 1; m = 2; fraction = 0.5   # Deuterium
[[KineticForces.ion_species]]
z = 1; m = 3; fraction = 0.5   # Tritium

Empty ion_species ⇒ single main ion from zi/mi — bit-for-bit unchanged. Explicit measured profiles via named HDF5 datasets (density = "n_D").

Validation

• 25 unit tests (test/runtests_multiion.jl): composition (exact single-ion reduction; 50/50 D-T Zeff; quasineutrality), TOML Vector{Any}→IonSpecies conversion, resolver (single-species reduction, D-T-e set, ν∝1/√m, one-of-{fraction,density}), explicit HDF5 profiles.
• Regression case solovev_kinetic_multiion: tracks the summed total + each per-species NTV end-to-end.
• End-to-end (ITER n=3 TC-24): KF path summed 4 species (D+T+C+e), total = exact per-species sum; Mode A self-consistent kinetic-matrix path completed with the summed matrices.
• fortran-physics-reviewer sign-off: collision-frequency z_s²/n_main form CORRECT (n_main is the faithful Fortran-matching choice, not n_e·Zeff); the self-consistent δW shift is additive-channel physics (ideal F/K/G added once downstream — no per-species double-count), not a summation bug.

Stacking note

This branch includes two bugfix commits that are their own PRs and should merge first; they deduplicate on rebase:

• #337 — Coulomb log natural-log fix (prerequisite for correct collisionality).
• #336 — IntegrationChunk arity fix (prerequisite for any Mode A / self-consistent kinetic run).

Not included (follow-ups)

• Scaling factors (density_factor, …) with a multi-ion list error out (documented).
• zpitch is a main-ion closure applied approximately to impurity/electron test species (annotated).

🤖 Generated with Claude Code

https://claude.ai/code/session_01SsLUydP2gJbE1tGoaHDiS1

[logan-nc]

@claude review this pr

[claude]

Review: Multi-main-ion NTV

Reviewed the diff (18 commits: KineticProfiles.jl, Compute.jl, CalculatedKineticMatrices.jl, KineticForcesStructs.jl, the new example/regression case, and runtests_multiion.jl).

Potential physics inconsistency

resolve_ntv_species's per-species collision frequency includes an explicit z_s² factor:

_nu(zsp, msp) = [... (zpitch[i]/3.5e17) * zsp^2 * n_main[i] * loglam[i] / (sqrt(Float64(msp)) * ...) ...]

but the pre-existing single-ion path in load_kinetic_profiles does not carry a zi² term:

nui[i] = T_i > 0 ? (zpitch/3.5e17) * n_i * ll / (sqrt(1.0*mi) * (T_i/1.602e-16)^1.5) : 0.0

For the flagship D/T case (z=1 for both species) this cancels (1²=1), so the headline result is unaffected, and the unit tests only exercise z=1 species. But for any non-hydrogenic main ion (e.g. He, z=2) the two paths would diverge by a factor of z². The resolve_ntv_species docstring even hedges this ("Reduces to load_kinetic_profiles for one z=1 main ion with fraction=1") rather than claiming an unconditional reduction, which suggests the discrepancy is known but not explained. Worth either (a) confirming with the physics reviewer which form is correct and fixing the other, or (b) adding an explicit docstring/comment note that the two formulas intentionally disagree for z≠1 and why.

Test coverage gap

The new regression case (solovev_kinetic_multiion) and its gpec.toml set kinetic_source = "fixed" / kinetic_factor = 0.0, so they only exercise the KineticForces ψ-quadrature species-summing path (resolve_ntv_species + combine_species_states). The other half of this PR — the self-consistent DCON kinetic-matrix accumulation in compute_calculated_kinetic_matrices (species loop wrapping the threaded ψ-loop, += accumulation into kw_flat/kt_flat) — is only verified per the PR description via a manual, non-committed ITER TC-24 run. Given this is a non-trivial change (nested Threads.@threads per species, block-diagonal accumulation), it would be worth adding a small regression/unit case with kinetic_source = "calculated" and a non-empty ion_species list so this path doesn't silently regress later.

Minor style

• The new ion_species field comment in KineticForcesControl is a 5-line block that largely duplicates the IonSpecies docstring. Consider trimming to 1-2 lines and pointing at IonSpecies, per the project's "concise comments" convention.
• combine_species_states's combined MethodResult drops torque_profile, records, panel_psis, and resonance_psis (left at their zero-value/nothing defaults). That's a reasonable choice, but it isn't called out in the docstring — worth a one-line note so a future reader of e.g. kinetic_forces/<method>/panel_psis in the summed group doesn't wonder why it's empty.

What looks solid

• multi_ion_composition reduces exactly to the existing single-ion Zeff formula (verified algebraically in the diff and by the new unit test) — good regression discipline for a formula change.
• The species === nothing fallback in compute_calculated_kinetic_matrices and the isempty(kf_ctrl.ion_species) fallback in main_from_inputs keep the single-ion path bit-for-bit unchanged, which is the right way to land a feature like this without disturbing existing users.
• The IonSpecies / TOML Vector{Any} → Vector{IonSpecies} convert overload, and the one-of-{fraction,density} validation, are handled correctly and are unit-tested.
• The separate IntegrationChunk arity fix addresses a real latent MethodError on a rarely-exercised serial+kinetic-singular-surface path — useful fix, currently inert since the chunk is discarded there.

Process note

Per CLAUDE.md, the regression harness should be run and its report included before merging any substantive change. The PR description shows unit test results and a manual ITER end-to-end check, but no regress --cases diiid_n1,solovev_n1,solovev_multi_n,solovev_kinetic_multiion --refs develop,local table. Given this PR also carries the Coulomb-log natural-log fix, a combined run would confirm nothing besides the already-noted diiid_n1 shift moved unexpectedly.

[logan-nc]

Review changes applied (5f9f40ee) + reviewer memory (b864d38d)

Physics

• ν_s now carries the test-particle z² (Krook deflection, Logan & Park 2013 Eq. 6). Confirmed as the correct form (deflection-off-background, not self-collision Z⁴) in a literature review. Arithmetically identical for z=1, so it only generalizes impurity/multi-species collisionality and makes the single-ion path the exact 1-species limit of the multi-ion path.
• Hard error on any negative density (bad input or cubic-resample overshoot) on nₑ, total nᵢ, each ion species, and the neutrality-closing impurity — replaces silent clamping.

Maintainability

• resolve_ntv_species returns concrete-eltype ResolvedNTVSpecies structs → type-stable species loop in the self-consistent kinetic-matrix build.
• Fixed combine_species_states docstring (linear-interp of dtdpsi; documents the dropped per-species fields).
• Shared physical constants + Coulomb-log helper across load_kinetic_profiles/resolve_ntv_species; multi-ion contract moved into struct docstrings; doc links normalized.
• runtests_multiion.jl made self-contained (using HDF5).

Verification

• Multi-ion suite: 28/28 passing. Formatted with pinned JuliaFormatter v1.0.62.
• Regression: numerically inert. diiid_n1 NTV torque FGAR + dW FGAR both OK (z=1 ⇒ z²=1). The multi-ion collision/composition formulas (_nu, n_imp, zeff, zpitch, nue) are byte-identical before/after this commit, so solovev_kinetic_multiion is unchanged. The ideal-MHD deltas seen in a whole-branch-vs-develop diff come from the stacked 66cc9dee ising fix (ForceFreeStates - BUG FIX - Pass direction arg to IntegrationChunk in kinetic EL crossing #336), not from these commits.