|
| 1 | +# Quantifies whether splining iota = 1/q instead of q would reduce the knot count needed |
| 2 | +# for a given edge accuracy (GitHub discussion around the auto-grid redesign; PR #179 |
| 3 | +# lineage owns any production iota work). For a grid ending at psihigh < 1, interpolation |
| 4 | +# error transforms as err_q ≈ err_iota · q², cancelling the flattening of the profile, so |
| 5 | +# no leading-order gain is expected — this script documents that with numbers on the |
| 6 | +# DIII-D-like example equilibrium. |
| 7 | +# |
| 8 | +# Usage: julia --project=. benchmarks/benchmark_q_vs_iota_edge.jl |
| 9 | +# Outputs (not committed): benchmarks/q_vs_iota_edge.csv, benchmarks/q_vs_iota_edge.png |
| 10 | + |
| 11 | +using Pkg; |
| 12 | +Pkg.activate(joinpath(@__DIR__, "..")) |
| 13 | +using GeneralizedPerturbedEquilibrium |
| 14 | +using FastInterpolations |
| 15 | +using Printf |
| 16 | +using Plots |
| 17 | + |
| 18 | +const GPE = GeneralizedPerturbedEquilibrium |
| 19 | +const EXAMPLE_DIR = joinpath(@__DIR__, "..", "examples", "DIIID-like_ideal_example") |
| 20 | + |
| 21 | +# Dense ldp reference equilibrium: treat its q(ψ) as ground truth |
| 22 | +function reference_q() |
| 23 | + _, eq_config, additional_input = GPE.build_inputs_from_toml(EXAMPLE_DIR) |
| 24 | + eq_config.grid_type = "ldp" |
| 25 | + eq_config.mpsi = 1024 |
| 26 | + equil = GPE.Equilibrium.setup_equilibrium(eq_config, additional_input) |
| 27 | + return equil, eq_config |
| 28 | +end |
| 29 | + |
| 30 | +# Rational surface ψ_s(q = m/n) and q'(ψ_s) recovered from a fitted spline by bisection |
| 31 | +function rational_locations(q_itp, dq_itp, psis, q_targets) |
| 32 | + out = Dict{Float64,Tuple{Float64,Float64}}() |
| 33 | + for qt in q_targets |
| 34 | + lo, hi = psis[1], psis[end] |
| 35 | + (q_itp(lo) - qt) * (q_itp(hi) - qt) > 0 && continue |
| 36 | + for _ in 1:200 |
| 37 | + mid = 0.5 * (lo + hi) |
| 38 | + (q_itp(lo) - qt) * (q_itp(mid) - qt) <= 0 ? (hi = mid) : (lo = mid) |
| 39 | + end |
| 40 | + ps = 0.5 * (lo + hi) |
| 41 | + out[qt] = (ps, dq_itp(ps)) |
| 42 | + end |
| 43 | + return out |
| 44 | +end |
| 45 | + |
| 46 | +function main() |
| 47 | + equil, eq_config = reference_q() |
| 48 | + xs_ref = equil.profiles.xs |
| 49 | + q_ref_itp = equil.profiles.q_spline |
| 50 | + dq_ref_itp = equil.profiles.q_deriv |
| 51 | + psilow, psihigh = xs_ref[1], xs_ref[end] |
| 52 | + |
| 53 | + q_lo = ceil(q_ref_itp(psilow)) |
| 54 | + q_hi = floor(q_ref_itp(psihigh)) |
| 55 | + q_targets = collect(q_lo:q_hi) # n=1 rationals |
| 56 | + ref_rats = rational_locations(q_ref_itp, dq_ref_itp, xs_ref, q_targets) |
| 57 | + |
| 58 | + psi_eval = collect(range(psilow, psihigh; length=4001)) |
| 59 | + edge_band = psi_eval .> 0.9 |
| 60 | + |
| 61 | + rows = String[] |
| 62 | + push!(rows, "N,repr,max_rel_q_err,edge_rel_q_err,edge_rel_qprime_err,max_psi_s_err,max_rel_q1_err") |
| 63 | + results = Dict{String,Vector{NTuple{2,Float64}}}("q" => [], "iota" => []) |
| 64 | + |
| 65 | + for N in (16, 32, 64, 128, 256) |
| 66 | + # Same log_asymptotic knot layout the auto grid uses for fixed mpsi |
| 67 | + cfg = deepcopy(eq_config) |
| 68 | + cfg.mpsi = N |
| 69 | + knots = GPE.Equilibrium._build_psi_grid(cfg, psilow, psihigh) |
| 70 | + q_nodes = [q_ref_itp(p) for p in knots] |
| 71 | + |
| 72 | + for repr in ("q", "iota") |
| 73 | + itp = repr == "q" ? cubic_interp(knots, q_nodes; extrap=ExtendExtrap()) : |
| 74 | + cubic_interp(knots, 1.0 ./ q_nodes; extrap=ExtendExtrap()) |
| 75 | + ditp = deriv1(itp) |
| 76 | + qf = repr == "q" ? (p -> itp(p)) : (p -> 1.0 / itp(p)) |
| 77 | + dqf = repr == "q" ? (p -> ditp(p)) : (p -> -ditp(p) / itp(p)^2) # q' = -ι'/ι² |
| 78 | + |
| 79 | + q_err = [abs(qf(p) - q_ref_itp(p)) / abs(q_ref_itp(p)) for p in psi_eval] |
| 80 | + dq_err = [abs(dqf(p) - dq_ref_itp(p)) / max(abs(dq_ref_itp(p)), 1e-10) for p in psi_eval] |
| 81 | + rats = rational_locations(qf, dqf, knots, q_targets) |
| 82 | + psi_s_err = maximum(abs(rats[qt][1] - ref_rats[qt][1]) for qt in keys(ref_rats); init=0.0) |
| 83 | + q1_err = maximum(abs(rats[qt][2] - ref_rats[qt][2]) / abs(ref_rats[qt][2]) for qt in keys(ref_rats); init=0.0) |
| 84 | + |
| 85 | + push!(rows, @sprintf("%d,%s,%.3e,%.3e,%.3e,%.3e,%.3e", |
| 86 | + N, repr, maximum(q_err), maximum(q_err[edge_band]), maximum(dq_err[edge_band]), psi_s_err, q1_err)) |
| 87 | + push!(results[repr], (Float64(N), maximum(q_err[edge_band]))) |
| 88 | + end |
| 89 | + end |
| 90 | + |
| 91 | + csv_path = joinpath(@__DIR__, "q_vs_iota_edge.csv") |
| 92 | + open(csv_path, "w") do io |
| 93 | + foreach(r -> println(io, r), rows) |
| 94 | + end |
| 95 | + println("Wrote ", abspath(csv_path)) |
| 96 | + foreach(println, rows) |
| 97 | + |
| 98 | + p = plot(; xscale=:log10, yscale=:log10, xlabel="knots N", ylabel="max relative q error, ψ > 0.9", |
| 99 | + title="q-spline vs ι-spline edge accuracy", legend=:topright, left_margin=12Plots.mm, bottom_margin=4Plots.mm) |
| 100 | + for (repr, pts) in results |
| 101 | + plot!(p, first.(pts), last.(pts); marker=:circle, lw=2, label="$repr-spline") |
| 102 | + end |
| 103 | + png_path = joinpath(@__DIR__, "q_vs_iota_edge.png") |
| 104 | + savefig(p, png_path) |
| 105 | + println("Wrote ", abspath(png_path)) |
| 106 | +end |
| 107 | + |
| 108 | +main() |
0 commit comments