KineticForces - REFACTOR - Restore the KineticForcesControl freeze and drop dead deck keys

A develop merge into this branch silently reverted PR #367's immutability
freeze on KineticForcesControl, because the multi-species loop mutated a
deepcopy of the run control. That made the branch carry an undocumented revert
of an upstream refactor and depend on the mutability upstream removed
("KineticForcesControl had no writers in src/ at all" — this loop was exactly
such a writer, and the only one left).

Build each species' control with the @kwdef keyword constructor instead, the
idiom the freeze commit documents and that check_psi_quadrature_convergence
already uses, and restore @kwdef struct. Threading species identity through
explicit arguments was rejected: ctrl is passed whole into integrate_psi_quadgk
and onward, so that would touch several nested signatures.

Also drop six keys from the multi-ion deck that develop removed when it
replaced the integrator-selection axis with integrator= (mat_flag, ode_flag,
force_wv_symmetry, use_parallel, parallel_threads, populate_dense_xi). They
only emitted deprecation warnings, but populate_dense_xi in particular implied
a guarantee it no longer provides — had integrator= been dropped on that
belief, the run would have fallen back to the riccati default, which does not
populate dense xi, and PE/NTV would have gone quietly wrong. Document the
PerSpecies layout in the h5 schema table.

Co-Authored-By: Claude Fable 5 <noreply@anthropic.com>

Author: logan-nc

docs/development/hdf5-conventions.md

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

examples/Solovev_kinetic_multiion_example/gpec.toml

@@ -49,8 +49,6 @@ reg_spot = 0.05                         # Regularization width for singular surf
 [ForceFreeStates]
 integrator = "forward"        # Dense ξ profiles required by PerturbedEquilibrium and kinetic torque diagnostics
 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
@@ -69,13 +67,9 @@ kinetic_factor = 0.0          # Scaling of kinetic matrices (0 = ideal path; >0
 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
+# Integration-range truncation (see ForceFreeStatesControl docstring for details)
 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.

src/GeneralizedPerturbedEquilibrium.jl

@@ -954,11 +954,11 @@ function run_kinetic_forces(
             states = KineticForces.KineticForcesState[]
             labels = String[]
             for sp in species
-                sctrl = deepcopy(kf_ctrl)
-                sctrl.zi = sp.z
-                sctrl.mi = sp.m
-                sctrl.electron = sp.electron
-                sctrl.ion_species = KineticForces.IonSpecies[]
+                # KineticForcesControl is immutable, so each species gets a fresh control
+                # built from the run's control with only its own identity overridden.
+                sctrl = KineticForces.KineticForcesControl(;
+                    (f => getfield(kf_ctrl, f) for f in fieldnames(KineticForces.KineticForcesControl))...,
+                    zi=sp.z, mi=sp.m, electron=sp.electron, ion_species=KineticForces.IonSpecies[])
                 st = KineticForces.KineticForcesState()
                 KineticForces.compute_torque_all_methods!(st, kf_intr, sctrl, result.equil, sp.profiles)
                 push!(states, st)

src/KineticForces/KineticForcesStructs.jl

@@ -80,8 +80,12 @@ Constructed via keyword arguments or from a TOML dict:
 ```julia
 ctrl = KineticForcesControl(; (Symbol(k) => v for (k, v) in inputs["KineticForces"])...)
 ```
+
+Immutable: vary a field by building a new control rather than assigning to one (the
+multi-species loop does this per species, and `check_psi_quadrature_convergence`'s test
+builds a second control for its differing tolerance).
 """
-@kwdef mutable struct KineticForcesControl
+@kwdef struct KineticForcesControl
     # Moment type
     moment::String = "pressure"     # "heat" or "pressure"