diff --git a/docs/src/error_fields.md b/docs/src/error_fields.md index c108780fd6..1cd5c3b47c 100644 --- a/docs/src/error_fields.md +++ b/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]`: diff --git a/src/ErrorFields/NTVLimits.jl b/src/ErrorFields/NTVLimits.jl index 4504a8b457..21d700914d 100644 --- a/src/ErrorFields/NTVLimits.jl +++ b/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)))" diff --git a/test/runtests_error_fields.jl b/test/runtests_error_fields.jl index 35f01a3c70..73b78fdfc1 100644 --- a/test/runtests_error_fields.jl +++ b/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)