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

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

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

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

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

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

examples/Solovev_kinetic_multiion_example/run_example.jl

@@ -0,0 +1,4 @@
+using Pkg;
+Pkg.activate(joinpath(@__DIR__, "../.."))
+using GeneralizedPerturbedEquilibrium
+GeneralizedPerturbedEquilibrium.main([dirname(@__FILE__)])

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