From 8d985115ad1ac3164ffe33113ce6d16afd7f2402 Mon Sep 17 00:00:00 2001 From: logan-nc Date: Sat, 24 Oct 2026 14:43:44 -0400 Subject: [PATCH 1/3] ErrorFields - BUGFIX! - Sum the NTV couplings over the resolved kinetic species The error-field NTV couplings computed every correction array's torque for the deck's main ion alone, while the KineticForces diagnostics of the same run summed over the resolved species. The run's species list now reaches efc_couplings, whose torque at every rotation shift is the same species sum with each species' own profiles shifted, and the per-species control construction is shared by both paths. Single-ion decks are unchanged; the reference rotation stays the full-composition main ion's. Co-Authored-By: Claude Fable 5.1 Claude-Session: https://claude.ai/code/session_0186CLfjeJeSusNTva5RMfGU --- docs/src/error_fields.md | 4 +- src/GeneralizedPerturbedEquilibrium.jl | 62 ++++++++++++++++---------- src/KineticForces/Compute.jl | 11 +++++ test/runtests_error_fields.jl | 18 ++++++++ test/runtests_multiion.jl | 9 ++++ 5 files changed, 79 insertions(+), 25 deletions(-) diff --git a/docs/src/error_fields.md b/docs/src/error_fields.md index 1aa1e337f5..9d5e41bbba 100644 --- a/docs/src/error_fields.md +++ b/docs/src/error_fields.md @@ -362,7 +362,9 @@ published zero-rotation limits, `C_c √(T_0/T_residual)` (SPARC, `residual_only evaluates each one's `C_c`, resonant fraction, and NTV torque per kAt² for its whole field and for its field with the dominant mode projected out — two plasma-response evaluations of the unit-current spectrum followed by the kinetic torque, which needs a `[KineticForces]` section — -and writes `ErrorFields/NTV/`. The spectrum is evaluated on the run's `[ForcingTerms]` boundary +and writes `ErrorFields/NTV/`. With `ion_species` or `electron = true` in that section, the torques +are the same species sum as the NTV diagnostics, and the reference rotation is the main ion's. The +spectrum is evaluated on the run's `[ForcingTerms]` boundary grid and normalized by the magnitude of the array's ampere-turns, `|nw| × max|I|`, so the winding sense and current pattern of the deck stay the phase reference; the torques are stored with their sign, and the limits consume the budget with their magnitude (the sign depends on the rotation diff --git a/src/GeneralizedPerturbedEquilibrium.jl b/src/GeneralizedPerturbedEquilibrium.jl index db5e75cfc4..afe4b85c0a 100755 --- a/src/GeneralizedPerturbedEquilibrium.jl +++ b/src/GeneralizedPerturbedEquilibrium.jl @@ -332,7 +332,7 @@ function main_from_inputs( end ef_start = time() - error_fields = run_error_fields(inputs, ffs_result, pe_state, preloaded_coil_sets, kf_ctrl, kinetic_profiles) + error_fields = run_error_fields(inputs, ffs_result, pe_state, preloaded_coil_sets, kf_ctrl, kinetic_profiles, kf_species) coil_sensitivities, monte_carlo, locking_risk, efc_couplings = error_fields === nothing ? (nothing, nothing, nothing, nothing) : error_fields if "ErrorFields" in keys(inputs) push!(runtimes, "error_fields" => time() - ef_start) @@ -1108,13 +1108,8 @@ function run_kinetic_forces( states = KineticForces.KineticForcesState[] labels = String[] for sp in species - # 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))..., - ion_charge=sp.z, ion_mass=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; + KineticForces.compute_torque_all_methods!(st, kf_intr, KineticForces.species_control(kf_ctrl, sp), result.equil, sp.profiles; axis_psi_c=axis_psi_c) push!(states, st) push!(labels, sp.label) @@ -1154,7 +1149,7 @@ function forcing_terms_control(inputs::Dict{String,Any}) end """ - run_error_fields(inputs, result, pe_state, preloaded_coil_sets, kf_ctrl, kinetic_profiles) -> (sensitivities, monte_carlo, risk, couplings) or nothing + run_error_fields(inputs, result, pe_state, preloaded_coil_sets, kf_ctrl, kinetic_profiles, species=nothing) -> (sensitivities, monte_carlo, risk, couplings) or nothing Linearize every coil set's resonant drive with respect to its rigid shifts and tilts and write `ErrorFields/CoilSensitivities/` when the deck carries an `[ErrorFields]` section. When the @@ -1163,8 +1158,9 @@ Carlo on the full-window dominant mode with the `[ErrorFields.MonteCarlo]` setti `ErrorFields/MonteCarlo/`; with an `[ErrorFields.scenario]` table as well, evaluate the locking risk (and the tolerance scan when `[ErrorFields.Risk]` names `scan_scales`) and write `ErrorFields/Risk/`; with `[ErrorFields.NTV]` naming correction arrays, evaluate their overlap and NTV -torque couplings per kilo-ampere-turn (needs the kinetic context) and write `ErrorFields/NTV/`. -Stages not requested return `nothing`. Needs the +torque couplings per kilo-ampere-turn (needs the kinetic context) and write `ErrorFields/NTV/`; +`species`, the run's resolved NTV species list, makes those torques the same species sum as the +KineticForces diagnostics. Stages not requested return `nothing`. Needs the perturbed-equilibrium state's singular-coupling matrix and coil-format forcing, and errors otherwise: a deck asking for error-field sensitivities without them is a misconfiguration, not a case to skip silently. Coil geometry is rebuilt from the deck unless a replay injected it. @@ -1175,7 +1171,8 @@ function run_error_fields( pe_state, preloaded_coil_sets::Union{Nothing,Vector{ForcingTerms.CoilSet}}, kf_ctrl::KineticForces.KineticForcesControl=KineticForces.KineticForcesControl(), - kinetic_profiles=nothing + kinetic_profiles=nothing, + species=nothing ) ("ErrorFields" in keys(inputs)) || return nothing @@ -1255,7 +1252,7 @@ function run_error_fields( missing = setdiff(ntv_ctrl.efc_coils, [cs.name for cs in efc_sets]) isempty(missing) || error("[ErrorFields.NTV] efc_coils not among the run's coil sets: $(join(missing, ", "))") ntv_start = time() - couplings = efc_couplings(result, efc_sets, rc, dom, cfg, kf_ctrl, kinetic_profiles; method=ntv_ctrl.method, verbose=ntv_ctrl.verbose, + couplings = efc_couplings(result, efc_sets, rc, dom, cfg, kf_ctrl, kinetic_profiles; species, method=ntv_ctrl.method, verbose=ntv_ctrl.verbose, rotation_scan=ntv_ctrl.rotation_scan, scan_points=ntv_ctrl.rotation_scan_points, scan_max_points=ntv_ctrl.rotation_scan_max_points, scan_tolerance=ntv_ctrl.rotation_scan_tolerance, span_factor=ntv_ctrl.rotation_span_factor, offset_factor=ntv_ctrl.rotation_offset_factor) @info "NTV couplings of $(length(couplings)) correction arrays in $(@sprintf("%.1f", time() - ntv_start)) s" @@ -1280,7 +1277,7 @@ function run_error_fields( end """ - efc_couplings(ffs, coil_sets, rc, dom, cfg, kf_ctrl, kinetic_profiles; mode=1, method="fgar", verbose=false) -> Vector{ErrorFields.EFCCoupling} + efc_couplings(ffs, coil_sets, rc, dom, cfg, kf_ctrl, kinetic_profiles; species=nothing, mode=1, method="fgar", verbose=false) -> Vector{ErrorFields.EFCCoupling} The couplings of each correction coil array in `coil_sets` per kilo-ampere-turn: its dominant-mode overlap and resonant fraction from its spectrum on the control surface, and its @@ -1292,7 +1289,12 @@ run's `[ForcingTerms]` coil configuration, so the boundary grid matches the one sensitivities were evaluated on; `rc` must be the coupling of `ffs`'s own solve and `kinetic_profiles` its kinetic context. The spectrum is normalized by the magnitude of the array's ampere-turns, `|nw| × max|I|`, keeping its winding sense and current pattern as the -phase reference; the torque is stored with its sign. +phase reference; the torque is stored with its sign. With `species`, the run's resolved NTV +species list, each torque is the sum over the species of their own torques — each with the +control of its identity alone and its own profiles, shifted in the rotation scan — exactly as +the KineticForces diagnostics sum them; the reference rotation stays that of the full-composition +main ion in `kinetic_profiles`. A species sum's cumulative profile is re-integrated from the +summed torque density, so its end point matches the summed total only to that integration's accuracy. """ function efc_couplings( ffs::ForceFreeStatesResult, @@ -1302,6 +1304,7 @@ function efc_couplings( cfg::ForcingTerms.CoilConfig, kf_ctrl::KineticForces.KineticForcesControl, kinetic_profiles; + species=nothing, mode::Int=1, method::AbstractString="fgar", rotation_scan::Bool=true, @@ -1320,8 +1323,8 @@ function efc_couplings( v = dom.right_singular_vectors[:, mode] quiet = PerturbedEquilibrium.PerturbedEquilibriumControl(; compute_response=true, compute_singular_coupling=false, verbose=verbose, write_outputs_to_HDF5=false) - # 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. + # Reference rotation and scan span of the full-composition main ion in the kinetic profiles, with + # the torque kernel's own diamagnetic frequencies (Logan & Park 2013 Eq. 7): ω_φ = ω_E + ω_*n + ω_*T. kp = kinetic_profiles chi1 = 2π * equil.psio chrg = kf_ctrl.ion_charge * Utilities.E_CHG @@ -1348,12 +1351,24 @@ function efc_couplings( KineticForces.set_perturbation_data!(kf_intr, pe_state, ffs, equil, ffs.metric) return kf_intr end - # 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() + # Torque and cumulative torque profile (resampled onto the kinetic ψ grid) of one perturbation at one + # rotation shift Δ, summed over the species when a list is given. + shifted(p, Δ) = Δ == 0 ? p : Equilibrium.shift_exb_rotation(p, Δ) + function torque_of(kf_intr, Δ) logger = verbose ? Base.CoreLogging.current_logger() : Base.CoreLogging.SimpleLogger(stderr, Base.CoreLogging.Warn) - Base.CoreLogging.with_logger(logger) do - KineticForces.compute_torque_all_methods!(kf_state, kf_intr, kf_ctrl, equil, profiles) + kf_state = Base.CoreLogging.with_logger(logger) do + if species === nothing + st = KineticForces.KineticForcesState() + KineticForces.compute_torque_all_methods!(st, kf_intr, kf_ctrl, equil, shifted(kp, Δ)) + st + else + states = map(species) do sp + st = KineticForces.KineticForcesState() + KineticForces.compute_torque_all_methods!(st, kf_intr, KineticForces.species_control(kf_ctrl, sp), equil, shifted(sp.profiles, Δ)) + st + end + KineticForces.combine_species_states(states) + end end r = kf_state.method_results[String(method)] T = real(r.total_torque) @@ -1393,9 +1408,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)))" diff --git a/src/KineticForces/Compute.jl b/src/KineticForces/Compute.jl index a4f8a4ec8c..5c7c24efde 100644 --- a/src/KineticForces/Compute.jl +++ b/src/KineticForces/Compute.jl @@ -263,6 +263,17 @@ function _accumulate_interp!(acc, grid, xs::AbstractVector{Float64}, ys::Abstrac return acc end +""" + species_control(ctrl::KineticForcesControl, sp) -> KineticForcesControl + +The control one species of the NTV sum is computed with: `ctrl` with only the identity of `sp` +(charge, mass, electron flag) in place of the deck's main ion and no species list of its own. +`KineticForcesControl` is immutable, so every species gets a fresh copy. +""" +species_control(ctrl::KineticForcesControl, sp) = + KineticForcesControl(; (f => getfield(ctrl, f) for f in fieldnames(KineticForcesControl))..., + ion_charge=sp.z, ion_mass=sp.m, electron=sp.electron, ion_species=IonSpecies[]) + """ combine_species_states(states) -> KineticForcesState diff --git a/test/runtests_error_fields.jl b/test/runtests_error_fields.jl index ccd9e098e3..e562ce4606 100644 --- a/test/runtests_error_fields.jl +++ b/test/runtests_error_fields.jl @@ -533,6 +533,24 @@ include("h5_metadata_check.jl") @test with_profiles[3].series_list[2][:y] == c.torque_residual_profile[:, 2] @test_throws ErrorException GPEC.efc_couplings(ffs, sets, rc, dom, cfg, GPEC.KineticForces.KineticForcesControl(), nothing) + # Species-resolved couplings: the deck's main ion resolved on its own reproduces the run's + # single-ion torque, and a species 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"])...) + kfile = joinpath(dir, "kinetic.dat") + kp = EQ.load_kinetic_profiles(kfile; zi=kf_ctrl.ion_charge, zimp=kf_ctrl.impurity_charge, mi=kf_ctrl.ion_mass, mimp=kf_ctrl.impurity_mass, + chi1=2π * ffs.equil.psio) + sp = EQ.resolve_ntv_species(kfile, [KF.IonSpecies(; z=kf_ctrl.ion_charge, m=kf_ctrl.ion_mass, fraction=1.0)]; + electron=false, zimp=kf_ctrl.impurity_charge, mimp=kf_ctrl.impurity_mass) + @test length(sp) == 2 && !sp[1].electron && sp[2].z == kf_ctrl.impurity_charge + coupling_of(species) = GPEC.efc_couplings(ffs, sets[3:3], rc, dom, cfg, kf_ctrl, kp; species, rotation_scan=false)[1] + ion, imp, both = coupling_of(sp[1:1]), coupling_of(sp[2:2]), coupling_of(sp) + @test ion.torque_full_per_kat2 ≈ c.torque_full_per_kat2 rtol = 1e-6 + @test ion.torque_residual_per_kat2 ≈ c.torque_residual_per_kat2 rtol = 1e-6 + @test both.torque_full_per_kat2 ≈ ion.torque_full_per_kat2 + imp.torque_full_per_kat2 rtol = 1e-12 + @test both.torque_residual_per_kat2 ≈ ion.torque_residual_per_kat2 + imp.torque_residual_per_kat2 rtol = 1e-12 + @test !EF.has_rotation_scan(both) && both.torque_full_profile[1, 1] == 0.0 + # Central differences: doubling the step moves the derivatives at O(h²). coarse = EF.compute_coil_sensitivities(sets, 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) diff --git a/test/runtests_multiion.jl b/test/runtests_multiion.jl index 95b0d58eee..c440c8efa3 100644 --- a/test/runtests_multiion.jl +++ b/test/runtests_multiion.jl @@ -33,6 +33,15 @@ using HDF5 @test isempty(KF.KineticForcesControl().ion_species) # default single-ion end + @testset "species_control" begin + # One species' control is the run's control with only its identity overridden and no species list. + ctrl = KF.KineticForcesControl(; ion_species=[KF.IonSpecies(; z=1, m=2, fraction=0.5), KF.IonSpecies(; z=1, m=3, fraction=0.5)], + electron=true, rtol_pitch=1e-3, fgar_flag=true) + sc = KF.species_control(ctrl, (; z=2, m=4, electron=false)) + @test (sc.ion_charge, sc.ion_mass, sc.electron) == (2, 4, false) && isempty(sc.ion_species) + @test all(getfield(sc, f) == getfield(ctrl, f) for f in fieldnames(KF.KineticForcesControl) if f ∉ (:ion_charge, :ion_mass, :electron, :ion_species)) + end + @testset "resolve_ntv_species" begin # Synthetic ASCII kinetic profile: psi, n_i(total), n_e, T_i, T_e, omega_E. psi = collect(0.0:0.05:1.0) From 787cf6acfbf9c8cce8cfda209740b6dfd0413726 Mon Sep 17 00:00:00 2001 From: logan-nc Date: Sat, 24 Oct 2026 14:51:52 -0400 Subject: [PATCH 2/3] KineticForces - TEST - Compare the species control's NaN tolerance defaults with isequal Co-Authored-By: Claude Fable 5.1 Claude-Session: https://claude.ai/code/session_0186CLfjeJeSusNTva5RMfGU --- test/runtests_multiion.jl | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/test/runtests_multiion.jl b/test/runtests_multiion.jl index c440c8efa3..9cfa4bc7c6 100644 --- a/test/runtests_multiion.jl +++ b/test/runtests_multiion.jl @@ -39,7 +39,8 @@ using HDF5 electron=true, rtol_pitch=1e-3, fgar_flag=true) sc = KF.species_control(ctrl, (; z=2, m=4, electron=false)) @test (sc.ion_charge, sc.ion_mass, sc.electron) == (2, 4, false) && isempty(sc.ion_species) - @test all(getfield(sc, f) == getfield(ctrl, f) for f in fieldnames(KF.KineticForcesControl) if f ∉ (:ion_charge, :ion_mass, :electron, :ion_species)) + # isequal: the energy-tolerance overrides default to NaN. + @test all(isequal(getfield(sc, f), getfield(ctrl, f)) for f in fieldnames(KF.KineticForcesControl) if f ∉ (:ion_charge, :ion_mass, :electron, :ion_species)) end @testset "resolve_ntv_species" begin From a305d61e5695f5eb40537ee080044230ea46d15a Mon Sep 17 00:00:00 2001 From: logan-nc Date: Sat, 24 Oct 2026 15:08:45 -0400 Subject: [PATCH 3/3] ErrorFields - REFACTOR - Factor the per-shift NTV torque of the couplings into ntv_torque_state The species sum at one rotation shift becomes a module-level function the couplings call, so the test can check it on the run's own perturbation instead of re-solving the plasma response. Co-Authored-By: Claude Fable 5.1 Claude-Session: https://claude.ai/code/session_0186CLfjeJeSusNTva5RMfGU --- src/GeneralizedPerturbedEquilibrium.jl | 42 +++++++++++++++++--------- test/runtests_error_fields.jl | 26 ++++++++-------- 2 files changed, 41 insertions(+), 27 deletions(-) diff --git a/src/GeneralizedPerturbedEquilibrium.jl b/src/GeneralizedPerturbedEquilibrium.jl index afe4b85c0a..c41b0894a0 100755 --- a/src/GeneralizedPerturbedEquilibrium.jl +++ b/src/GeneralizedPerturbedEquilibrium.jl @@ -1276,6 +1276,31 @@ function run_error_fields( return sens, monte_carlo, risk, couplings 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 deck's single ion on `kinetic_profiles` when `species` is `nothing`, otherwise the sum +over the resolved species, each with the control of its identity alone and its own profiles +shifted, as `run_kinetic_forces` sums them. Unshifted profiles are used as they are, so `Δω = 0` +is the same evaluation as the torque diagnostics. `axis_psi_c` is the near-axis validity boundary +passed to the torque. +""" +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, KineticForces.species_control(kf_ctrl, sp), equil, shifted(sp.profiles); axis_psi_c) + st + end + return KineticForces.combine_species_states(states) +end + """ efc_couplings(ffs, coil_sets, rc, dom, cfg, kf_ctrl, kinetic_profiles; species=nothing, mode=1, method="fgar", verbose=false) -> Vector{ErrorFields.EFCCoupling} @@ -1351,24 +1376,11 @@ function efc_couplings( KineticForces.set_perturbation_data!(kf_intr, pe_state, ffs, equil, ffs.metric) return kf_intr end - # Torque and cumulative torque profile (resampled onto the kinetic ψ grid) of one perturbation at one - # rotation shift Δ, summed over the species when a list is given. - shifted(p, Δ) = Δ == 0 ? p : Equilibrium.shift_exb_rotation(p, Δ) + # Torque and cumulative torque profile (resampled onto the kinetic ψ grid) of one perturbation at one rotation shift. function torque_of(kf_intr, Δ) logger = verbose ? Base.CoreLogging.current_logger() : Base.CoreLogging.SimpleLogger(stderr, Base.CoreLogging.Warn) kf_state = Base.CoreLogging.with_logger(logger) do - if species === nothing - st = KineticForces.KineticForcesState() - KineticForces.compute_torque_all_methods!(st, kf_intr, kf_ctrl, equil, shifted(kp, Δ)) - st - else - states = map(species) do sp - st = KineticForces.KineticForcesState() - KineticForces.compute_torque_all_methods!(st, kf_intr, KineticForces.species_control(kf_ctrl, sp), equil, shifted(sp.profiles, Δ)) - st - end - KineticForces.combine_species_states(states) - end + ntv_torque_state(kf_intr, kf_ctrl, equil, kp, species, Δ) end r = kf_state.method_results[String(method)] T = real(r.total_torque) diff --git a/test/runtests_error_fields.jl b/test/runtests_error_fields.jl index e562ce4606..680362be89 100644 --- a/test/runtests_error_fields.jl +++ b/test/runtests_error_fields.jl @@ -533,23 +533,25 @@ include("h5_metadata_check.jl") @test with_profiles[3].series_list[2][:y] == c.torque_residual_profile[:, 2] @test_throws ErrorException GPEC.efc_couplings(ffs, sets, rc, dom, cfg, GPEC.KineticForces.KineticForcesControl(), nothing) - # Species-resolved couplings: the deck's main ion resolved on its own reproduces the run's - # single-ion torque, and a species list's torque is the sum of its members' own torques. + # Species-resolved torque of the run's own perturbation, on the diagnostics' axis boundary: the + # deck's 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"])...) kfile = joinpath(dir, "kinetic.dat") kp = EQ.load_kinetic_profiles(kfile; zi=kf_ctrl.ion_charge, zimp=kf_ctrl.impurity_charge, mi=kf_ctrl.ion_mass, mimp=kf_ctrl.impurity_mass, chi1=2π * ffs.equil.psio) - sp = EQ.resolve_ntv_species(kfile, [KF.IonSpecies(; z=kf_ctrl.ion_charge, m=kf_ctrl.ion_mass, fraction=1.0)]; - electron=false, zimp=kf_ctrl.impurity_charge, mimp=kf_ctrl.impurity_mass) - @test length(sp) == 2 && !sp[1].electron && sp[2].z == kf_ctrl.impurity_charge - coupling_of(species) = GPEC.efc_couplings(ffs, sets[3:3], rc, dom, cfg, kf_ctrl, kp; species, rotation_scan=false)[1] - ion, imp, both = coupling_of(sp[1:1]), coupling_of(sp[2:2]), coupling_of(sp) - @test ion.torque_full_per_kat2 ≈ c.torque_full_per_kat2 rtol = 1e-6 - @test ion.torque_residual_per_kat2 ≈ c.torque_residual_per_kat2 rtol = 1e-6 - @test both.torque_full_per_kat2 ≈ ion.torque_full_per_kat2 + imp.torque_full_per_kat2 rtol = 1e-12 - @test both.torque_residual_per_kat2 ≈ ion.torque_residual_per_kat2 + imp.torque_residual_per_kat2 rtol = 1e-12 - @test !EF.has_rotation_scan(both) && both.torque_full_profile[1, 1] == 0.0 + resolve(list) = EQ.resolve_ntv_species(kfile, list; electron=false, zimp=kf_ctrl.impurity_charge, mimp=kf_ctrl.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) = GPEC.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 + @test torque(resolve([KF.IonSpecies(; z=kf_ctrl.ion_charge, m=kf_ctrl.ion_mass, fraction=1.0)])) ≈ single rtol = 1e-6 + dt = resolve([KF.IonSpecies(; z=kf_ctrl.ion_charge, m=2, fraction=0.5), KF.IonSpecies(; z=kf_ctrl.ion_charge, 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, rc, ffs.equil, cfg,