From 77807c61707f29e2ac63b284bc2a2509c03e17f5 Mon Sep 17 00:00:00 2001 From: logan-nc Date: Mon, 9 Nov 2026 09:09:33 -0500 Subject: [PATCH] Utilities - BUGFIX - Use the Lin-Liu & Miller upper bound in trapped_fraction trapped_fraction used h = B_min/B_max and the inverted ratio ^2/, which overestimates the effective trapped fraction by 25-40 % (0.7745 vs exact 0.6129 at epsilon = 0.2). It now returns the Lin-Liu & Miller 1995 upper bound f_tu with h = /B_max and /^2, within 2 % of the exact value (above it, as an upper bound) for epsilon between 0.05 and 0.4. B_min is kept in the signature but is no longer used. Co-Authored-By: Claude Opus 5.5 Claude-Session: https://claude.ai/code/session_01BgNxCBwLNHQmVq7iXMCTae --- src/Utilities/NeoclassicalResistivity.jl | 14 ++++++++------ test/runtests_utilities.jl | 19 +++++++++++++++++++ 2 files changed, 27 insertions(+), 6 deletions(-) diff --git a/src/Utilities/NeoclassicalResistivity.jl b/src/Utilities/NeoclassicalResistivity.jl index a4b15dde18..3240033da0 100644 --- a/src/Utilities/NeoclassicalResistivity.jl +++ b/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 diff --git a/test/runtests_utilities.jl b/test/runtests_utilities.jl index b4a5a92717..ac7ce86dca 100644 --- a/test/runtests_utilities.jl +++ b/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