Merge branch 'develop' into bugfix/kf-periodic-vpar-spline

Author: logan-nc

.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

.claude/agent-memory/fortran-physics-reviewer/kinetic_ntv_map.md

@@ -31,3 +31,60 @@ the Fortran `pentrc/` sources (file mapping in fortran_correspondence_map.md).
 - **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². The two are NOT
+  numerically equal (Zeff=1.5, C6, n_i/n_e=0.9 → zpitch·n_i = 1.225·n_e vs n_e·Zeff = 1.5·n_e, ~18%
+  apart); that gap IS the intended momentum-restoring design difference that zpitch(Zeff) carries.
+  n_main vs n_e·Zeff is a deliberate PENTRC-fidelity choice, exact in the z=1 main-ion domain.
+- 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.

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

docs/development/hdf5-conventions.md

@@ -17,7 +17,7 @@ These rules govern `gpec.h5` (and any future GPEC-produced HDF5 output); harness
 
 - **Groups are CamelCase at every level** (`ForceFreeStates/`, `PerSurface/`, `GalerkinIntegration/`).
 - **Datasets (leaves) are snake_case** (`eigenmode_energies`, `delta_prime_matrix`). Established physics symbols keep their natural case (`E`, `F`, `Q_root`, `pest3_Delta`, `2piF`).
-- **Data-driven tokens are stored verbatim**: coil-set names under `Input/RawInputs/Coils/`, KineticForces method tokens (`fgar`, …), scan indices (`Surface_<k>`, `psi_<i>`).
+- **Data-driven tokens are stored verbatim**: coil-set names under `Input/RawInputs/Coils/`, KineticForces method tokens (`fgar`, …), NTV species labels under `KineticForces/PerSpecies/` (`ion_z1_m2`, `impurity_z6_m12`, `electron` — charge and mass identify the species, with a numeric discriminator appended only if a run repeats a `(z, m)` pair), scan indices (`Surface_<k>`, `psi_<i>`).
 - **Word-valued names and boolean flags**: multi-word dataset names are snake_case English (`resonance_psi`, `trajectory_offsets`, `layer_widths`), never CamelCase — CamelCase is reserved for groups. A boolean flag is named for the state it asserts when true, with an `is_` prefix only where the bare word would read as a noun or collide with a data family: `is_rational` (bare `rational` would clash with the `rational_*` coordinate family) versus `enabled`, `truncated`, `no_root`, which already read as predicates.
 - **Literature capitalization for physics symbols**: names match the standard literature — `D_I`, `D_R`, `Delta_prime`, `tau_R`, `tau_A` (not `di`, `dr`, `delta_prime`, `taur`); lowercase stays where the literature is lowercase (`alpha`, `q`, `beta*`, `delta_s`).
 - **Scalar equilibrium parameters spell out the physics**: `R_axis`, `Z_axis`, `B_T_axis`, `a_mean`, `aspect_ratio`, `I_p`, `q_edge`, `beta_N`, `delta_upper`/`delta_lower`. The qualifier is a trailing subscript (`_axis`, `_edge`, `_wall`, `_min`, `_max`, `_upper`, `_lower`, `_extremum`), and a numbered literature definition keeps its number as the last subscript (`beta_p_1`, `l_i_2`). Fortran-era contractions (`bt0`, `amean`, `crnt`, `qa`, `betan`, `li1`) survive only as `Equilibrium.EquilibriumParameters` struct fields: `EQUIL_H5_NAMES` in `src/HDF5Schema.jl` maps each field to its dataset name, and `EQUIL_H5_SKIP` drops fields that duplicate another dataset or echo a control flag.
@@ -44,7 +44,7 @@ Top level (10 groups):
 | `LocalStability/` | Mercier `D_I`, resistive interchange `D_R`, `ballooning_Delta_prime` on `psi`; the ballooning α boundary on `ballooning_psi` |
 | `SingularSurfaces/` | Per-rational-surface data: `rational_psi`/`rational_q`/`rational_m`/`rational_n`, GGJ coefficients, `Delta_prime_matrix`/`Delta_prime_raw`/`Delta_coil`/`pest3_A`/`pest3_B`/`pest3_Gamma` (Riccati or Galerkin alike), `Kinetic/` |
 | `PerturbedEquilibrium/` | `ForcingModes/`, `Response/`, `ResponseMatrices/`, `SingularCoupling/`, `Energies/`, control-surface spectra |
-| `KineticForces/` | `<method>/` (torque/energy profiles, `EnergyIntegrals/`, `KineticMatrices/`) |
+| `KineticForces/` | `<method>/` (torque/energy profiles, `EnergyIntegrals/`, `KineticMatrices/`); multi-ion runs add `PerSpecies/<species>/<method>/` with the same per-method layout, summing to the top-level total |
 | `Tearing/` | `PerSurface/` (+ `DpMatrix/`), `Roots/`, `LayerWidths/`, `Diagnostics/{ValidRoots,Poles,FilteredRoots}`, `Scan/Surface_<k>/` |
 | `SurfaceGeometries/` | `{Plasma,Wall}/{x,y,z}` point clouds |
 

docs/src/kinetic_forces.md

@@ -16,7 +16,10 @@ extension:
     `chi_phi` (toroidal momentum diffusivity ``\chi_\phi``). Each dataset
     carries a `units` attribute; the root carries `schema_version` and
     `provenance`. Densities are m⁻³, temperatures eV, frequencies rad/s,
-    diffusivities m²/s. Write files with `write_kinetic_h5`.
+    diffusivities m²/s. Additional datasets named `n_*` (e.g. `n_D`, `n_T`)
+    are named per-species densities for multi-ion runs; names outside `n_*`
+    are reserved for future schema fields and ignored. Write files with
+    `write_kinetic_h5` (which round-trips the per-species densities).
   - **ASCII (`.gpeckf`/`.kin`/`.dat`)** — legacy six-column whitespace table
     `psi_n n_i n_e T_i[eV] T_e[eV] omega_E`, retained for backward
     compatibility. Header rows are skipped.
@@ -25,6 +28,52 @@ The NTV calculation consumes `n_i, n_e, T_i, T_e, omega_E`; `chi_e`/`chi_phi`,
 when present, are carried for the resistive-layer (SLAYER) analysis and ignored
 here.
 
+## Multi-ion runs
+
+A plasma may declare an arbitrary list of main-ion species; the NTV is computed
+per species under one shared full-composition ``Z_\mathrm{eff}`` and summed
+(``\tau = \sum_s \tau_s``) over the main ions, the quasineutrality-closing
+impurity, and (with `electron = true`) the electrons. This applies to both NTV
+paths: the post-PE ψ-quadrature diagnostic and the self-consistent
+`kinetic_source = "calculated"` matrices.
+
+```toml
+[KineticForces]
+kinetic_file = "kinetic.h5"    # n_i column/dataset = TOTAL main-ion density
+electron = true                # add electron NTV in addition to the ion species
+zimp = 6                       # impurity charge (closes quasineutrality)
+mimp = 12                      # impurity mass
+
+[[KineticForces.ion_species]]
+z = 1
+m = 2
+fraction = 0.5                 # this species' share of the total n_i
+[[KineticForces.ion_species]]
+z = 1
+m = 3
+density = "n_T"                # or: explicit n_* dataset from the HDF5 kinetic file
+```
+
+Each species sets exactly one of `fraction` (share of the file's total `n_i`)
+or `density` (a named `n_*` dataset). An all-fraction list must sum to 1; in a
+mixed list the impurity content is set by the file's `n_i`/`n_e` deficit, not
+by a fraction shortfall, and fractions may sum below (never above) 1. Every
+main-ion charge must satisfy `z < zimp`. An empty `ion_species` list runs the
+single main ion from `zi`/`mi` (with `electron = true` still adding the
+electron species — the electron flag always means *in addition to* the ions).
+
+The summed total is written to `KineticForces/<method>/` exactly as in a single-ion
+run, and each species' own contribution to `KineticForces/PerSpecies/<label>/<method>/`
+— labels are `ion_z<z>_m<m>`, `impurity_z<z>_m<m>`, and `electron` (a numeric suffix
+is appended only if a run repeats a `(z, m)` pair). Every per-species group carries the
+same datasets and metadata as the total, so `dTdpsi`, `T`, and `EnergyIntegrals/` are
+available per species. The summed cumulative torque profile is a diagnostic (linear
+interpolation onto the union grid); the summed `total_torque` scalar is the exact
+Gauss-Kronrod value.
+
+The profile-scaling knobs below are not supported together with a multi-ion
+`ion_species` list (they error).
+
 ## Profile Scaling Knobs
 
 Seven scaling factors are available on `KineticForcesControl` to modify kinetic

examples/Solovev_kinetic_calculated_example/gpec.toml

@@ -1,7 +1,9 @@
 # Solovev analytical equilibrium — n=1 ideal stability with calculated kinetic matrices.
-# Exercises the real KineticForces (NTV) physics path at full strength; kinetic regression
-# fixture (et[1]≈1.894-1.525i). The equilibrium is generated analytically from the embedded
-# [SOL_INPUT] section, with n·T tied to the Solovev P₀≈4.27e4 Pa via the local kinetic.dat.
+# Exercises the real KineticForces (NTV) physics path at full strength, with the main ion
+# declared as a 50/50 D-T mix so the self-consistent kinetic matrices are accumulated over
+# the multi-ion species loop (pure plasma, n_i = n_e, so no impurity species arises). The
+# equilibrium is generated analytically from the embedded [SOL_INPUT] section, with n·T tied
+# to the Solovev P₀≈4.27e4 Pa via the local kinetic.dat.
 
 [Equilibrium]
 eq_type = "sol"                         # Type of the input 2D equilibrium file
@@ -49,6 +51,19 @@ ucrit = 1e3                   # Maximum fraction of solutions allowed before re-
 write_outputs_to_HDF5 = true  # Write outputs to the HDF5 file
 verbose = false               # Enable verbose logging
 
+[KineticForces]
+kinetic_file = "kinetic.dat"  # Kinetic profile file: psi_n, n_i, n_e, T_i, T_e, omega_E columns
+
+[[KineticForces.ion_species]]
+z = 1                         # Species charge (fundamental units)
+m = 2                         # Species mass (proton masses) — Deuterium
+fraction = 0.5                # This species' share of the total main-ion density n_i
+
+[[KineticForces.ion_species]]
+z = 1                         # Species charge (fundamental units)
+m = 3                         # Species mass (proton masses) — Tritium
+fraction = 0.5                # This species' share of the total main-ion density n_i
+
 # Solovev analytic equilibrium parameters (eq_type = "sol"); see SolovevConfig in src/Equilibrium.
 [SOL_INPUT]
 mr = 128       # Number of radial grid zones

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