Utilities - BUGFIX! - Use the Lin-Liu & Miller upper bound in trapped_fraction

src/Utilities/NeoclassicalResistivity.jl

@@ -4,7 +4,7 @@
 #   - eq 18b (nu*_e):         Sauter et al. 1999, Eq. (18b)
 #   - eq 13 (F_33 Sauter):    Sauter et al. 1999, Eqs. (13a)-(13b)
 #   - eq 17 (F_33 Redl):      Redl et al. 2021, Eqs. (17)-(18)
-#   - f_t (Lin-Liu & Miller): Phys. Plasmas 2, 1666 (1995), Eq. (6)
+#   - f_t (Lin-Liu & Miller): Phys. Plasmas 2, 1666 (1995), Eq. (4)
 #   - NRL Coulomb log:        NRL Plasma Formulary 2009
 
 """
@@ -178,21 +178,23 @@ end
 """
     trapped_fraction(avg_B, avg_Bsq, B_min, B_max) -> Float64
 
-Lin-Liu & Miller 1995, Phys. Plasmas **2**, 1666, Eq. (6):
+Lin-Liu & Miller 1995, Phys. Plasmas **2**, 1666, Eq. (4), upper bound on the effective trapped fraction
+(about 3 % above the exact value as ε → 0, closer at larger ε):
 
 ```
-f_t = 1 − ⟨B⟩² / ⟨B²⟩ · (1 − √(1 − h) · (1 + h/2)),   h = B_min / B_max
+f_tu = 1 − ⟨B²⟩ / ⟨B⟩² · (1 − √(1 − h) · (1 + h/2)),   h = ⟨B⟩ / B_max
 ```
 
-Arguments are flux-surface averages and extrema of |B| on the surface.
+Arguments are flux-surface averages and the maximum of |B| on the surface.
+`B_min` is accepted for call-site compatibility and is unused.
 """
 function trapped_fraction(avg_B::Real, avg_Bsq::Real,
     B_min::Real, B_max::Real)
     B_max > 0 || throw(ArgumentError("trapped_fraction: B_max must be > 0"))
     avg_Bsq > 0 || throw(ArgumentError("trapped_fraction: avg_Bsq must be > 0"))
-    h = clamp(B_min / B_max, 0.0, 1.0)
+    h = clamp(avg_B / B_max, 0.0, 1.0)
     factor = 1.0 - sqrt(1.0 - h) * (1.0 + 0.5 * h)
-    ft = 1.0 - (avg_B^2 / avg_Bsq) * factor
+    ft = 1.0 - (avg_Bsq / avg_B^2) * factor
     return clamp(ft, 0.0, 1.0)
 end
 

test/runtests_utilities.jl

@@ -1,3 +1,5 @@
+import QuadGK
+
 @testset "Utilities Unit Tests" begin
     @testset "adaptive_sample resolves a narrow resonance and brackets a zero crossing" begin
         adaptive_sample = GeneralizedPerturbedEquilibrium.Utilities.adaptive_sample
@@ -72,4 +74,21 @@
         GeneralizedPerturbedEquilibrium.Utilities.get_complex_coeffs!(out, fc, 10, 1)
         @test out[3] == c2  # Mode 2 is at index 3 (0-indexed mode)
     end
+
+    @testset "trapped_fraction is a tight upper bound on the exact effective trapped fraction" begin
+        trapped_fraction = GeneralizedPerturbedEquilibrium.Utilities.NeoclassicalResistivity.trapped_fraction
+        for eps in (0.001, 0.05, 0.2, 0.4)
+            theta = range(0, 2pi; length=2001)[1:(end-1)]
+            B = 1 ./ (1 .+ eps .* cos.(theta))
+            w = 1 ./ B
+            avg(g) = sum(g .* w) / sum(w)
+            B_max = maximum(B)
+            # Exact: f_t = 1 − (3/4)⟨B²⟩ ∫₀^{1/B_max} λ dλ / ⟨√(1 − λB)⟩
+            integral, _ = QuadGK.quadgk(lam -> lam / avg(sqrt.(max.(1 .- lam .* B, 0.0))), 0.0, 1 / B_max; rtol=1e-8)
+            f_exact = 1 - 0.75 * avg(B .^ 2) * integral
+            f_tu = trapped_fraction(avg(B), avg(B .^ 2), minimum(B), B_max)
+            @test f_tu >= f_exact
+            @test f_tu <= 1.03 * f_exact  # 1.5/1.46 as eps -> 0 (Lin-Liu & Miller Eqs. 14-16)
+        end
+    end
 end