diff --git a/benchmarks/benchmark_delta_crit_radial_label.jl b/benchmarks/benchmark_delta_crit_radial_label.jl new file mode 100644 index 000000000..7c5d87849 --- /dev/null +++ b/benchmarks/benchmark_delta_crit_radial_label.jl @@ -0,0 +1,168 @@ +# Radial-label sensitivity of the SLAYER critical-Δ threshold. +# +# The tearing threshold is Δ'_rs > Δ_crit, with Δ'_rs the outer Δ' converted to the r_s +# reference through K^(2μ), K = r_s·dψ_N/dr. The Connor et al. 2015 Eq. 59 (`:toroidal`) +# critical-Δ scales with the radial label through the same K, so its threshold margin should +# be nearly label-invariant; the cylindrical `:rfitzp` formula is not covariant under a +# relabeling. Two cases: +# - the DIII-D-like SLAYER example (shaped): margin Δ'_rs/Δ_crit per label and branch; +# - the TJ-analytic circular ε scan: on a circular equilibrium every label coincides to +# O(ε²), so the label spread and the toroidal/rfitzp ratio must both approach their +# large-aspect-ratio limits (0 and 1) as ε → 0; a residual is a bug, not a convention. +# Each case computes a fresh Riccati Δ' matrix (no HDF5 output, no PE, no SLAYER stage). +# +# Usage: julia --project=. benchmarks/benchmark_delta_crit_radial_label.jl [--no-diiid] [--no-tj] +# Outputs go to benchmarks/delta_crit_radial_label/ (not committed). +using Printf, TOML, Plots +using GeneralizedPerturbedEquilibrium +const GPE = GeneralizedPerturbedEquilibrium +using GeneralizedPerturbedEquilibrium.Equilibrium: read_kinetic_file +using GeneralizedPerturbedEquilibrium.InnerLayer: build_slayer_inputs +using GeneralizedPerturbedEquilibrium.Utilities: KineticProfiles +using FastInterpolations: cubic_interp + +const EXAMPLE_D3D = joinpath(@__DIR__, "..", "examples", "DIIID-like_SLAYER_example") +const EXAMPLE_LAR = joinpath(@__DIR__, "..", "examples", "LAR_epsilon_scan") +const OUT = joinpath(@__DIR__, "delta_crit_radial_label") +const LABELS = (:midplane, :flux, :volume) +mkpath(OUT) + +# Outer-region solve only, from a run directory holding gpec.toml. +function outer_solve(dir) + inputs, eq_config, additional_input = GPE.build_inputs_from_toml(dir) + delete!(inputs, "SLAYER") + inputs["ForceFreeStates"]["write_outputs_to_HDF5"] = false + inputs["ForceFreeStates"]["force_termination"] = true + ffs = GPE.main_from_inputs(inputs, eq_config, additional_input, dir, "benchmark").ffs + return ffs.equil, ffs.surfaces, ffs.delta_prime.matrix +end + +# Per-surface rows (K, μ, Δ'_rs, Δ_crit for both branches) for one label. +function label_rows(equil, sings, dp, profiles, lab; kw...) + p_rf = build_slayer_inputs(equil, sings, profiles; dc_type=:rfitzp, rs_method=lab, kw...) + p_tor = build_slayer_inputs(equil, sings, profiles; dc_type=:toroidal, rs_method=lab, kw...) + dp_rs = real.(GPE.Tearing.Runner.delta_prime_to_rs_reference(dp, p_tor)) + return [(; label=lab, m=a.m, n=a.n, mn="$(a.m)/$(a.n)", K=a.k_ref, mu=a.mu_mercier, + dp_rs=dp_rs[k, k], dc_rf=a.dc_tmp, dc_tor=b.dc_tmp) for (k, (a, b)) in enumerate(zip(p_rf, p_tor))] +end + +function print_rows(rows) + @printf("%-9s %5s %7s %7s %8s %9s %9s %10s %10s %9s %9s\n", + "label", "m/n", "K", "mu", "K^2mu", "dp_rs", "dc_rf", "dc_tor", "dp/dc_rf", "dp/dc_tor", "tor/rf") + for r in rows + @printf("%-9s %5s %7.4f %7.4f %8.4f %9.3f %9.3f %10.3f %10.4f %9.4f %9.4f\n", + r.label, r.mn, r.K, r.mu, r.K^(2r.mu), r.dp_rs, r.dc_rf, r.dc_tor, r.dp_rs / r.dc_rf, r.dp_rs / r.dc_tor, r.dc_tor / r.dc_rf) + end +end + +# Label spread (max/min − 1) of the margin Δ'_rs/Δ_crit per surface, for both branches. +function margin_spread(rows, mn) + rs = filter(r -> r.mn == mn, rows) + m_rf = [abs(r.dp_rs / r.dc_rf) for r in rs] + m_tor = [abs(r.dp_rs / r.dc_tor) for r in rs] + return maximum(m_rf) / minimum(m_rf) - 1, maximum(m_tor) / minimum(m_tor) - 1 +end + +# --------------------------------------------------------------------------- +# DIII-D-like SLAYER example (shaped) +# --------------------------------------------------------------------------- +if !("--no-diiid" in ARGS) + println("=== DIII-D-like SLAYER example ===") + profile_file = TOML.parsefile(joinpath(EXAMPLE_D3D, "gpec.toml"))["SLAYER"]["profile_file"] + equil, sings, dp = outer_solve(EXAMPLE_D3D) + kin = read_kinetic_file(joinpath(EXAMPLE_D3D, profile_file)) + npsi = length(kin.psi) + profiles = KineticProfiles(; psi=kin.psi, n_e=kin.n_e, T_e=kin.T_e, T_i=kin.T_i, + omega=(kin.omega_E === nothing ? zeros(npsi) : kin.omega_E), omega_e=zeros(npsi), omega_i=zeros(npsi)) + _chi(v) = (v !== nothing && any(!=(0.0), v)) ? (let itp = cubic_interp(kin.psi, v); ψ -> Float64(itp(ψ)) end) : 1.0 + kw = (; chi_perp=_chi(kin.chi_e), chi_tor=_chi(kin.chi_phi)) + println("Δ' (ψ_N reference) diagonal: ", [round(real(dp[i, i]); digits=3) for i in 1:length(sings)]) + rows = reduce(vcat, [label_rows(equil, sings, dp, profiles, lab; kw...) for lab in LABELS]) + rows = filter(r -> r.m in (2, 3, 4), rows) + print_rows(rows) + mns = unique([r.mn for r in rows]) + println("\nLabel spread of the threshold margin |Δ'_rs/Δ_crit| (max/min over labels − 1):") + @printf("%5s %10s %10s\n", "m/n", "rfitzp", "toroidal") + for mn in mns + s_rf, s_tor = margin_spread(rows, mn) + @printf("%5s %10.4f %10.4f\n", mn, s_rf, s_tor) + end + x = 1:length(mns) + p = plot(; layout=(1, 2), size=(1000, 400), left_margin=8Plots.mm, bottom_margin=5Plots.mm, titlefontsize=10) + for (i, (dc, ttl)) in enumerate(((:dc_rf, "rfitzp"), (:dc_tor, "toroidal (Eq. 59, r_s ref.)"))) + plot!(p[i]; xticks=(x, mns), xlabel="rational surface m/n", ylabel="Δ'_rs / Δ_crit", + title="DIII-D-like threshold margin, $ttl", legend=(i == 1 ? :topright : false)) + for (j, lab) in enumerate(LABELS) + vals = [let r = only(filter(r -> r.mn == mn && r.label == lab, rows)); r.dp_rs / getfield(r, dc) end for mn in mns] + scatter!(p[i], x .+ 0.12 * (j - 2), vals; ms=6, label=String(lab)) + end + hline!(p[i], [1.0]; color=:black, ls=:dash, label="") + end + for ext in ("png", "pdf") + f = joinpath(OUT, "diiid_threshold_margin_by_label.$ext") + savefig(p, f) + println("saved: ", abspath(f)) + end +end + +# --------------------------------------------------------------------------- +# TJ-analytic circular ε scan +# --------------------------------------------------------------------------- +if !("--no-tj" in ARGS) + println("\n=== TJ-analytic circular ε scan ===") + # Baseline from the LAR ε-scan example; only lar_r0 (ε), the grid, and the on-axis + # pressure are overridden. pc = 0.01 (the LAR_resistive_match_test value) gives a + # finite D_R so Δ_crit is not numerically tiny; the margin ratios are D_R-independent. + base = TOML.parsefile(joinpath(EXAMPLE_LAR, "gpec.toml")) + epsilons = [0.05, 0.1, 0.2, 0.3] + # Synthetic flat-density, parabolic-temperature kinetic profiles (no kinetic file ships + # with the analytic equilibrium); they set χ∥ through the W_d loop but cancel in the ratios. + psi_k = collect(0.0:0.1:1.0) + nk = length(psi_k) + profiles = KineticProfiles(; psi=psi_k, n_e=fill(3.0e19, nk), T_e=2000.0 .* (1 .- 0.8 .* psi_k), + T_i=2000.0 .* (1 .- 0.8 .* psi_k), omega=zeros(nk), omega_e=zeros(nk), omega_i=zeros(nk)) + scan = [] + for eps in epsilons + run_dir = mktempdir(; prefix="gpec_dcrit_label_tj_") + cfg = deepcopy(base) + cfg["TJ_ANALYTIC_INPUT"]["lar_r0"] = cfg["TJ_ANALYTIC_INPUT"]["lar_a"] / eps + cfg["TJ_ANALYTIC_INPUT"]["pc"] = 0.01 + cfg["Equilibrium"]["mpsi"] = 128 + cfg["Equilibrium"]["mtheta"] = 256 + open(joinpath(run_dir, "gpec.toml"), "w") do io + TOML.print(io, cfg) + end + equil, sings, dp = outer_solve(run_dir) + rows = reduce(vcat, [label_rows(equil, sings, dp, profiles, lab; chi_perp=1.0, chi_tor=1.0) for lab in LABELS]) + println("\nε = $eps Δ' (ψ_N reference) diagonal: ", [round(real(dp[i, i]); digits=3) for i in 1:length(sings)]) + print_rows(rows) + for mn in unique([r.mn for r in rows]) + s_rf, s_tor = margin_spread(rows, mn) + r_mid = only(filter(r -> r.mn == mn && r.label == :midplane, rows)) + r_flx = only(filter(r -> r.mn == mn && r.label == :flux, rows)) + push!(scan, (; eps, mn, s_rf, s_tor, tor_rf_mid=r_mid.dc_tor / r_mid.dc_rf, tor_rf_flux=r_flx.dc_tor / r_flx.dc_rf)) + end + end + println("\nε scan summary (spread = max/min over labels − 1 of |Δ'_rs/Δ_crit|):") + @printf("%6s %5s %12s %12s %14s %14s\n", "eps", "m/n", "spread rf", "spread tor", "tor/rf mid", "tor/rf flux") + for s in scan + @printf("%6.3f %5s %12.4f %12.4f %14.4f %14.4f\n", s.eps, s.mn, s.s_rf, s.s_tor, s.tor_rf_mid, s.tor_rf_flux) + end + mns = unique([s.mn for s in scan]) + p = plot(; layout=(1, 2), size=(1000, 400), left_margin=8Plots.mm, bottom_margin=5Plots.mm, titlefontsize=10) + plot!(p[1]; xlabel="ε = a/R₀", ylabel="label spread of |Δ'_rs/Δ_crit|", title="TJ circular: label spread of the margin", legend=:topleft) + plot!(p[2]; xlabel="ε = a/R₀", ylabel="Δ_crit,toroidal / Δ_crit,rfitzp", title="TJ circular: toroidal / rfitzp", legend=:topleft) + for mn in mns + ss = filter(s -> s.mn == mn, scan) + plot!(p[1], [s.eps for s in ss], [s.s_rf for s in ss]; lw=2, marker=:circle, label="rfitzp $mn") + plot!(p[1], [s.eps for s in ss], [s.s_tor for s in ss]; lw=2, marker=:diamond, ls=:dash, label="toroidal $mn") + plot!(p[2], [s.eps for s in ss], [s.tor_rf_mid for s in ss]; lw=2, marker=:circle, label="midplane $mn") + plot!(p[2], [s.eps for s in ss], [s.tor_rf_flux for s in ss]; lw=2, marker=:diamond, ls=:dash, label="flux $mn") + end + hline!(p[2], [1.0]; color=:black, ls=:dot, label="") + for ext in ("png", "pdf") + f = joinpath(OUT, "tj_epsilon_scan_by_label.$ext") + savefig(p, f) + println("saved: ", abspath(f)) + end +end diff --git a/benchmarks/benchmark_toroidal_delta_crit.jl b/benchmarks/benchmark_toroidal_delta_crit.jl new file mode 100644 index 000000000..2b5880aa7 --- /dev/null +++ b/benchmarks/benchmark_toroidal_delta_crit.jl @@ -0,0 +1,183 @@ +# Verification and figures for the r_s-referenced Connor et al. 2015 Eq. 59 toroidal +# critical-Δ geometric factor (`toroidal_dgeo`): +# 1. large-aspect-ratio convergence on TJ-analytic circular equilibria, dgeo/√(n s r_s/R₀) +# vs ψ_N for several ε, for the correct Λ = ψ_t'² ι'/2π and for the Fortran STRIDE form +# that carries one power of ψ_t' (dimensional, off by ψ_t'^(-1/2)); +# 2. scale invariance under B₀ → B₀/2 and (a, R₀) → 2(a, R₀) at fixed ε; +# 3. an independent recomputation of ⟨B²⟩, ⟨|∇ψ_N|²⟩ and V' from R(ψ,θ), Z(ψ,θ), F(ψ), psio +# alone (no DCON metric elements), compared with ResistGeometry; +# 4. Δ_crit per rational surface on the DIII-D-like SLAYER example, rfitzp vs toroidal. +# +# Usage: julia --project=. benchmarks/benchmark_toroidal_delta_crit.jl +# Outputs go to benchmarks/toroidal_delta_crit/ (not committed). +using Printf, TOML, Plots +using GeneralizedPerturbedEquilibrium +using GeneralizedPerturbedEquilibrium.Equilibrium: TJAnalyticConfig, EquilibriumConfig, setup_equilibrium, read_kinetic_file +using GeneralizedPerturbedEquilibrium.ForceFreeStates: resist_geometry, ForceFreeStatesInternal, ForceFreeStatesControl, + sing_lim!, sing_find!, resist_eval_all! +using GeneralizedPerturbedEquilibrium.InnerLayer: toroidal_dgeo, r_based_shear, surface_minor_radius, surface_da_dpsi, + build_slayer_inputs +using GeneralizedPerturbedEquilibrium.Utilities: KineticProfiles +using FastInterpolations +using FastInterpolations: cubic_interp, Series, PeriodicBC, integrate + +const EXAMPLE_D3D = joinpath(@__DIR__, "..", "examples", "DIIID-like_SLAYER_example") +const OUT = joinpath(@__DIR__, "toroidal_delta_crit") +mkpath(OUT) + +function tj_equil(; eps=0.2, a=1.0, B0=12.0, mpsi=128, mtheta=256) + tj = TJAnalyticConfig(lar_r0=a / eps, lar_a=a, qc=1.5, qa=3.6, pc=0.001, mu=2.0, B0=B0, ma=128, mtau=128) + eq = EquilibriumConfig(eq_type="tj_analytic", psilow=0.01, psihigh=0.995, mpsi=mpsi, mtheta=mtheta, etol=1e-7) + return setup_equilibrium(eq, tj) +end + +# (dgeo, Fortran-form dgeo, rfitzp factor √(n s r_s/R₀), ResistGeometry) at one surface. +function dgeo_at(pe, psi, n) + chi1 = 2π * pe.psio + q = pe.profiles.q_spline(psi) + q1 = pe.profiles.q_deriv(psi) + rg = resist_geometry(pe, psi, q1) + rs = surface_minor_radius(pe, psi) + da = surface_da_dpsi(pe, psi) + dg = toroidal_dgeo(; chi1=chi1, v1=rg.v1_local, q=q, q1=q1, n=n, + avg_bsq=rg.avg_bsq, avg_dpsisq=rg.avg_dpsisq, k_ref=rs / da) + dg_fortran = dg / sqrt(q * chi1 / rg.v1_local) + lar = sqrt(n * r_based_shear(rs, q, q1, da) * rs / pe.ro) + return dg, dg_fortran, lar, rg +end + +# 1. Large-aspect-ratio convergence +println("1. Large-aspect-ratio convergence (n = 1)") +psis = collect(0.05:0.025:0.95) +p1 = plot(; xlabel="ψ_N", ylabel="dgeo / √(n s r_s/R₀)", title="Connor Eq. 59, r_s reference (Λ = ψ_t'² ι'/2π)", + legend=:topleft, xlims=(0, 1), left_margin=8Plots.mm, bottom_margin=5Plots.mm, titlefontsize=10) +p2 = plot(; xlabel="ψ_N", ylabel="dgeo / √(n s r_s/R₀)", title="Fortran STRIDE form (Λ with ψ_t'¹)", + legend=:topleft, xlims=(0, 1), left_margin=8Plots.mm, bottom_margin=5Plots.mm, titlefontsize=10) +hline!(p1, [1.0]; color=:black, ls=:dash, label="rfitzp") +hline!(p2, [1.0]; color=:black, ls=:dash, label="rfitzp") +for eps in (0.05, 0.1, 0.2, 0.3) + pe = tj_equil(; eps=eps) + rc = Float64[] + rf = Float64[] + for psi in psis + dg, dgf, lar, _ = dgeo_at(pe, psi, 1) + push!(rc, dg / lar) + push!(rf, dgf / lar) + end + @printf(" ε = %.2f ratio(correct) min/max = %.4f / %.4f ratio(Fortran) min/max = %.3f / %.3f\n", + eps, minimum(rc), maximum(rc), minimum(rf), maximum(rf)) + plot!(p1, psis, rc; lw=2, label=@sprintf("ε = %.2f", eps)) + plot!(p2, psis, rf; lw=2, label=@sprintf("ε = %.2f", eps)) +end +for ext in ("png", "pdf") + f = joinpath(OUT, "lar_convergence.$ext") + savefig(plot(p1, p2; layout=(1, 2), size=(1100, 420)), f) + println("saved: ", abspath(f)) +end + +# 2. Scale invariance +println("\n2. Scale invariance at ε = 0.2 (dgeo must be identical; the Fortran form scales as ψ_t'^(-1/2))") +ref, halfB, twice = tj_equil(), tj_equil(; B0=6.0), tj_equil(; a=2.0) +@printf("%6s %12s %12s %12s | %12s %12s %12s\n", "psi", "dgeo(ref)", "dgeo(B0/2)", "dgeo(2a,2R)", "fort(ref)", "fort(B0/2)", "fort(2a,2R)") +for psi in (0.3, 0.6, 0.9) + d0, f0 = dgeo_at(ref, psi, 1) + d1, f1 = dgeo_at(halfB, psi, 1) + d2, f2 = dgeo_at(twice, psi, 1) + @printf("%6.2f %12.6f %12.6f %12.6f | %12.6f %12.6f %12.6f\n", psi, d0, d1, d2, f0, f1, f2) +end + +# 3. Independent metric: ψ_N is axisymmetric, so |∇ψ_N|² = |∂_θ(R,Z)|²/J₂² with J₂ = R_ψ Z_θ − R_θ Z_ψ; +# B² = (F/R)² + psio²|∇ψ_N|²/R² with F = F_spline/(2π); volume element 2πR|J₂| dψ dθ. +function independent_averages(pe, psi) + ys = pe.rzphi_ys + F = pe.profiles.F_spline(psi) / (2π) + w = zeros(length(ys)) + b2 = zeros(length(ys)) + g2 = zeros(length(ys)) + for (i, th) in enumerate(ys) + f1 = pe.rzphi_rsquared((psi, th)) + f2 = pe.rzphi_offset((psi, th)) + f1p = pe.rzphi_rsquared((psi, th); deriv=DerivOp(1, 0)) + f1t = pe.rzphi_rsquared((psi, th); deriv=DerivOp(0, 1)) + f2p = pe.rzphi_offset((psi, th); deriv=DerivOp(1, 0)) + f2t = pe.rzphi_offset((psi, th); deriv=DerivOp(0, 1)) + rfac = sqrt(f1) + eta = 2π * (th + f2) + R = pe.ro + rfac * cos(eta) + rp, rt = f1p / (2rfac), f1t / (2rfac) + ep, et = 2π * f2p, 2π * (1 + f2t) + Rp = rp * cos(eta) - rfac * sin(eta) * ep + Zp = rp * sin(eta) + rfac * cos(eta) * ep + Rt = rt * cos(eta) - rfac * sin(eta) * et + Zt = rt * sin(eta) + rfac * cos(eta) * et + J2 = Rp * Zt - Rt * Zp + gpsi2 = (Rt^2 + Zt^2) / J2^2 + w[i] = 2π * R * abs(J2) + g2[i] = gpsi2 + b2[i] = (F / R)^2 + pe.psio^2 * gpsi2 / R^2 + end + s = integrate(cubic_interp(ys, Series(hcat(w, w .* b2, w .* g2)); bc=PeriodicBC())) + return (; v1=s[1], avg_bsq=s[2] / s[1], avg_dpsisq=s[3] / s[1]) +end + +function metric_table(pe, psis, label) + println("\n3. Independent metric check: ", label) + @printf("%6s %14s %14s %10s | %14s %14s %10s | %12s %12s %10s\n", + "psi", " RG", " indep", "rel", "<|dpsi|2> RG", "indep", "rel", "v1 RG", "v1 indep", "rel") + for psi in psis + rg = resist_geometry(pe, psi, pe.profiles.q_deriv(psi)) + ia = independent_averages(pe, psi) + @printf("%6.4f %14.6e %14.6e %10.2e | %14.6e %14.6e %10.2e | %12.5e %12.5e %10.2e\n", + psi, rg.avg_bsq, ia.avg_bsq, abs(ia.avg_bsq / rg.avg_bsq - 1), + rg.avg_dpsisq, ia.avg_dpsisq, abs(ia.avg_dpsisq / rg.avg_dpsisq - 1), + rg.v1_local, ia.v1, abs(ia.v1 / rg.v1_local - 1)) + end +end +metric_table(ref, (0.3, 0.6, 0.9), "TJ circular ε = 0.2") + +# 4. DIII-D-like example: Δ_crit per rational surface +inputs = TOML.parsefile(joinpath(EXAMPLE_D3D, "gpec.toml")) +equil = setup_equilibrium(EquilibriumConfig(inputs["Equilibrium"], EXAMPLE_D3D), nothing) +ctrl = ForceFreeStatesControl(; (Symbol(k) => v for (k, v) in inputs["ForceFreeStates"])...) +intr = ForceFreeStatesInternal(; dir_path=EXAMPLE_D3D) +intr.nlow = ctrl.nn_low +intr.nhigh = ctrl.nn_high +intr.npert = 1 +sing_lim!(intr, ctrl, equil) +sing_find!(intr, equil) +resist_eval_all!(intr, equil) +sings = intr.sing +metric_table(equil, [s.psifac for s in sings], "DIII-D-like rational surfaces") + +kin = read_kinetic_file(joinpath(EXAMPLE_D3D, inputs["SLAYER"]["profile_file"])) +npsi = length(kin.psi) +profiles = KineticProfiles(; psi=kin.psi, n_e=kin.n_e, T_e=kin.T_e, T_i=kin.T_i, + omega=(kin.omega_E === nothing ? zeros(npsi) : kin.omega_E), omega_e=zeros(npsi), omega_i=zeros(npsi)) +_chi(v) = (v !== nothing && any(!=(0.0), v)) ? (let itp = cubic_interp(kin.psi, v); ψ -> Float64(itp(ψ)) end) : 1.0 +kw = (; chi_perp=_chi(kin.chi_e), chi_tor=_chi(kin.chi_phi)) +p_rf = build_slayer_inputs(equil, sings, profiles; dc_type=:rfitzp, kw...) +p_tor = build_slayer_inputs(equil, sings, profiles; dc_type=:toroidal, kw...) +println("\n4. DIII-D-like SLAYER example, per rational surface (midplane label)") +@printf("%5s %7s %7s %8s %8s %10s %9s %11s %11s %11s %8s\n", + "m/n", "psi_N", "q", "r_s", "s_r", "D_R", "dgeo", "√(nsr/R)", "dc_rfitzp", "dc_toroid", "ratio") +for (s, a, b) in zip(sings, p_rf, p_tor) + lar = sqrt(a.n * a.sval_r * a.rs / a.R0) + @printf("%5s %7.4f %7.3f %8.4f %8.4f %10.4f %9.4f %11.4f %11.4f %11.4f %8.4f\n", + "$(a.m)/$(a.n)", s.psifac, s.q, a.rs, a.sval_r, a.dr_val, b.dgeo_val, lar, a.dc_tmp, b.dc_tmp, b.dc_tmp / a.dc_tmp) +end +keep = [i for (i, a) in enumerate(p_rf) if a.m in (2, 3, 4) && a.n == 1] +labels = ["$(p_rf[i].m)/$(p_rf[i].n)" for i in keep] +x = 1:length(keep) +fig = plot(; xticks=(x, labels), xlabel="rational surface m/n", ylabel="Δ_crit (r_s reference)", + title="DIII-D-like, Δ_crit at the 2/1, 3/1, 4/1 surfaces (midplane label)", legend=:topleft, + left_margin=8Plots.mm, bottom_margin=5Plots.mm, titlefontsize=10, size=(560, 420), xlims=(0.5, length(keep) + 0.5)) +scatter!(fig, x .- 0.08, [p_rf[i].dc_tmp for i in keep]; ms=7, label="rfitzp") +scatter!(fig, x .+ 0.08, [p_tor[i].dc_tmp for i in keep]; ms=7, marker=:diamond, label="toroidal (Eq. 59, r_s ref.)") +for (xx, i) in zip(x, keep) + annotate!(fig, xx + 0.08, p_tor[i].dc_tmp, text(@sprintf(" ×%.2f", p_tor[i].dc_tmp / p_rf[i].dc_tmp), 8, :left)) +end +for ext in ("png", "pdf") + f = joinpath(OUT, "diiid_delta_crit_234.$ext") + savefig(fig, f) + println("saved: ", abspath(f)) +end diff --git a/src/ForceFreeStates/ResistEval.jl b/src/ForceFreeStates/ResistEval.jl index 04953b6c0..7719c4f5b 100644 --- a/src/ForceFreeStates/ResistEval.jl +++ b/src/ForceFreeStates/ResistEval.jl @@ -42,6 +42,7 @@ supporting flux-surface averages. | `M` | Mass factor | | `avg_bsq_over_dpsisq` | ⟨B²/|∇ψ|²⟩ — needed for τ_R | | `avg_bsq` | ⟨B²⟩ — needed for τ_R | +| `avg_dpsisq` | ⟨|∇ψ|²⟩ — needed for the toroidal critical-Δ factor | | `avg_B` | ⟨B⟩ — needed for Lin-Liu-Miller f_t | | `B_max`, `B_min` | θ-extrema of B on the surface [T] | | `f_trap` | Lin-Liu & Miller 1995 trapped-particle fraction | @@ -69,6 +70,7 @@ struct ResistGeometry M::Float64 avg_bsq_over_dpsisq::Float64 avg_bsq::Float64 + avg_dpsisq::Float64 avg_B::Float64 B_max::Float64 B_min::Float64 @@ -121,11 +123,11 @@ function resist_geometry(equil::Equilibrium.PlasmaEquilibrium, v2 = profiles.dVdpsi_deriv(psi_f) q = profiles.q_spline(psi_f) - # Build the 6 GGJ θ-integrands plus a 7th (B) for the neoclassical - # resistivity f_t calculation, and accumulate running extrema of + # Build the 6 GGJ θ-integrands plus B (neoclassical resistivity f_t) and + # |∇ψ|² (toroidal critical-Δ factor), and accumulate running extrema of # (B, R) for Lin-Liu-Miller f_t and the local ε. ntheta = length(equil.rzphi_ys) - ff = zeros(Float64, ntheta, 7) + ff = zeros(Float64, ntheta, 8) B_max = -Inf B_min = Inf R_max = -Inf @@ -163,6 +165,7 @@ function resist_geometry(equil::Equilibrium.PlasmaEquilibrium, ff[itheta, 5] = bsq ff[itheta, 6] = dpsisq / bsq ff[itheta, 7] = B_here + ff[itheta, 8] = dpsisq @views ff[itheta, :] .*= jac / v1 end @@ -190,7 +193,7 @@ function resist_geometry(equil::Equilibrium.PlasmaEquilibrium, return ResistGeometry( E_coef, F_coef, G_coef, H_coef, K_coef, M_coef, - avg[1], avg[5], + avg[1], avg[5], avg[8], avg_B, B_max, B_min, f_trap, R_major, eps_local, p, p1, v1, ) diff --git a/src/InnerLayer/InnerLayer.jl b/src/InnerLayer/InnerLayer.jl index 43d1fc214..78d2df4c4 100644 --- a/src/InnerLayer/InnerLayer.jl +++ b/src/InnerLayer/InnerLayer.jl @@ -25,7 +25,7 @@ import .GGJ: delta_convergence, solution_profile, asymptotic_profile, q4_surface import .SLAYER: SLAYERModel, SLAYERParameters, slayer_parameters, r_based_shear import .SLAYER: riccati_del_s, slayer_layer_thickness, LayerWidths -import .SLAYER: surface_minor_radius, surface_da_dpsi, build_slayer_inputs +import .SLAYER: surface_minor_radius, surface_da_dpsi, radial_label, build_slayer_inputs, toroidal_dgeo export InnerLayerModel, InnerLayerParameters, InnerLayerResponse, solve_inner, solve_inner_profile export GGJ, GGJModel, GGJParameters @@ -37,6 +37,6 @@ export delta_convergence, solution_profile, asymptotic_profile, q4_surface_bench export SLAYER, SLAYERModel, SLAYERParameters, slayer_parameters, r_based_shear export riccati_del_s, slayer_layer_thickness, LayerWidths -export surface_minor_radius, surface_da_dpsi, build_slayer_inputs +export surface_minor_radius, surface_da_dpsi, radial_label, build_slayer_inputs, toroidal_dgeo end # module InnerLayer diff --git a/src/InnerLayer/SLAYER/LayerInputs.jl b/src/InnerLayer/SLAYER/LayerInputs.jl index 96177903d..dba3c85a5 100644 --- a/src/InnerLayer/SLAYER/LayerInputs.jl +++ b/src/InnerLayer/SLAYER/LayerInputs.jl @@ -9,14 +9,14 @@ # Geometry extraction: # - Minor radius at the outboard midplane (θ = 0) via # `equil.rzphi_rsquared((ψ, 0.0))`. -# - `da/dψ` via central finite difference on the same bicubic. +# - `da/dψ` from the interpolant's own analytic ψ-derivative. # - r-based magnetic shear via `r_based_shear(rs, q, q1, da/dψ)` (defined # in LayerParameters.jl). using ..Utilities: KineticProfiles using ...Utilities.NeoclassicalResistivity: NeoResistivityModel, SpitzerModel, coulomb_log_e, nu_star_e -using FastInterpolations: DerivOp +using FastInterpolations: DerivOp, integrate, cubic_interp, cumulative_integrate, ExtendExtrap """ surface_minor_radius(equil, psi; theta=0.0) -> Float64 @@ -32,33 +32,157 @@ function surface_minor_radius(equil, psi::Real; theta::Real=0.0) end """ - surface_da_dpsi(equil, psi; theta=0.0, h=1e-5) -> Float64 + surface_da_dpsi(equil, psi; theta=0.0) -> Float64 -Central finite-difference approximation of `d(minor radius)/dψ` at `psi`. -Falls back to one-sided differences near the flux-coordinate boundaries -(0 or 1). +Analytic ψ-derivative of the minor radius at `psi` and poloidal angle +`theta`, taken from the `rzphi_rsquared` interpolant's own ψ-derivative as +`da/dψ = (∂r²/∂ψ) / (2a)`. Valid wherever the interpolant is, including +under extrapolation past the ψ grid. Diverges at the magnetic axis, where +`a ~ √ψ`; callers evaluating near `ψ = 0` must check `isfinite`. """ -function surface_da_dpsi(equil, psi::Real; theta::Real=0.0, h::Real=1e-5) - psi_f = Float64(psi) - # Clamp to safe sampling range within (0, 1) - eps_edge = 10 * h - lo = psi_f - h - hi = psi_f + h - if lo < eps_edge - # one-sided forward - a0 = surface_minor_radius(equil, max(psi_f, eps_edge); theta=theta) - a1 = surface_minor_radius(equil, max(psi_f, eps_edge) + h; theta=theta) - return (a1 - a0) / h - elseif hi > 1.0 - eps_edge - # one-sided backward - a0 = surface_minor_radius(equil, min(psi_f, 1.0 - eps_edge) - h; theta=theta) - a1 = surface_minor_radius(equil, min(psi_f, 1.0 - eps_edge); theta=theta) - return (a1 - a0) / h +function surface_da_dpsi(equil, psi::Real; theta::Real=0.0) + return _da_dpsi_at_theta(equil, Float64(psi), Float64(theta)) +end + +# d(√r²)/dψ at one (ψ, θ) from the interpolant's own ψ-derivative. Shared by all radial-label +# conventions so none carries its own stencil. +@inline function _da_dpsi_at_theta(equil, psi::Float64, theta::Float64) + r_sq = equil.rzphi_rsquared((psi, theta)) + a = sqrt(max(r_sq, 0.0)) + a > 0 || return Inf # magnetic axis: a ~ √ψ, so da/dψ genuinely diverges + return equil.rzphi_rsquared((psi, theta); deriv=DerivOp(1, 0)) / (2a) +end + +""" + radial_label(equil; rs_method=:midplane, theta=0.0) -> (r_at, dr_dpsi_at) + +Build the pair of closures `r_at(ψ)` and `dr_dpsi_at(ψ)` defining one radial +label for the layer stack. Both closures must come from the same label, +because the r-based shear `(r/q)dq/dr`, `τ_R = μ₀r²/η`, `τ_E = r²/χ`, the +`d_β/r` normalization, the Δ' reference-length factor `k_ref = r_s/(da/dψ)`, +and any metre-to-ψ width conversion all have to live in one coordinate; a +mismatched `r` and `dr/dψ` silently corrupts every one of them. All +derivatives are analytic — no label carries a finite-difference stencil. + +# Labels + + - `:midplane` -- outboard-midplane chord from the magnetic axis at `theta` + (historical default), natural for comparison with midplane diagnostics. + - `:halfwidth` -- midplane half-chord, the mean of the outboard and inboard + chords at `θ = 0` and `θ = 0.5`; shift-free. Coincides with the flux label + on circular equilibria but is its own convention on shaped ones. + - `:fsa` -- θ-mean surface radius, a 128-point midpoint mean of the local + minor radius; the closest geometric approximation to `:flux` at interior + surfaces of shaped equilibria. + - `:volume` -- cylinder-equivalent label `√(V(ψ)/(2π²R₀))`, the + Rutherford-literature convention. + - `:flux` -- toroidal-flux label. Fitzpatrick, Nucl. Fusion (2025), + Eq. 30: `dψ_p/dr = B₀ r g/q` integrates to `ψ_t = B₀r²/2`, so + `r = √(2ψ_t/B₀)` with `ψ_t = psio·∫₀^ψ (q/g) dψ′` and `g = F/(B₀R₀)`. + Defined from flux alone, it carries no circular-cross-section assumption, + and its derivative `dr/dψ ∝ q` grows toward a separatrix where the + geometric labels' `da/dψ` collapses. + +On shaped equilibria the labels agree at low-q surfaces and diverge strongly +near the edge, where the slab-layer matching is label-ambiguous regardless of +choice. The label is selected programmatically; it is not exposed via TOML. +""" +function radial_label(equil; rs_method::Symbol=:midplane, theta::Real=0.0) + theta_f = Float64(theta) + + _flux_r, _flux_dr = if rs_method === :flux + b0f = Float64(equil.params.b0) + R0f = Float64(equil.ro) + psiof = Float64(equil.psio) + xs_f = collect(Float64, equil.profiles.xs) + # g = F/(B0 R0) departs from 1 by ~3% on a DIII-D-like deck and ~0.8% on a circular + # one, so carry it rather than assuming g = 1: r² = 2∫(q/g)dψ_p/B0. + _g_at(x) = Float64(equil.profiles.F_spline(x)) / (2π * b0f * R0f) + qg = [Float64(equil.profiles.q_spline(x)) / _g_at(x) for x in xs_f] + Phi = collect(Float64, cumulative_integrate(cubic_interp(xs_f, qg))) + r_knots = sqrt.(max.(2 .* psiof .* Phi ./ b0f, 0.0)) + rspl = cubic_interp(xs_f, r_knots; extrap=ExtendExtrap()) + (ψ -> Float64(rspl(Float64(ψ))), + ψ -> psiof * Float64(equil.profiles.q_spline(Float64(ψ))) / + (b0f * _g_at(Float64(ψ)) * max(Float64(rspl(Float64(ψ))), eps()))) else - a_plus = surface_minor_radius(equil, psi_f + h; theta=theta) - a_minus = surface_minor_radius(equil, psi_f - h; theta=theta) - return (a_plus - a_minus) / (2h) + (nothing, nothing) end + + _a_at(ψ, θ) = sqrt(max(equil.rzphi_rsquared((Float64(ψ), Float64(θ))), 0.0)) + + _rs_at(ψ) = + if rs_method === :fsa + N = 128 + s = 0.0 + @inbounds for k in 1:N + s += _a_at(ψ, (k - 0.5) / N) + end + s / N + elseif rs_method === :halfwidth + 0.5 * (_a_at(ψ, 0.0) + _a_at(ψ, 0.5)) + elseif rs_method === :volume + V = integrate(equil.profiles.dVdpsi_spline, 1e-4, Float64(ψ)) + sqrt(max(V, 0.0) / (2π^2 * equil.ro)) + elseif rs_method === :flux + _flux_r(ψ) + else + surface_minor_radius(equil, ψ; theta=theta_f) + end + + _da_dpsi_at(ψ) = + if rs_method === :fsa + N = 128 + s = 0.0 + @inbounds for k in 1:N + s += _da_dpsi_at_theta(equil, Float64(ψ), (k - 0.5) / N) + end + s / N + elseif rs_method === :halfwidth + 0.5 * (_da_dpsi_at_theta(equil, Float64(ψ), 0.0) + _da_dpsi_at_theta(equil, Float64(ψ), 0.5)) + elseif rs_method === :volume + Float64(equil.profiles.dVdpsi_spline(ψ)) / (4π^2 * equil.ro * max(_rs_at(ψ), eps())) + elseif rs_method === :flux + _flux_dr(ψ) + else + _da_dpsi_at_theta(equil, Float64(ψ), theta_f) + end + + return (_rs_at, _da_dpsi_at) +end + +""" + toroidal_dgeo(; chi1, v1, q, q1, n, avg_bsq, avg_dpsisq, k_ref) -> Float64 + +Geometric factor of the toroidal critical-Δ, Connor, Ham, Hastie & Liu 2015 +(PPCF 57 065001) Eq. 59, `V_s·(α²Λ²/(⟨B²⟩⟨|∇V|²⟩))^{1/4}`, converted from the +paper's `Y = (V−V_s)/V_s` reference to the `x̂ = (r−r_s)/r_s` reference shared +by the slab layer, the `rfitzp` critical-Δ, and the reference-length-converted +outer Δ'. In GPEC quantities, with `V` the flux-surface volume, `ψ_N` the +normalized poloidal flux, and `'` = d/dV: + + - `α = 2πn/χ'` with `χ' = chi1/v1` (`χ` the full poloidal flux, + `chi1 = 2π·psio`, `v1 = dV/dψ_N`) + - `Λ = ψ_t'²·(ι/2π)'` with `ψ_t' = q·chi1/v1` (toroidal flux) and + `(ι/2π)' = −q1/(q²·v1)` (`q1 = dq/dψ_N`) + - `⟨|∇V|²⟩ = v1²·⟨|∇ψ_N|²⟩`, both averages being normalized flux-surface + averages (⟨1⟩ = 1) + - the reference conversion `r_s·(dV/dr)/V_s = k_ref·v1/V_s` with + `k_ref = r_s·dψ_N/dr`, so `V_s` cancels and is never integrated. + +At large aspect ratio this reduces to `√(n·s·r_s/R₀)` with `s = (r_s/q)·dq/dr`, +so `dc_type=:toroidal` coincides with `:rfitzp` there; the paper's own Eq. 61 +is recovered after dividing by `r_s`. +""" +function toroidal_dgeo(; chi1::Real, v1::Real, q::Real, q1::Real, n::Integer, + avg_bsq::Real, avg_dpsisq::Real, k_ref::Real) + v1 != 0 || throw(ArgumentError("toroidal_dgeo: dV/dψ must be non-zero")) + q != 0 || throw(ArgumentError("toroidal_dgeo: q must be non-zero")) + alpha = 2π * n * v1 / chi1 + psit1 = q * chi1 / v1 + lambda = psit1^2 * (-q1 / (q^2 * v1)) + grad_v_sq = v1^2 * avg_dpsisq + return k_ref * v1 * (alpha^2 * lambda^2 / (avg_bsq * grad_v_sq))^0.25 end """ @@ -106,14 +230,21 @@ profiles, without an intermediate file round-trip. NOT the Mercier index `D_I = E + F + H − 1/4`. The two differ by `(H − 1/2)²`, which is non-trivial on shaped equilibria (~factor 3 on DIII-D); this code uses the physically correct `D_R`. - - `dgeo_val` -- Connor 2015 (PPCF 57 065001) Eq. 59 geometric factor - used by `dc_type=:toroidal`. When `nothing` (default), an error is - raised if `dc_type=:toroidal` is also requested — the auto-derived - formula additionally needs ⟨|∇ψ|²⟩ FSA which `ResistGeometry` - doesn't currently expose. Pass a scalar / vector / callable to use - a prescribed value. (For `dc_type=:rfitzp` and `:lar`, dgeo_val is - not consulted.) + - `dgeo_val` -- Connor et al. 2015 (PPCF 57 065001) Eq. 59 geometric + factor of the toroidal critical-Δ, in the `r_s` reference (see + [`toroidal_dgeo`](@ref)). When `nothing` (default), it is derived + per-surface from the equilibrium through the surface's `ResistGeometry` + (`sing.restype`, populated by `ForceFreeStates.resist_eval_all!`); an + error is raised if `dc_type=:toroidal` is requested on a surface without + one. Pass a scalar / callable to use a prescribed value. Only + `dc_type=:toroidal` consumes it. - `dc_type` -- `:none` (default), `:lar`, `:rfitzp`, or `:toroidal`. + - `rs_method` -- radial label defining `r_s` for the whole layer stack: + `:midplane` (default), `:halfwidth`, `:fsa`, `:volume`, or `:flux`. See + [`radial_label`](@ref) for the definitions. S, the r-based shear, W_d, + and the Δ' reference-length factor `k_ref` all follow the choice + together, so every option is self-consistent. Not exposed via TOML — + programmatic use only. - `theta` -- poloidal angle at which to measure minor radius (default `0.0`, outboard midplane). - `resistivity_model` -- `SauterNeoModel()` (default), `RedlNeoModel()`, @@ -156,38 +287,7 @@ function build_slayer_inputs(equil, sings, profiles::KineticProfiles; Float64(bt(ψ)) end - # Minor-radius extractor: `:midplane` = outboard-midplane chord - # (original behavior); `:fsa` = θ-mean of √rzphi_rsquared, the - # flux-surface-averaged minor radius. - _rs_at(ψ) = - if rs_method === :fsa - integrand(θ) = sqrt(equil.rzphi_rsquared((Float64(ψ), Float64(θ)))) - N = 128 - s = 0.0 - @inbounds for k in 1:N - s += integrand((k - 0.5) / N) - end - s / N - else - surface_minor_radius(equil, ψ; theta=theta) - end - _da_dpsi_at(ψ) = - if rs_method === :fsa - # central finite difference on _rs_at - h = 1e-5 - lo = ψ - h - hi = ψ + h - eps_edge = 10h - if lo < eps_edge - (_rs_at(max(ψ, eps_edge) + h) - _rs_at(max(ψ, eps_edge))) / h - elseif hi > 1.0 - eps_edge - (_rs_at(min(ψ, 1.0 - eps_edge)) - _rs_at(min(ψ, 1.0 - eps_edge) - h)) / h - else - (_rs_at(ψ + h) - _rs_at(ψ - h)) / (2h) - end - else - surface_da_dpsi(equil, ψ; theta=theta) - end + _rs_at, _da_dpsi_at = radial_label(equil; rs_method=rs_method, theta=theta) # Per-surface ω_*e, ω_*i (diamagnetic frequencies) from spline # derivatives. When `compute_omega_star=true` we override any ω_*e/ω_*i @@ -268,27 +368,44 @@ function build_slayer_inputs(equil, sings, profiles::KineticProfiles; _eval(dr_val, psi) end - # dgeo_val: only used by dc_type=:toroidal (the Connor-Hastie- - # Helander 2015 formula). Auto-derivation requires ⟨|∇ψ|²⟩ FSA - # which the current `ResistGeometry` doesn't expose; for now we - # require an explicit value if the toroidal dc_type is selected. + # Reference-length conversion inputs for the outer Δ': K = r_s·(dψ_N/dr)|_s + # and μ = √(−D_I) with D_I = E + F + H − 1/4 (Glasser-Greene-Johnson 1975). + # A Mercier-unstable surface (D_I ≥ 0) has no real exponent; μ → 0 there, + # which makes the conversion factor K^(2μ) → 1 continuously (Δ' left raw). + k_ref_k = da_dpsi == 0.0 ? 1.0 : rs / da_dpsi + + # dgeo_val: Connor et al. 2015 Eq. 59 geometric factor in the r_s reference + # (see `toroidal_dgeo`), derived whenever the surface carries a ResistGeometry; + # only dc_type=:toroidal consumes it. dgeo_val_k = if dgeo_val === nothing - dc_type === :toroidal && + if rg !== nothing + toroidal_dgeo(; chi1=chi1, v1=rg.v1_local, q=q, q1=q1, n=n_res, + avg_bsq=rg.avg_bsq, avg_dpsisq=rg.avg_dpsisq, k_ref=k_ref_k) + elseif dc_type === :toroidal throw( ArgumentError( - "build_slayer_inputs: dc_type=:toroidal " * - "needs `dgeo_val` (Connor 2015 PPCF 57 " * - "065001 Eq. 59 geometric factor). " * - "Auto-derivation from equilibrium not " * - "yet implemented; pass a scalar / vector " * - "/ callable explicitly." + "build_slayer_inputs: dc_type=:toroidal with " * + "dgeo_val=nothing requires `sing.restype` populated " * + "by ForceFreeStates.resist_eval_all!. " * + "Surface k=$k has restype=nothing." ) ) - 0.0 + else + 0.0 + end else _eval(dgeo_val, psi) end + mu_k = if rg === nothing + @warn("build_slayer_inputs: sing.restype not populated; using the " * + "slab Mercier exponent μ = 1/2 for the Δ' reference-length " * + "conversion at all such surfaces.", maxlog=1) + 0.5 + else + sqrt(max(-(rg.E + rg.F + rg.H - 0.25), 0.0)) + end + out[k] = slayer_parameters(; n_e=prof.n_e, t_e=prof.T_e, t_i=prof.T_i, omega=prof.omega, omega_e=ω_e_use, omega_i=ω_i_use, @@ -304,7 +421,9 @@ function build_slayer_inputs(equil, sings, profiles::KineticProfiles; f_trap=f_trap_kw, nu_e_star=nu_e_star_kw, R_major_eff=R_major_eff, - lnLambda_form=lnLambda_form + lnLambda_form=lnLambda_form, + k_ref=k_ref_k, + mu_mercier=mu_k ) end return out diff --git a/src/InnerLayer/SLAYER/LayerParameters.jl b/src/InnerLayer/SLAYER/LayerParameters.jl index 5dbf9163f..0d32c01e5 100644 --- a/src/InnerLayer/SLAYER/LayerParameters.jl +++ b/src/InnerLayer/SLAYER/LayerParameters.jl @@ -41,11 +41,20 @@ de-normalization. The parametrization uses `P_perp`, `P_tor`, and | `bt` | Toroidal field [T] | | `sval_r` | r-based magnetic shear r_s · (dq/dr) / q (Fitzpatrick convention) | | `dr_val` | Resistive interchange D_R = E + F + H² (critical-Δ input; auto-derived from GGJ coefficients unless overridden) | -| `dgeo_val` | Connor-Hastie-Helander 2015 Eq. 59 geometric factor (0 unless supplied) | +| `dgeo_val` | Connor et al. 2015 Eq. 59 toroidal critical-Δ geometric factor in the r_s reference (see `toroidal_dgeo`) | | `eta` | Parallel resistivity entering τ_R = μ₀r_s²/η [Ω·m] | | `d_beta` | Beta-weighted ion length scale c_β · d_i [m] | | `dc_tmp` | Critical-Δ offset from chi_parallel matching | | `dc_type` | Selector for `dc_tmp` formula | +| `k_ref` | Reference-length ratio K = r_s · (dψ_N/dr) at this surface (1 = no Δ' conversion) | +| `mu_mercier` | Mercier exponent μ = √(−D_I) governing the Δ' reference-length conversion (1/2 = slab/cylindrical value) | + +`k_ref` and `mu_mercier` feed the ψ_N → r_s reference-length conversion of +the outer Δ' matrix (see `delta_prime_to_rs_reference` in the Tearing +runner): the slab layer, its `dc_tmp` critical-Δ, and the `S^(1/3)` Δ(Q) +scale are all referenced to unit `x̂ = (r−r_s)/r_s`, while the outer BVP Δ' +is referenced to unit Δψ_N. Hand-built parameters default to `k_ref = 1`, +which makes the conversion the identity. The complex normalized growth rate `Q = ω + iγ` is **not** stored here; it is passed as a separate argument to `solve_inner`. @@ -85,6 +94,10 @@ Base.@kwdef struct SLAYERParameters <: InnerLayerParameters # Critical-Δ offset dc_tmp::Float64 = 0.0 dc_type::Symbol = :none + + # Reference-length conversion inputs for the outer Δ' (ψ_N → r_s-based x̂) + k_ref::Float64 = 1.0 + mu_mercier::Float64 = 0.5 end # Allowed dc_type values for the critical-Δ offset. `:none` is the default @@ -117,7 +130,7 @@ function r_based_shear(rs::Real, q::Real, dq_dpsi::Real, da_dpsi::Real) end # Internal: solve the Wd self-consistency loop for the chi_parallel-based -# critical Δ (Connor-Hastie-Helander 2015). Returns dc_tmp as a Float64. +# critical Δ (Connor et al. 2015, PPCF 57 065001). Returns dc_tmp as a Float64. function _solve_dc_tmp(; dc_type::Symbol, dr_val::Real, dgeo_val::Real, chi_perp::Real, t_e::Real, zeff::Real, tau_ee::Real, rs::Real, R0::Real, sval_r::Real, n_tor::Integer, @@ -199,9 +212,17 @@ parametrization (P_perp/P_tor/D_norm; the older magnetic/electron Prandtl - `zeff` -- effective charge - `chi_perp`, `chi_tor` -- perpendicular / toroidal heat diffusivity [m²/s] - `m`, `n` -- poloidal / toroidal mode numbers at the surface - - `dr_val`, `dgeo_val` -- inputs for the critical-Δ formula + - `dr_val`, `dgeo_val` -- inputs for the critical-Δ formula: the resistive + interchange index `D_R` and the Connor et al. 2015 Eq. 59 geometric factor + in the `r_s` reference (`toroidal_dgeo`) - `dc_type` -- one of `:none`, `:lar`, `:rfitzp`, `:toroidal` - `ising` -- singular-surface index for traceability + - `k_ref` -- reference-length ratio K = r_s·(dψ_N/dr) at the surface, + used by the Tearing runner to convert the ψ_N-referenced outer Δ' to + the r_s-referenced convention this layer works in (default `1.0`, + i.e. no conversion; `build_slayer_inputs` fills it from the equilibrium) + - `mu_mercier` -- Mercier exponent μ = √(−D_I) for the same conversion + (default `0.5`, the slab/cylindrical value at D_I = −1/4) # Resistivity kwargs @@ -252,7 +273,9 @@ function slayer_parameters(; f_trap::Union{Real,Nothing}=nothing, nu_e_star::Union{Real,Nothing}=nothing, R_major_eff::Union{Real,Nothing}=nothing, - lnLambda_form::Symbol=:nrl) + lnLambda_form::Symbol=:nrl, + k_ref::Real=1.0, + mu_mercier::Real=0.5) # Coulomb logarithm shared by the resistivity closure and τ_ee. lnLamb = coulomb_log_e(n_e, t_e; form=lnLambda_form) @@ -358,6 +381,7 @@ function slayer_parameters(; rs=rs, R0=R0, bt=bt, sval_r=sval_r, dr_val=dr_val, dgeo_val=dgeo_val, eta=eta, d_beta=d_beta, - dc_tmp=dc_tmp, dc_type=dc_type + dc_tmp=dc_tmp, dc_type=dc_type, + k_ref=k_ref, mu_mercier=mu_mercier ) end diff --git a/src/InnerLayer/SLAYER/SLAYER.jl b/src/InnerLayer/SLAYER/SLAYER.jl index ece897c15..201e1c3c2 100644 --- a/src/InnerLayer/SLAYER/SLAYER.jl +++ b/src/InnerLayer/SLAYER/SLAYER.jl @@ -56,7 +56,7 @@ include("LayerInputs.jl") export SLAYERModel, SLAYERParameters, slayer_parameters export r_based_shear export riccati_del_s, slayer_layer_thickness, LayerWidths -export surface_minor_radius, surface_da_dpsi, build_slayer_inputs +export surface_minor_radius, surface_da_dpsi, radial_label, build_slayer_inputs, toroidal_dgeo export NeoResistivityModel, SpitzerModel, SpitzerHarmModel, SauterNeoModel, RedlNeoModel end # module SLAYER diff --git a/src/Tearing/Dispersion/SurfaceCoupling.jl b/src/Tearing/Dispersion/SurfaceCoupling.jl index 271162c1a..0d57207ae 100644 --- a/src/Tearing/Dispersion/SurfaceCoupling.jl +++ b/src/Tearing/Dispersion/SurfaceCoupling.jl @@ -53,10 +53,14 @@ end surface_coupling(model::SLAYERModel, params::SLAYERParameters, dp_diag::Number; dc::Real=0.0) -> SurfaceCoupling -SLAYER convenience constructor. `scale` is set to `params.lu^(1/3)` so that -the dimensionless Δ from `riccati_f` is mapped to outer ψ-units before -subtraction from the Δ' diagonal. `tauk` is taken from `params.tauk` for use -by `MultiSurfaceCoupling` Q rescaling. +SLAYER convenience constructor. `scale` is set to `params.lu^(1/3)`, which +maps the dimensionless inner-layer Δ from `riccati_f` to the r_s-referenced +outer convention (S^(1/3) with S = τ_R/τ_H built on r_s; the Lundquist +number, `dc`, and the layer all share the `x̂ = (r−r_s)/r_s` reference +length). `dp_diag` must already be in that same r_s reference — the Tearing +runner converts the ψ_N-referenced BVP Δ' via `delta_prime_to_rs_reference` +before building couplings. `tauk` is taken from `params.tauk` for use by +`MultiSurfaceCoupling` Q rescaling. """ function surface_coupling(model::SLAYERModel, params::SLAYERParameters, dp_diag::Number; dc::Real=0.0) diff --git a/src/Tearing/Runner/Control.jl b/src/Tearing/Runner/Control.jl index 1a12ff7e6..e95c69a62 100644 --- a/src/Tearing/Runner/Control.jl +++ b/src/Tearing/Runner/Control.jl @@ -41,7 +41,8 @@ constructor. - `dr_val`, `dgeo_val` -- critical-Δ formula inputs. `nothing` (default) auto-derives them from the equilibrium: `dr_val` from the resistive interchange index `D_R = E + F + H²` at each surface, `dgeo_val` from the - toroidal geometric factor (required only by `dc_type=:toroidal`). Supply a + Connor et al. 2015 Eq. 59 toroidal geometric factor in the `r_s` reference + (consumed only by `dc_type=:toroidal`). Supply a scalar only to override the auto-derivation; an explicit `0.0` disables the critical-Δ offset (Δ_crit ≡ 0) - `theta_sample` -- poloidal angle at which to sample minor radius diff --git a/src/Tearing/Runner/HDF5Output.jl b/src/Tearing/Runner/HDF5Output.jl index ca4287147..fe0a6fd4b 100644 --- a/src/Tearing/Runner/HDF5Output.jl +++ b/src/Tearing/Runner/HDF5Output.jl @@ -90,11 +90,14 @@ const TEARING_H5_ANNOTATIONS = [ "PerSurface/sval_r" => (; long_name="r-based magnetic shear r_s·(dq/dr)/q (Fitzpatrick convention)", dims=("surface",)), "PerSurface/D_R" => (; long_name="resistive interchange D_R = E + F + H² for the critical-Δ formula (auto-derived from GGJ coefficients unless overridden)", dims=("surface",)), - "PerSurface/D_geo" => (; long_name="Connor-Hastie-Helander 2015 Eq. 59 geometric factor (0 unless supplied)", dims=("surface",)), + "PerSurface/D_geo" => (; long_name="Connor et al. 2015 Eq. 59 toroidal critical-Δ geometric factor in the r_s reference (0 when no ResistGeometry)", dims=("surface",)), "PerSurface/eta" => (; long_name="parallel resistivity at each surface", units="Ohm*m", dims=("surface",)), "PerSurface/d_beta" => (; long_name="β-weighted ion drift scale d_β", units="m", dims=("surface",)), "PerSurface/D_c_offset" => (; long_name="critical-Δ offset from χ_∥/χ_⊥ matching (Connor-Hastie-Helander 2015 Eq. 59)", dims=("surface",)), "PerSurface/D_c_type" => (; long_name="per-surface D_c prescription label", dims=("surface",)), + "PerSurface/k_ref" => (; long_name="reference-length ratio K = r_s·(dψ_N/dr) at each surface", dims=("surface",)), + "PerSurface/mu_mercier" => (; long_name="Mercier exponent μ = √(−D_I) at each surface", dims=("surface",)), + "PerSurface/delta_prime_conversion" => (; long_name="ψ_N → r_s reference-length factor K^(2μ) applied to each Δ' diagonal", dims=("surface",)), "PerSurface/E" => (; long_name="Glasser-Greene-Johnson coefficient E per surface", dims=("surface",)), "PerSurface/F" => (; long_name="Glasser-Greene-Johnson coefficient F per surface", dims=("surface",)), "PerSurface/G" => (; long_name="Glasser-Greene-Johnson coefficient G per surface", dims=("surface",)), @@ -103,7 +106,7 @@ const TEARING_H5_ANNOTATIONS = [ "PerSurface/M" => (; long_name="Glasser-Greene-Johnson coefficient M per surface", dims=("surface",)), "PerSurface/tau_A" => (; long_name="Alfvén time τ_A per surface (GGJ layer parameters)", units="s", dims=("surface",)), "PerSurface/dVdpsi" => (; long_name="dV/dψ_N at each surface", units="m^3", dims=("surface",)), - "PerSurface/Delta_prime_matrix" => (; long_name="full complex Δ' matrix coupling the rational surfaces", dims=("surface_row", "surface_col")), + "PerSurface/Delta_prime_matrix" => (; long_name="full complex Δ' matrix coupling the rational surfaces, as used in the matching (SLAYER path: r_s-referenced via K^(2μ); the ψ_N-referenced BVP matrix is SingularSurfaces/Delta_prime_matrix)", dims=("surface_row", "surface_col")), "Roots/Q_root" => (; long_name="complex dispersion-root normalized frequency Q (NaN = no root)", dims=("surface",)), "Roots/omega" => (; long_name="mode rotation angular frequency ω = Re(Q)/τ_k of each root", units="rad/s", dims=("surface",)), @@ -184,6 +187,11 @@ function _write_per_surface!(g, params::AbstractVector{SLAYERParameters}, ps["D_geo"] = Float64[p.dgeo_val for p in params] # Store dc_type per-surface as string array ps["D_c_type"] = String[String(p.dc_type) for p in params] + # Reference-length conversion applied to Δ' (ψ_N → r_s): K, μ, and the + # diagonal factor K^(2μ) actually multiplying each Δ'_kk. + ps["k_ref"] = Float64[p.k_ref for p in params] + ps["mu_mercier"] = Float64[p.mu_mercier for p in params] + ps["delta_prime_conversion"] = Float64[p.k_ref^(2 * p.mu_mercier) for p in params] ps["Delta_prime_matrix"] = dp_matrix return nothing diff --git a/src/Tearing/Runner/Result.jl b/src/Tearing/Runner/Result.jl index 0cffc21c9..d81508989 100644 --- a/src/Tearing/Runner/Result.jl +++ b/src/Tearing/Runner/Result.jl @@ -20,7 +20,9 @@ downstream inspection and HDF5 output. factor q of each analyzed surface, aligned with `params`. Empty when the analysis was built from bare parameters (`run_slayer_from_inputs` without the surface list), in which case the HDF5 writer skips them. - - `dp_matrix` -- outer-region Δ' matrix used in the analysis + - `dp_matrix` -- outer-region Δ' matrix used in the analysis. + SLAYER path: r_s-referenced (the ψ_N BVP matrix transformed by + `delta_prime_to_rs_reference`); GGJ path: the ψ_N matrix unchanged - `Q_root` -- tearing eigenvalue(s) in normalized Q * length `nsurfaces` in `:uncoupled` mode * length `1` in `:coupled` mode (global eigenvalue normalized by diff --git a/src/Tearing/Runner/run_slayer.jl b/src/Tearing/Runner/run_slayer.jl index 2a29c63a8..650761b3c 100644 --- a/src/Tearing/Runner/run_slayer.jl +++ b/src/Tearing/Runner/run_slayer.jl @@ -151,6 +151,39 @@ function _build_surface_coupling(model::GGJModel, params::GGJParameters, return surface_coupling(model, params, dp_diag) end +# --------------------------------------------------------------------- +# Reference-length conversion of the outer Δ' for the slab layer +# --------------------------------------------------------------------- +""" + delta_prime_to_rs_reference(dp_matrix, params) -> Matrix{ComplexF64} + +Convert the outer-region Δ' matrix from its ψ_N reference length (the +STRIDE/BVP convention: Frobenius coefficients normalized per unit Δψ_N) to +the r_s-based `x̂ = (r − r_s)/r_s` reference the slab layer works in +(Fitzpatrick 2023 convention — the same reference used by `dc_tmp` and by +the `S^(1/3)` Δ(Q) scale, both built on r_s). + +Near surface `k` the tearing-parity outer solution is +`A_L·|x|^(1/2−μ) + A_S·|x|^(1/2+μ)` with `μ = √(−D_I)` (Glasser, Wang & +Park 2016, Eq. 26). Rescaling the radial variable `x_ψ = K·x̂` with +`K = r_s·(dψ_N/dr)|_s` maps the coefficients as `Â_L = A_L·K^(1/2−μ)` and +`Â_S = A_S·K^(1/2+μ)`, so the response matrix (small coefficient at surface +`i` per unit large coefficient at surface `j`) transforms as + + Δ̂_ij = K_i^(1/2+μ_i) · Δ'_ij · K_j^(μ_j−1/2) + +whose diagonal is `K^(2μ)·Δ'_kk`; at D_I = −1/4 (μ = 1/2) this reduces to +the textbook `Δ̂ = r_s·Δ'_phys`. `K` and `μ` are carried per surface in +`SLAYERParameters.k_ref` / `.mu_mercier`; hand-built parameters default to +`k_ref = 1`, making the conversion the identity. +""" +function delta_prime_to_rs_reference(dp_matrix::AbstractMatrix, + params::AbstractVector) + dl = [p.k_ref^(0.5 + p.mu_mercier) for p in params] + dr = [p.k_ref^(p.mu_mercier - 0.5) for p in params] + return Diagonal(dl) * Matrix{ComplexF64}(dp_matrix) * Diagonal(dr) +end + # --------------------------------------------------------------------- # Core analysis entry point that takes pre-built parameters. # --------------------------------------------------------------------- @@ -194,6 +227,14 @@ function run_slayer_from_inputs(params::AbstractVector{<:InnerLayerParameters}, "(build_slayer_inputs for SLAYER, build_ggj_inputs for GGJ).") ) + # Slab-layer path: convert Δ' from its ψ_N reference length to the + # r_s-based convention shared by the layer Δ(Q) and the critical-Δ (see + # `delta_prime_to_rs_reference`). GGJ is genuinely toroidal/ψ-based + # (its `rescale_delta` handles inner→outer units natively) — no conversion. + if !_is_ggj(model) + dp = delta_prime_to_rs_reference(dp, params) + end + # The coupled determinant uses the reduced m×m (tearing-only) form, which # drops the interchange channel. For GGJ that channel carries the Glasser # interchange stabilization, so coupled-GGJ results omit real physics. diff --git a/test/runtests_slayer_inputs.jl b/test/runtests_slayer_inputs.jl index 55d6124e9..9bc429b54 100644 --- a/test/runtests_slayer_inputs.jl +++ b/test/runtests_slayer_inputs.jl @@ -3,7 +3,7 @@ using GeneralizedPerturbedEquilibrium.Equilibrium using GeneralizedPerturbedEquilibrium.Utilities using GeneralizedPerturbedEquilibrium.InnerLayer - using GeneralizedPerturbedEquilibrium.ForceFreeStates: SingType + using GeneralizedPerturbedEquilibrium.ForceFreeStates: SingType, resist_geometry using TOML # Load the Solovev analytic equilibrium shipped with the examples. @@ -42,8 +42,8 @@ @test r1 > 0 end - @testset "surface_da_dpsi: FD agrees with numerical derivative" begin - # Reference via a tighter FD + @testset "surface_da_dpsi: agrees with FD reference" begin + # The analytic ψ-derivative must reproduce a tight FD of the minor radius. for psi in (0.1, 0.4, 0.7) h_ref = 1e-4 r_p = surface_minor_radius(equil, psi + h_ref) @@ -53,16 +53,24 @@ end end - @testset "surface_da_dpsi: one-sided near boundaries" begin - # Near ψ=0 and ψ=1, the function falls back to one-sided FD and - # should still produce a finite positive number (minor radius is - # still increasing). + @testset "surface_da_dpsi: finite near the boundaries" begin + # The analytic form needs no clamping: near ψ=0 and ψ=1 it still returns a + # finite positive number (large near the axis, where a ~ √ψ). d_near_axis = surface_da_dpsi(equil, 1e-6) d_near_edge = surface_da_dpsi(equil, 1.0 - 1e-6) @test isfinite(d_near_axis) && d_near_axis > 0 @test isfinite(d_near_edge) && d_near_edge > 0 end + @testset "radial_label: analytic derivatives match FD" begin + for rsm in (:midplane, :halfwidth, :fsa, :volume, :flux) + r_at, dr_at = radial_label(equil; rs_method=rsm) + h = 1e-5 + ref = (r_at(0.5 + h) - r_at(0.5 - h)) / (2h) + @test dr_at(0.5) ≈ ref rtol = 1e-3 + end + end + @testset "build_slayer_inputs: returns correct per-surface data" begin sings = [_mk_sing(psi=0.3, q=2.0, q1=1.5, m=2, n=1), _mk_sing(psi=0.6, q=3.0, q1=2.5, m=3, n=1)] @@ -158,6 +166,25 @@ @test sl_var[1].P_perp ≈ sl_s[1].P_perp * 6.0 / 2.0 rtol = 1e-10 end + @testset "build_slayer_inputs: rs_method radial labels are self-consistent" begin + sings = [_mk_sing(psi=0.5, q=2.4, q1=1.2, m=2, n=1)] + got = Dict{Symbol,Any}() + for rsm in (:midplane, :halfwidth, :fsa, :volume, :flux) + sl = build_slayer_inputs(equil, sings, profiles; bt=2.0, dr_val=0.0, rs_method=rsm) + got[rsm] = sl[1] + @test isfinite(sl[1].rs) && sl[1].rs > 0 + @test isfinite(sl[1].k_ref) && sl[1].k_ref > 0 + @test isfinite(sl[1].sval_r) + end + # The outboard-shifted axis compresses the outboard chord, so the shift-free + # half-chord is at least the outboard axis-to-edge distance. + @test got[:halfwidth].rs >= got[:midplane].rs + # Labels genuinely differ (each self-consistent set has its own rs, S, k_ref) + @test got[:fsa].rs != got[:midplane].rs + @test got[:volume].lu != got[:midplane].lu + @test got[:flux].rs != got[:midplane].rs + end + @testset "build_slayer_inputs: dc_type propagates and dr_val activates offset" begin sings = [_mk_sing(psi=0.5, q=2.4, q1=1.2, m=2, n=1)] @@ -178,6 +205,49 @@ @test isfinite(sl_rf[1].dc_tmp) end + @testset "build_slayer_inputs: toroidal dgeo_val derived from ResistGeometry" begin + psi_s, q_s, q1_s = 0.5, 2.4, 1.2 + sing = _mk_sing(psi=psi_s, q=q_s, q1=q1_s, m=2, n=1) + + # Without a ResistGeometry the toroidal factor cannot be derived. + @test_throws ArgumentError build_slayer_inputs(equil, [sing], profiles; + bt=2.0, dc_type=:toroidal, dr_val=0.01) + + sing.restype = resist_geometry(equil, psi_s, q1_s) + p = build_slayer_inputs(equil, [sing], profiles; + bt=2.0, dc_type=:toroidal, dr_val=0.01)[1] + @test isfinite(p.dgeo_val) && p.dgeo_val > 0 + @test isfinite(p.dc_tmp) && p.dc_tmp < 0 + + # Matches a direct evaluation of the r_s-referenced Eq. 59 factor. + rg = sing.restype + rs = surface_minor_radius(equil, psi_s) + k_ref = rs / surface_da_dpsi(equil, psi_s) + dgeo_ref = toroidal_dgeo(; chi1=2π * equil.psio, v1=rg.v1_local, q=q_s, q1=q1_s, n=1, + avg_bsq=rg.avg_bsq, avg_dpsisq=rg.avg_dpsisq, k_ref=k_ref) + @test p.dgeo_val ≈ dgeo_ref rtol = 1e-12 + @test p.k_ref ≈ k_ref rtol = 1e-12 + + # An explicit dgeo_val still overrides the derivation. + p_fix = build_slayer_inputs(equil, [sing], profiles; + bt=2.0, dc_type=:toroidal, dr_val=0.01, dgeo_val=0.3)[1] + @test p_fix.dgeo_val == 0.3 + + # The factor is derived for every dc_type once a ResistGeometry is present. + p_rf = build_slayer_inputs(equil, [sing], profiles; + bt=2.0, dc_type=:rfitzp, dr_val=0.01)[1] + @test p_rf.dgeo_val ≈ dgeo_ref rtol = 1e-12 + + # The radial label enters the factor only through k_ref, the same ratio that + # converts Δ' to the r_s reference, so dgeo/k_ref is label-invariant. + for lab in (:flux, :volume) + p_lab = build_slayer_inputs(equil, [sing], profiles; + bt=2.0, dc_type=:toroidal, dr_val=0.01, rs_method=lab)[1] + @test p_lab.k_ref != p.k_ref + @test p_lab.dgeo_val / p_lab.k_ref ≈ p.dgeo_val / p.k_ref rtol = 1e-10 + end + end + @testset "build_slayer_inputs: empty sings returns empty vector" begin sl = build_slayer_inputs(equil, SingType[], profiles; bt=2.0) @test sl isa Vector{SLAYERParameters} diff --git a/test/runtests_slayer_runner.jl b/test/runtests_slayer_runner.jl index 0eacc18ed..fea52dc46 100644 --- a/test/runtests_slayer_runner.jl +++ b/test/runtests_slayer_runner.jl @@ -153,6 +153,41 @@ @test_throws ArgumentError run_slayer_from_inputs(params, bad_dp, c) end + @testset "delta_prime_to_rs_reference: ψ_N → r_s conversion" begin + # Two surfaces with distinct K and μ; the transform is + # Δ̂_ij = K_i^(1/2+μ_i) · Δ'_ij · K_j^(μ_j−1/2). + K1, mu1 = 0.9, 0.54 + K2, mu2 = 1.2, 0.60 + p1 = SLAYERParameters(; tau=1.0, lu=1e7, c_beta=0.1, D_norm=2.0, + P_perp=20.0, P_tor=10.0, Q_e=-1.0, Q_i=0.5, iota_e=2 / 3, + tauk=1e-4, tau_r=1.0, delta_n=1.0, rs=0.4, R0=1.7, bt=2.0, + sval_r=1.0, eta=2.5e-8, d_beta=4e-3, m=2, n=1, ising=1, + k_ref=K1, mu_mercier=mu1) + p2 = SLAYERParameters(; tau=1.0, lu=1e7, c_beta=0.1, D_norm=2.0, + P_perp=20.0, P_tor=10.0, Q_e=-1.0, Q_i=0.5, iota_e=2 / 3, + tauk=1e-4, tau_r=1.0, delta_n=1.0, rs=0.5, R0=1.7, bt=2.0, + sval_r=1.0, eta=2.5e-8, d_beta=4e-3, m=3, n=1, ising=2, + k_ref=K2, mu_mercier=mu2) + dp = ComplexF64[10.0+1im 2.0-0.5im; 3.0+0im 1.5+2im] + out = Runner.delta_prime_to_rs_reference(dp, [p1, p2]) + # Diagonal carries the scalar K^(2μ) + @test out[1, 1] ≈ K1^(2mu1) * dp[1, 1] + @test out[2, 2] ≈ K2^(2mu2) * dp[2, 2] + # Off-diagonals carry the split row/column factors + @test out[1, 2] ≈ K1^(0.5 + mu1) * K2^(mu2 - 0.5) * dp[1, 2] + @test out[2, 1] ≈ K2^(0.5 + mu2) * K1^(mu1 - 0.5) * dp[2, 1] + # At the slab point μ = 1/2 the diagonal factor is exactly K + pslab = SLAYERParameters(; tau=1.0, lu=1e7, c_beta=0.1, D_norm=2.0, + P_perp=20.0, P_tor=10.0, Q_e=-1.0, Q_i=0.5, iota_e=2 / 3, + tauk=1e-4, tau_r=1.0, delta_n=1.0, rs=0.4, R0=1.7, bt=2.0, + sval_r=1.0, eta=2.5e-8, d_beta=4e-3, m=2, n=1, ising=1, + k_ref=0.8, mu_mercier=0.5) + dp1 = ComplexF64[5.0+0im;;] + @test Runner.delta_prime_to_rs_reference(dp1, [pslab])[1, 1] ≈ 0.8 * 5.0 + # Default parameters (k_ref = 1) give the identity regardless of μ + @test Runner.delta_prime_to_rs_reference(dp, [_mk_params(), _mk_params()]) ≈ dp + end + @testset "run_slayer_from_inputs: coupled mode finds known root" begin # Build a 2-surface problem with a known coupled root by construction. p1 = _mk_params(; rs=0.5, lu=1.0e7, tauk=1.0e-4, Q_e=-1.0, Q_i=0.5, diff --git a/test/runtests_tj_analytic.jl b/test/runtests_tj_analytic.jl index 5bbcb25d2..91d134c8d 100644 --- a/test/runtests_tj_analytic.jl +++ b/test/runtests_tj_analytic.jl @@ -90,4 +90,55 @@ using GeneralizedPerturbedEquilibrium.Equilibrium: TJAnalyticConfig, Equilibrium R_out = R0 + 1.05 # plasma LCFS is at R ≈ R0 + 0.94 @test inp.psi_in((R_out, 0.0)) < 0 end + + @testset "toroidal critical-Δ factor → √(n s r_s/R₀) at ε = 0.05" begin + # Connor et al. 2015 Eq. 59 in the r_s reference must reduce to the + # rfitzp factor on a large-aspect-ratio circular equilibrium. + using GeneralizedPerturbedEquilibrium.ForceFreeStates: resist_geometry + using GeneralizedPerturbedEquilibrium.InnerLayer: toroidal_dgeo, r_based_shear, + surface_minor_radius, surface_da_dpsi + tj = TJAnalyticConfig(lar_r0 = 1.0 / 0.05, lar_a = 1.0, + qc = 1.5, qa = 3.6, pc = 0.001, mu = 2.0, B0 = 12.0, + ma = 64, mtau = 64) + eq = EquilibriumConfig(eq_type = "tj_analytic", + psilow = 0.01, psihigh = 0.995, + mpsi = 64, mtheta = 128, etol = 1e-7) + pe = setup_equilibrium(eq, tj) + chi1 = 2π * pe.psio + for (psi_s, n) in ((0.3, 1), (0.6, 2)) + q = pe.profiles.q_spline(psi_s) + q1 = pe.profiles.q_deriv(psi_s) + rg = resist_geometry(pe, psi_s, q1) + rs = surface_minor_radius(pe, psi_s) + da = surface_da_dpsi(pe, psi_s) + s_r = r_based_shear(rs, q, q1, da) + dgeo = toroidal_dgeo(; chi1=chi1, v1=rg.v1_local, q=q, q1=q1, n=n, + avg_bsq=rg.avg_bsq, avg_dpsisq=rg.avg_dpsisq, k_ref=rs / da) + @test dgeo ≈ sqrt(n * s_r * rs / pe.ro) rtol = 1e-2 + end + end + + @testset "toroidal critical-Δ factor is dimensionless (scale invariance)" begin + # The r_s-referenced Eq. 59 factor must be invariant under B₀ → B₀/2 and + # (a, R₀) → 2(a, R₀) at fixed ε; a missing power of ψ_t' in Λ breaks this. + using GeneralizedPerturbedEquilibrium.ForceFreeStates: resist_geometry + using GeneralizedPerturbedEquilibrium.InnerLayer: toroidal_dgeo, surface_minor_radius, surface_da_dpsi + function _dgeo(a, B0, psi) + tj = TJAnalyticConfig(lar_r0 = a / 0.2, lar_a = a, + qc = 1.5, qa = 3.6, pc = 0.001, mu = 2.0, B0 = B0, + ma = 64, mtau = 64) + eq = EquilibriumConfig(eq_type = "tj_analytic", + psilow = 0.01, psihigh = 0.995, + mpsi = 64, mtheta = 128, etol = 1e-7) + pe = setup_equilibrium(eq, tj) + q1 = pe.profiles.q_deriv(psi) + rg = resist_geometry(pe, psi, q1) + rs = surface_minor_radius(pe, psi) + return toroidal_dgeo(; chi1=2π * pe.psio, v1=rg.v1_local, q=pe.profiles.q_spline(psi), q1=q1, n=1, + avg_bsq=rg.avg_bsq, avg_dpsisq=rg.avg_dpsisq, k_ref=rs / surface_da_dpsi(pe, psi)) + end + d_ref = _dgeo(1.0, 12.0, 0.5) + @test _dgeo(1.0, 6.0, 0.5) ≈ d_ref rtol = 1e-6 + @test _dgeo(2.0, 12.0, 0.5) ≈ d_ref rtol = 1e-6 + end end