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

[logan-nc]

Release note

• Audience: users
• Numerical impact: solovev_kinetic_calculated and its collisionless nuzero variant are re-baselined because their shared deck now declares a 50/50 D-T mix (Re(et[1]) +1.51%, Im(et[1]) +0.93%); diiid_n1, solovev_n1 and solovev_kinetic_ntv are unchanged (harness @ 0e21828)
• Migration: electron = true now adds electrons in addition to the ion species rather than computing them instead of the ions, and per-species NTV output moved from flat top-level kinetic_forces_ groups into the KineticForces/PerSpecies/ subtree.

A plasma can now declare an arbitrary list of main-ion species, each with its own charge, mass and density, and the NTV is computed per species on one shared full-composition Zeff and summed over the main ions, the quasineutrality-closing impurity and the electrons. This applies to both the post-PE ψ-quadrature diagnostic and the self-consistent kinetic matrices, and replaces the previous workaround of running each ion separately, which gave every run the wrong Zeff.

Regression report

regress --cases diiid_n1,solovev_n1,solovev_kinetic_ntv,solovev_kinetic_nuzero,solovev_kinetic_calculated,solovev_kinetic_multiion --refs develop,local (develop @ 14a3c9d6)

diiid_n1                     47 unchanged
solovev_n1                   21 unchanged
solovev_kinetic_ntv           6 unchanged
solovev_kinetic_nuzero        5 changed, 9 unchanged
solovev_kinetic_calculated    6 changed, 8 unchanged
solovev_kinetic_multiion      new case; develop ref fails as expected
                              ("Example directory not found"), local values recorded

solovev_kinetic_calculated      develop         local           Diff
total energy Re(et[1])          1.845905e+00    1.873742e+00    2.784e-02 (1.51%)
total energy Im(et[1])         -1.421059e+00   -1.434284e+00    1.322e-02 (0.93%)
plasma energy Re(ep[1])        -8.593965e+00   -8.563835e+00    3.013e-02 (0.35%)
vacuum energy Re(ev[1])         1.043987e+01    1.043758e+01    2.293e-03 (0.02%)

solovev_kinetic_nuzero (collisionless override of the same deck)
total energy Re(et[1])          2.245440e+00    2.243336e+00    2.103e-03 (0.09%)
total energy Im(et[1])         -2.034869e+00   -2.026028e+00    8.841e-03 (0.43%)

Every changed row is confined to the two cases sharing the calculated deck, which this PR deliberately switched to a 50/50 D-T mix so the multi-ion kinetic-matrix species loop is pinned end-to-end. diiid_n1's NTV torque and dW rows are exactly 0.0e+00, so the single-ion paths are untouched and the movement is the fixture change alone.

Notes for reviewers

Since the last review round this branch also merged develop and restored develop's #367 immutability freeze on KineticForcesControl, which an earlier develop merge had silently reverted here because the species loop mutated a deepcopy of the run control. Each species' control is now built with the @kwdef keyword constructor — the idiom that freeze commit documents. The multi-ion example produces bit-identical per-species torques before and after that rewrite.

---

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.
• Unified electron flag: electron = true always adds the electron NTV in addition to the ion species, on both the single-ion and multi-ion paths (a single-ion + electron run resolves as {main ion, impurity, electron} summed). The Fortran PENTRC electron-instead-of-ions behavior was an IO limitation and is deliberately not reproduced. electron = false with an empty ion_species list is bit-for-bit unchanged.
• Input validation: fraction lists can never sum above 1 (shares of n_i); an all-fraction list must sum to 1 (in mixed fraction/density lists the impurity content is set by the file's n_i/n_e deficit); every main-ion z < zimp (zpitch-closure pole at Zeff = zimp, with a near-pole @warn); hard error on any negative density.

API

[KineticForces]
kinetic_file = "kinetic.dat"   # n_i column = TOTAL main-ion density
electron = true                # adds electrons to the species sum
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 + electron = false ⇒ single main ion from zi/mi — bit-for-bit unchanged. Explicit measured profiles via named HDF5 datasets (density = "n_D"; datasets must be named n_* and round-trip through write_kinetic_h5). Docs: new "Multi-ion runs" section in docs/src/kinetic_forces.md.

Validation

• Unit tests (test/runtests_multiion.jl, 46): composition (exact single-ion reduction; 50/50 D-T Zeff; quasineutrality), TOML conversion, resolver (single-species reduction, D-T-e set, ν∝1/√m, one-of-{fraction,density}, unified electron semantics), explicit HDF5 profiles + write/read roundtrip, all validation guards, union-grid summing incl. duplicate-node guard.
• Regression cases: solovev_kinetic_multiion (ψ-quadrature path: summed total + per-species NTV) and solovev_kinetic_calculated now runs 50/50 D-T (self-consistent kinetic-matrix species loop; deliberate baseline change et[1]: 1.846−1.421i → 1.874−1.434i, with the collisionless nuzero override variant shifting accordingly).
• 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.
• Two independent review rounds applied (claude-bot + krystophny + a 3-agent physics/clean-code/performance audit): test-particle z² Krook form confirmed (deflection-off-background; exact z=1 reduction), resonant-n_s vs collision-n_main separation confirmed, additive matrix accumulation confirmed (no ideal F/K/G double-count), :static threading on the kinetic-matrix ψ-loop.
• Regression report vs develop: see PR comment.

Behavior changes (deliberate)

1. electron = true with empty ion_species now computes ion + impurity + electron NTV summed (was: electron only). No shipped example or regression case used the old mode.
2. solovev_kinetic_calculated baseline re-blessed for the D-T fixture change.

🤖 Generated with Claude Code

[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.

[logan-nc]

@krystophny could you review this (for TC24)?

[krystophny]

@logan-nc Three things before this goes in.

The fraction path in resolve_ntv_species never checks that the fractions sum to one. If they don't, the shortfall ends up in the quasineutrality impurity, which pushes zeff toward zimp, and zpitch has a pole there. A typo like 0.5/0.4 gives no error, just a collisionality that runs away. Either require sum == 1 within a tolerance, or say what a partial sum means.

write_kinetic_h5 doesn't write species_densities, but read_kinetic_h5 reads it, and docs/src/kinetic_forces.md tells people to author files with write_kinetic_h5. So the documented way of making a kinetic file can't produce one the new explicit-profile input accepts. Needs the put() loop.

The resolve_ntv_species docstring still says explicit per-species profiles are "not yet wired". They are, a few lines below.

The hard error on the scaling factors is fine.

[matt-pharr]

@logan-nc I am going to refrain from reviewing this until we have our group sit-downs next week. Unfortunately I have a lot of these PR review requests building up and cannot get them done fast enough with my thesis writing right now.

[logan-nc]

Review round 2 applied (six commits) + independent 3-agent audit

All three of @krystophny's items, the remaining claude-bot items, and the findings of a fresh independent review (physics / clean-code / performance agents over the full diff) are addressed. Point-by-point:

@krystophny's items

1. Fraction-sum validation — implemented as two hard errors (resolve_ntv_species):

• any list (mixed included) with Σ fraction > 1 + 1e-6 errors — fractions are shares of the total n_i, so over-allocation silently inflated every ν_s before;
• an all-fraction list additionally requires |Σ − 1| ≤ 1e-6.

A blanket Σ==1 check would have broken supported mixed fraction+density lists, where the legitimate Fortran semantics apply: the impurity content is set by the file's n_i/n_e deficit, not by a fraction shortfall (docstring now states this). Two new guards land in the same commit, both motivated by your zpitch-pole observation: every main-ion z < zimp is enforced (the closure has a pole at Zeff = zimp and goes negative beyond), and a @warn fires when max(Zeff) comes within 10% of zimp. multi_ion_composition also no longer parks the ne ≤ 0 fallback exactly on the pole. Unit-tested (0.9/1.2/mixed/z≥zimp/near-pole warning).

2. write_kinetic_h5 species densities — fixed with a roundtrip test (write → read → resolve_ntv_species consumes the file). Per-species datasets must be named n_*; the reader now correspondingly harvests only n_*-prefixed non-standard datasets, so future schema additions can't be silently swept in as densities.

3. Stale "not yet wired" docstring — fixed; the docstring now describes the implemented named-dataset lookup + resample behavior.

claude-bot items

• Mode A coverage: the existing solovev_kinetic_calculated fixture now declares a 50/50 D-T main-ion mix, pinning the self-consistent kinetic-matrix species accumulation end-to-end in the harness. Deliberate baseline change: et[1] moves 1.894−1.525i → 1.874−1.434i (T's higher mass reduces the kinetic damping contribution).
• ion_species comment duplication: already resolved at HEAD (single-line field comment + docstring); no change.

Independent audit findings applied

• Unified electron semantics (decision by @logan-nc): electron = true now always adds the electron NTV in addition to the ion species — a single-ion + electron run resolves through the species machinery as {main ion, quasineutrality impurity, electron}, summed. The old single-ion behavior (electron instead of ions) was an artifact of Fortran PENTRC's one-species-per-run IO, not physics. electron = false + empty ion_species is bit-for-bit unchanged, and no shipped example/regression case used the old electron-only mode.
• Threading: the kinetic-matrix ψ-loop indexed per-thread buffers by threadid() under default dynamic scheduling — a task-migration hazard multiplied by the new species loop. Now :static (matching the sibling quadrature loop), with per-species scalars hoisted so the closure captures concrete types.
• Docs: new "Multi-ion runs" section in docs/src/kinetic_forces.md (species table, fraction vs density, n_* schema, per-species output groups, electron flag); ν_s citation now points at the PENTRC numerical form rather than implying Logan & Park Eq. (6) carries the prefactor; the zpitch·n_main field term's z=1 validity domain is stated; t_cumulative documented as a diagnostic (endpoint ≠ exact summed Gauss-Kronrod total); example TOMLs fully annotated.
• Physics audit verdict on the rest: composition/ν forms, resonant-vs-collision density separation, and the additive kinetic-matrix accumulation (no F/K/G double-count) all confirmed sound; the develop merge conflict resolution verified clean.

Regression

Full run: regress --cases diiid_n1,solovev_n1,solovev_kinetic_ntv,solovev_kinetic_nuzero,solovev_kinetic_calculated,solovev_kinetic_multiion --refs develop,local (develop @ 1308b85e, single-thread, pinned manifest):

Case	Result
diiid_n1	48 unchanged (NTV torque FGAR + dW FGAR exactly 0.0 diff)
solovev_n1	21 unchanged
solovev_kinetic_ntv	6 unchanged (ψ-quadrature NTV bit-identical)
solovev_kinetic_calculated	6 changed — deliberate D-T fixture change (below)
solovev_kinetic_nuzero	7 changed — same fixture (reuses the calculated deck via [overrides])
solovev_kinetic_multiion	new case: no develop baseline (example dir doesn't exist on develop); local values recorded

The two changed cases are the one intended baseline change (the calculated deck now declares a 50/50 D-T main-ion mix to pin the Mode A species loop):

solovev_kinetic_calculated        develop         local           Diff
total energy Re(et[1])            1.845905e+00    1.873742e+00    2.784e-02 (1.51%)
total energy Im(et[1])           -1.421059e+00   -1.434284e+00    1.323e-02 (0.93%)
plasma energy Re(ep[1])          -8.593964e+00   -8.563835e+00    3.013e-02 (0.35%)
vacuum energy Re(ev[1])           1.043987e+01    1.043758e+01    2.293e-03 (0.02%)

solovev_kinetic_nuzero (collisionless override of the same deck)
total energy Re(et[1])            2.245439e+00    2.243336e+00    2.103e-03 (0.09%)
total energy Im(et[1])           -2.034869e+00   -2.026028e+00    8.841e-03 (0.43%)

Attribution was verified by isolation: re-running solovev_kinetic_nuzero against 61912a0b (this branch immediately before this review round) reproduces develop's values bit-for-bit — the branch's feature code is numerically inert on single-ion paths, and the only mover in this round is the fixture change itself.

[logan-nc]

@matt-pharr has been assigned to many reviews, so I'll hand this off to @ebursch since he is already developing things related to torques.

@ebursch this needs a human sanity check of the diff - sorry it's a bit lengthy but luckily the physics isn't hard. Ultimately, we are just trying to loop over multiple species in the NTV calculation.

[logan-nc]

Regression re-baselined against current develop (b3abe074)

The table I posted earlier is superseded: develop has since moved 49 commits (notably #364 refactor/hdf5-metadata), and this branch now merges it. Re-run with --refs develop,local, single thread, pinned manifest.

Case	Result
diiid_n1	48 unchanged
solovev_n1	21 unchanged
solovev_kinetic_ntv	6 unchanged (ψ-quadrature NTV bit-identical)
solovev_kinetic_calculated	6 changed — deliberate D-T fixture change
solovev_kinetic_nuzero	7 changed — same fixture (reuses the calculated deck via [overrides])
solovev_kinetic_multiion	new case; develop ref fails as expected (Example directory not found), local values recorded

Every changed row is confined to the two cases that share the calculated deck, which this PR deliberately switched to a 50/50 D-T main-ion mix so the multi-ion kinetic-matrix species loop is pinned end-to-end:

solovev_kinetic_calculated        develop         local           Diff
total energy Re(et[1])            1.845905e+00    1.873742e+00    2.784e-02 (1.51%)
total energy Im(et[1])           -1.421059e+00   -1.434284e+00    1.323e-02 (0.93%)
plasma energy Re(ep[1])          -8.593964e+00   -8.563835e+00    3.013e-02 (0.35%)
vacuum energy Re(ev[1])           1.043987e+01    1.043758e+01    2.293e-03 (0.02%)
ODE steps (total)                 1800            1801            1 (0.06%)

solovev_kinetic_nuzero (collisionless override of the same deck)
total energy Re(et[1])            2.245439e+00    2.243336e+00    2.103e-03 (0.09%)
total energy Im(et[1])           -2.034869e+00   -2.026028e+00    8.841e-03 (0.43%)
plasma energy Re(ep[1])          -8.485145e+00   -8.479030e+00    6.115e-03 (0.07%)

Attribution is pinned: diiid_n1's NTV torque FGAR and dW FGAR rows are exactly 0.0e+00, and solovev_kinetic_ntv is entirely unchanged, so the single-ion paths are untouched and the movement is the fixture change alone.

HDF5 schema conformance

Per-species NTV output now nests under KineticForces/PerSpecies/<label>/<method>/ per docs/development/hdf5-conventions.md (the flat kinetic_forces_<label> top-level groups violated the CamelCase group rule). Verified on a real 3-species run:

KineticForces/
  fgar/                    total_torque = 1.317868e-04
  PerSpecies/
    ion_z1_m2/fgar/        +7.794827e-05
    ion_z1_m3/fgar/        +8.206067e-05
    electron/fgar/         -2.822214e-05

Species sum to the total bitwise; every per-species dataset carries long_name + units; dTdpsi is native ComplexF64. runtests_h5_schema.jl 14/14, runtests_kinetic.jl 277/277, runtests_multiion.jl 52/52.

[krystophny]

Thanks, Nik. I rechecked the requested points and the multi-ion regression. We accidentally included an unnecessary CI hash refresh, reverted it, and everything is now fine to merge from my side.