Repo - BUGFIX! - Fix numerical edge cases in q extrema, rational surfaces, error-field correction and the 2D vacuum modes

src/Equilibrium/Derived/GlobalParameters.jl

@@ -171,18 +171,14 @@ function equilibrium_qfind!(equil::PlasmaEquilibrium)
         c = q_d2(x0)
         d = q_d3(x0)
 
-        if d != 0.0
-            xcrit = -c / d
-            dx2 = xcrit^2 - 2b / d
-            if dx2 ≥ 0
-                dx = sqrt(dx2)
-                for delta in (dx, -dx)
-                    x = xcrit - delta
-                    if 0 ≤ x < xmax
-                        ψ = x0 + x
-                        push!(psiexl, ψ)
-                        push!(qexl, q_spline(ψ))
-                    end
+        disc = c^2 - 2b * d
+        if disc ≥ 0
+            s = -(c + copysign(sqrt(disc), c)) / 2
+            for x in (disc > 0 ? (2s / d, b / s) : (2s / d,))
+                if 0 ≤ x < xmax
+                    ψ = x0 + x
+                    push!(psiexl, ψ)
+                    push!(qexl, q_spline(ψ))
                 end
             end
         end

src/ErrorFields/Correction.jl

@@ -86,7 +86,7 @@ function correction_requirement(ctx::ResonantDriveContext, source::CoilOverlap,
     a = all(k -> isfinite(k) && k > 0, kat) ? kat : ones(length(kat))
     δ = [o.delta for o in arrays] ./ a
     f_mn = sum(abs2, δ) > 0 ? (-source.delta .* conj.(δ) ./ sum(abs2, δ)) ./ a : fill(NaN + NaN * im, length(δ))
-    f_ls = -((R ./ a') \ r_source) ./ a
+    f_ls = -(pinv(R ./ a') * r_source) ./ a
     r_ls = r_source .+ R * f_ls
     n_src = norm(source.spectrum)
     cosines = [(n = norm(o.spectrum); n > 0 && n_src > 0 ? abs(dot(o.spectrum, source.spectrum)) / (n * n_src) : NaN) for o in arrays]

src/ErrorFields/NTVLimits.jl

@@ -289,7 +289,7 @@ function correction_current(δ_ef::Real, c::EFCCoupling; delta_threshold::Real,
         a == 0 && return excess / c.delta_per_kat
         disc = c.delta_per_kat^2 - 4a * excess
         disc < 0 && return NaN
-        I = (c.delta_per_kat - sqrt(disc)) / (2a)
+        I = 2excess / (c.delta_per_kat + sqrt(disc))
         # The budget must not be exhausted at the root: a negative threshold is no correction (OMFIT masked these).
         return abs(c.torque_residual_per_kat2) * I^2 < torque_budget ? I : NaN
     end
@@ -416,7 +416,7 @@ function correction_current_upper(δ_ef::Real, c::EFCCoupling; delta_threshold::
         # Over-correction side: C·I − δ_ef = target·(1 − t_res·I²/T_0), i.e. a·I² + C·I − (δ_ef + target) = 0
         # with the same a as the under-correction quadratic; its positive root is the upper edge.
         a = target * t_res / torque_budget
-        upper = (-c.delta_per_kat + sqrt(c.delta_per_kat^2 + 4a * (δ_ef + target))) / (2a)
+        upper = 2(δ_ef + target) / (c.delta_per_kat + sqrt(c.delta_per_kat^2 + 4a * (δ_ef + target)))
         return min(upper, sqrt(torque_budget / t_res))
     end
     g(I) = c.delta_per_kat * I - δ_ef - target * threshold_factor(c, I; torque_budget, rotation_exponent, omega_reference, model)

src/ErrorFields/Risk.jl

@@ -436,7 +436,8 @@ end
 The tolerance multiplier at which the locking risk reaches `target_percent`, by linear
 interpolation in the logarithms of both the scale and the risk between the two scan points
 that bracket the target (the scan need not be monotonic; the first bracketing pair from the
-smallest scale is used). Returns `NaN` when no pair brackets the target. `corrected` reads the
+smallest scale is used). A pair that starts at scale 0 is interpolated linearly in the scale.
+Returns `NaN` when no pair brackets the target. `corrected` reads the
 error-field-corrected curve.
 """
 function allowable_tolerance(scan::ToleranceScan, target_percent::Real; corrected::Bool=false)
@@ -448,9 +449,11 @@ function allowable_tolerance(scan::ToleranceScan, target_percent::Real; correcte
         (p1 > 0 && p2 > 0) || continue
         l1, l2 = log10(p1), log10(p2)
         (min(l1, l2) <= lt <= max(l1, l2)) || continue
-        l1 == l2 && return scan.tolerance_scale[i]
+        s1, s2 = scan.tolerance_scale[i], scan.tolerance_scale[i+1]
+        lt == l1 && return s1
+        lt == l2 && return s2
         f = (lt - l1) / (l2 - l1)
-        return 10^(log10(scan.tolerance_scale[i]) + f * (log10(scan.tolerance_scale[i+1]) - log10(scan.tolerance_scale[i])))
+        return s1 > 0 ? 10^(log10(s1) + f * (log10(s2) - log10(s1))) : s1 + f * (s2 - s1)
     end
     return NaN
 end

src/ForceFreeStates/Surfaces/Finding.jl

@@ -17,6 +17,7 @@ function _find_rational_surfaces(equil::Equilibrium.PlasmaEquilibrium, nlow::Int
         # Loop over extrema of q, find all rational values in between
         for iex in 2:equil.params.mextrema
             dq = equil.params.qextrema_q[iex] - equil.params.qextrema_q[iex-1]
+            dq == 0 && continue
             m = trunc(Int, n * equil.params.qextrema_q[iex-1])
             if dq > 0
                 m += 1

src/Vacuum/Vacuum.jl

@@ -58,6 +58,8 @@ end
 
     mpert = length(inputs.m_modes)
     mlow = inputs.m_modes[1]
+    inputs.m_modes == mlow:(mlow+mpert-1) ||
+        throw(ArgumentError("2D vacuum response needs contiguous ascending poloidal modes, got $(inputs.m_modes)"))
     num_points_surf = inputs.mtheta
 
     fill!(vac_data.wv, 0)

test/runtests_error_fields.jl

@@ -370,6 +370,10 @@ include("h5_metadata_check.jl")
             @test norm(two.resonant_field_least_squares_t) < 1e-8 * norm(two.resonant_field_source_t)
             @test two.least_squares_rank == 2 && one.least_squares_rank == 1
             @test two.current_factor_least_squares ≈ [-1.0, 0.0] atol = 1e-6
+            # Identical arrays leave R rank-deficient; the minimum-norm solution splits the current evenly.
+            dup = EF.correction_requirement(ctx_plot, src, [src, src])
+            @test dup.least_squares_rank == 1
+            @test dup.current_factor_least_squares ≈ [-0.5, -0.5] atol = 1e-6
             # With two arrays the dominant-mode condition is one equation in two unknowns: the
             # minimum-current member of its family cancels the mode exactly and carries less current
             # than the exact solution (−1, 0), which is also in the family.

test/runtests_ntv_limits.jl

@@ -35,6 +35,10 @@ using HDF5
         @test EF.correction_current(3δt, c; delta_threshold=δt, torque_budget=T0, safety_factor=2.0, ntv=false) ≈ δt / c.delta_per_kat
         c0 = EF.EFCCoupling("x", c.delta_per_kat, 40.0, 0.05, 0.0)
         @test EF.correction_current(3δt, c0; delta_threshold=δt, torque_budget=T0) ≈ 2δt / c.delta_per_kat
+        # A residual torque too weak to resolve in C² − 4a·excess still gives the linear current, not zero.
+        cw = EF.EFCCoupling("w", c.delta_per_kat, 40.0, 0.05, 1.0e-20)
+        @test EF.correction_current(3δt, cw; delta_threshold=δt, torque_budget=T0) ≈ 2δt / c.delta_per_kat rtol = 1e-10
+        @test EF.correction_current_upper(3δt, cw; delta_threshold=δt, torque_budget=T0) ≈ 4δt / c.delta_per_kat rtol = 1e-10
         # Beyond the correctable limit there is no valid root.
         lim = EF.max_correctable_overlap(c; delta_threshold=δt, torque_budget=T0)
         @test isnan(EF.correction_current(1.01 * lim.with_ntv, c; delta_threshold=δt, torque_budget=T0))

test/runtests_risk.jl

@@ -137,6 +137,12 @@ using Statistics
         @test isapprox(s25, 2 / 3; rtol=0.05)
         @test isnan(EF.allowable_tolerance(scan, 99.0))     # never reached within the scan
         @test_throws ArgumentError EF.allowable_tolerance(scan, 0.0)
+        # A scan from scale 0 (the as-designed risk) interpolates that interval linearly and returns sampled points exactly.
+        z = fill(0.0, 3)
+        scan0 = EF.ToleranceScan([0.0, 1.0, 2.0], [0.5, 2.0, 10.0], [0.5, 2.0, 10.0], z, z, 0.5)
+        @test EF.allowable_tolerance(scan0, 0.5) == 0.0
+        @test EF.allowable_tolerance(scan0, 2.0) == 1.0
+        @test EF.allowable_tolerance(scan0, 1.0) ≈ 0.5
         @test_throws ArgumentError EF.tolerance_scan(table, ts, sets, mc_ctrl, sharp, scen; scales=Float64[])
     end
 

test/runtests_sing.jl

@@ -172,6 +172,25 @@ using FastInterpolations: cubic_interp, CubicFit, LinearBinarySearch, Series, Ex
         end
     end
 
+    @testset "q extrema of a quadratic profile and flat segments" begin
+        FFS = GeneralizedPerturbedEquilibrium.ForceFreeStates
+        equil = load_equilibrium_from_gpec(joinpath(@__DIR__, "test_data", "regression_solovev_ideal_example", "gpec.toml"))
+        # On a uniform grid a quadratic q is fit exactly, so most cells have q‴ = 0; put the minimum on a knot.
+        xs = collect(range(0.0, 1.0; length=129))
+        ψmin = xs[52]
+        q_quad = @. 1.8 + 3 * (xs - ψmin)^2
+        equil.profiles = GeneralizedPerturbedEquilibrium.Equilibrium.ProfileSplines(xs, ones(129), zeros(129), ones(129), q_quad)
+        GeneralizedPerturbedEquilibrium.Equilibrium.equilibrium_qfind!(equil)
+        @test equil.params.qmin ≈ 1.8 atol = 1e-10
+        @test any(ψ -> isapprox(ψ, ψmin; atol=1e-10), equil.params.qextrema_psi)
+        @test count(s -> s.m == 2, FFS._find_rational_surfaces(equil, 1, 1)) == 2
+        # Two extrema at the same rational q bound no isolated crossing, and the search must still end.
+        equil.params.qextrema_psi = [xs[1], xs[end]]
+        equil.params.qextrema_q = [2.0, 2.0]
+        equil.params.mextrema = 2
+        @test isempty(FFS._find_rational_surfaces(equil, 1, 1))
+    end
+
     @testset "sing_find reverse shear" begin
         FFS = GeneralizedPerturbedEquilibrium.ForceFreeStates
         equil = load_equilibrium_from_gpec(joinpath(@__DIR__, "test_data", "regression_solovev_ideal_example", "gpec.toml"))

test/runtests_vacuum.jl

@@ -61,6 +61,12 @@ using GeneralizedPerturbedEquilibrium.Vacuum: ConformalWall, EllipticalWall, Dee
             @test ModDeeWall().center_offset == 0.0
         end
 
+        @testset "2D response rejects non-contiguous modes" begin
+            inputs = VacuumInput(; mtheta_in=5, nzeta_in=1, x=[1.0, 1.1, 1.2, 1.1, 1.0], z=[0.0, -0.1, 0.0, 0.1, 0.0], ν=zeros(5),
+                mtheta=5, m_modes=[1, 3], n_modes=[1], nzeta=1)
+            @test_throws ArgumentError GeneralizedPerturbedEquilibrium.Vacuum.compute_vacuum_response(inputs, nothing)
+        end
+
         @testset "PlasmaGeometry" begin
             @testset "from VacuumInput" begin
                 inputs = VacuumInput(;