ErrorFields - BUGFIX! - Sum the NTV couplings over the resolved kinetic species

docs/src/error_fields.md

@@ -90,6 +90,7 @@ rotation_offset_factor = 2.0
 The tilt pivot matters for multi-filament winding packs: `"conductor"` rotates each filament about its own centre, `"set"` rotates the pack rigidly about its common centre, which is what an engineering axis-line tolerance constrains.
 A coil set that is not an error-field source (a correction array or anything in `exclude_coils`) is swept like the others but left out of the as-designed field, the Monte Carlo and the scans; a tolerance that names it is an error.
 The NTV limits need a `[KineticForces]` section, since each array's torque per kAt² is one kinetic torque evaluation of its unit-current spectrum, with and without the dominant mode projected out.
+With several ion species or `electron = true`, that torque is the same species sum as the NTV diagnostics; the reference rotation is the main ion's.
 
 The tolerances are engineering data, reviewed independently of any run, so they live in their own file, echoed into `Input/RawInputs/ErrorFields/tolerance_toml_raw` for replay; every key is checked against the schema.
 Its tables are `[ErrorFields.defaults]` (fallbacks for every coil block), one `[[ErrorFields.coil]]` per coil set that moves on its own, `[[ErrorFields.coherent_group]]` for coil sets that move together, `[ErrorFields.correctability]` and `[ErrorFields.other_field]`:

src/ErrorFields/NTVLimits.jl

@@ -464,6 +464,29 @@ function efc_current_curve(c::EFCCoupling; delta_threshold::Real, torque_budget:
         max_correctable_overlap(c; kw(torque_budget)...)...)
 end
 
+"""
+    ntv_torque_state(kf_intr, kf_ctrl, equil, kinetic_profiles, species, Δω; axis_psi_c=0.0) -> KineticForcesState
+
+The NTV torque of the perturbation held by `kf_intr` with the E×B rotation rigidly shifted by `Δω`
+(rad/s): the main ion on `kinetic_profiles` when `species` is `nothing`, otherwise the sum over the
+resolved species (`KineticForces.ntv_species`), each on its own shifted profiles. `Δω = 0` uses the
+profiles unshifted.
+"""
+function ntv_torque_state(kf_intr, kf_ctrl::KineticForces.KineticForcesControl, equil, kinetic_profiles, species, Δω::Real; axis_psi_c::Float64=0.0)
+    shifted(p) = Δω == 0 ? p : Equilibrium.shift_exb_rotation(p, Δω)
+    if species === nothing
+        st = KineticForces.KineticForcesState()
+        KineticForces.compute_torque_all_methods!(st, kf_intr, kf_ctrl, equil, shifted(kinetic_profiles); axis_psi_c)
+        return st
+    end
+    states = map(species) do sp
+        st = KineticForces.KineticForcesState()
+        KineticForces.compute_torque_all_methods!(st, kf_intr, kf_ctrl, equil, shifted(sp.profiles); axis_psi_c, z=sp.z, m=sp.m, electron=sp.electron)
+        st
+    end
+    return KineticForces.combine_species_states(states)
+end
+
 """
     efc_couplings(resp, coil_sets, rc, dom, cfg, kf_ctrl; mode=1, method="fgar", rotation_scan=true, scan_points=9,
                   scan_max_points=21, scan_tolerance=0.05, span_factor=1.0, offset_factor=2.0,
@@ -474,8 +497,9 @@ dominant-mode overlap and resonant fraction from its spectrum on the control sur
 NTV torque (kinetic method `method`, with sign) for the whole field and for the field with mode `mode`
 of `dom` projected out, tabulated against the rotation shift when `rotation_scan`. `cfg` is the run's
 `[ForcingTerms]` coil configuration and `rc` must be the coupling of `resp`; the torque uses the kinetic
-profiles attached to `resp.ffs.equil`. The spectrum is normalized by `|nw| × max|I|`, keeping the
-array's winding sense and current pattern as the phase reference.
+profiles attached to `resp.ffs.equil`, summed over the resolved species when the run has several or
+`kf_ctrl.electron`; the reference rotation is the main ion's. The spectrum is normalized by `|nw| × max|I|`,
+keeping the array's winding sense and current pattern as the phase reference.
 """
 function efc_couplings(
     resp::PerturbedEquilibrium.PlasmaResponse,
@@ -506,6 +530,7 @@ function efc_couplings(
     # Reference rotation and scan span from the kinetic profiles, with the torque kernel's own
     # diamagnetic frequencies (Logan & Park 2013 Eq. 7): ω_φ = ω_E + ω_*n + ω_*T.
     kp = KineticForces.ntv_profiles(kf_ctrl, equil)
+    species = KineticForces.ntv_species(kf_ctrl, equil)
     chi1 = 2π * equil.psio
     chrg = KineticForces.main_ion(equil)[1] * Utilities.E_CHG
     ψk = kp.xs
@@ -532,11 +557,10 @@ function efc_couplings(
     end
     # One quiet logger for the whole scan, so a `maxlog` warning prints once per scan rather than once per point.
     logger = verbose ? Base.CoreLogging.current_logger() : Base.CoreLogging.SimpleLogger(stderr, Base.CoreLogging.Warn)
-    # Torque and cumulative torque profile (resampled onto the kinetic ψ grid) of one perturbation at one rotation.
-    function torque_of(kf_intr, profiles)
-        kf_state = KineticForces.KineticForcesState()
-        Base.CoreLogging.with_logger(logger) do
-            KineticForces.compute_torque_all_methods!(kf_state, kf_intr, kf_ctrl, equil, profiles)
+    # Torque and cumulative torque profile (resampled onto the kinetic ψ grid) of one perturbation at one rotation shift.
+    function torque_of(kf_intr, Δ)
+        kf_state = Base.CoreLogging.with_logger(logger) do
+            ntv_torque_state(kf_intr, kf_ctrl, equil, kp, species, Δ)
         end
         r = kf_state.method_results[String(method)]
         T = real(r.total_torque)
@@ -576,9 +600,8 @@ function efc_couplings(
         profiles_full = Dict{Float64,Vector{Float64}}()
         profiles_res = Dict{Float64,Vector{Float64}}()
         function torques(Δ)
-            profiles = Δ == 0 ? kp : Equilibrium.shift_exb_rotation(kp, Δ)
-            T_f, p_f = torque_of(kf_full, profiles)
-            T_r, p_r = torque_of(kf_res, profiles)
+            T_f, p_f = torque_of(kf_full, Δ)
+            T_r, p_r = torque_of(kf_res, Δ)
             profiles_full[Δ] = p_f
             profiles_res[Δ] = p_r
             verbose && @info "    $(cs.name) Δω = $(@sprintf("%+.3e", Δ)) rad/s: torque $(@sprintf("%.3e", T_f)) N·m per kAt² (residual $(@sprintf("%.3e", T_r)))"

test/runtests_error_fields.jl

@@ -568,6 +568,26 @@ include("display_check.jl")
             @test_throws "needs kinetic profiles" EF.efc_couplings(pe.response, sets, rc, dom, cfg, GPEC.KineticForces.KineticForcesControl())
             pe.response.ffs.equil.kinetic = attached
 
+            # Species-resolved torque of the run's own perturbation, on the diagnostics' axis boundary: the
+            # main ion resolved on its own reproduces the single-ion torque, and a D-T list's torque is the
+            # sum of its members' own torques.
+            KF, EQ = GPEC.KineticForces, GPEC.Equilibrium
+            kf_ctrl = KF.KineticForcesControl(; (Symbol(k) => v for (k, v) in inputs["KineticForces"])...)
+            kin = ffs.equil.kinetic
+            kp = KF.ntv_profiles(kf_ctrl, ffs.equil)
+            resolve(list) = EQ.resolve_ntv_species(kin, list; electron=false, zimp=kin.impurity_charge, mimp=kin.impurity_mass)
+            kf_intr = KF.KineticForcesInternal(ffs.equil; verbose=false)
+            KF.set_perturbation_data!(kf_intr, pe, ffs, ffs.equil, ffs.metric)
+            psi_c = KF.axis_validity_boundary(kf_ctrl, nothing, kp, ffs.equil, Float64[sng.psifac for sng in ffs.surfaces])
+            torque(species) = EF.ntv_torque_state(kf_intr, kf_ctrl, ffs.equil, kp, species, 0.0; axis_psi_c=psi_c).method_results["fgar"].total_torque
+            single = torque(nothing)
+            @test isfinite(single) && single != 0
+            zi, mi = KF.main_ion(ffs.equil)
+            @test torque(resolve([EQ.IonSpecies(; z=zi, m=mi, fraction=1.0)])) ≈ single rtol = 1e-6
+            dt = resolve([EQ.IonSpecies(; z=zi, m=2, fraction=0.5), EQ.IonSpecies(; z=zi, m=3, fraction=0.5)])
+            @test length(dt) >= 2 && !any(s.electron for s in dt)
+            @test torque(dt) ≈ sum(torque(dt[i:i]) for i in eachindex(dt)) rtol = 1e-12
+
             # Central differences: doubling the step moves the derivatives at O(h²).
             coarse = EF.compute_coil_sensitivities(sets[1:1], rc, ffs.equil, cfg,
                 EF.ErrorFieldsControl(; fd_step_shift_m=2e-3, fd_step_tilt_deg=0.2); psi=ffs.psilim, b_t0=sens.b_t0)