From 357116b8ade86d004f3871e1ad61d406883efa8f Mon Sep 17 00:00:00 2001 From: d-burg Date: Thu, 3 Sep 2026 13:02:50 -0400 Subject: [PATCH 01/16] =?UTF-8?q?InnerLayer.SLAYER=20-=20FEATURE!=20-=20De?= =?UTF-8?q?rive=20the=20toroidal=20critical-=CE=94=20geometric=20factor=20?= =?UTF-8?q?in=20the=20r=5Fs=20reference?= MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit The Connor et al. 2015 (PPCF 57 065001) Eq. 59 geometric factor of the toroidal critical-Δ (dc_type=:toroidal) is now derived per surface from the equilibrium instead of requiring a user-supplied dgeo_val. `toroidal_dgeo` evaluates V_s·(α²Λ²/(⟨B²⟩⟨|∇V|²⟩))^{1/4} with Λ = ψ_t'² ι'/2π and converts it 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 Δ'; the factor r_s(dV/dr)/V_s = k_ref·v1/V_s cancels V_s. ResistGeometry gains the ⟨|∇ψ_N|²⟩ average it needs. At large aspect ratio the factor reduces to √(n s r_s/R₀), so :toroidal coincides with :rfitzp there; the factor is dimensionless and scales with the radial label exactly as the converted Δ' does. The PerSurface/D_geo dataset now carries derived values instead of zeros. Tests: derivation wiring and label invariance (LayerInputs), large-aspect-ratio limit and scale invariance (TJ analytic). Benchmarks regenerate the convergence, verification, label-sensitivity and TJ ε-scan figures. Co-Authored-By: Claude Fable 5.1 --- .../benchmark_delta_crit_radial_label.jl | 168 ++++++++++++++++ benchmarks/benchmark_toroidal_delta_crit.jl | 183 ++++++++++++++++++ .../verify_toroidal_delta_crit_symbolic.py | 109 +++++++++++ src/ForceFreeStates/Surfaces/ResistEval.jl | 11 +- src/InnerLayer/InnerLayer.jl | 4 +- src/InnerLayer/SLAYER/LayerInputs.jl | 94 ++++++--- src/InnerLayer/SLAYER/LayerParameters.jl | 8 +- src/InnerLayer/SLAYER/SLAYER.jl | 2 +- src/Tearing/Runner/Control.jl | 3 +- src/Tearing/Runner/HDF5Output.jl | 2 +- test/runtests_slayer_inputs.jl | 45 ++++- test/runtests_tj_analytic.jl | 51 +++++ 12 files changed, 639 insertions(+), 41 deletions(-) create mode 100644 benchmarks/benchmark_delta_crit_radial_label.jl create mode 100644 benchmarks/benchmark_toroidal_delta_crit.jl create mode 100644 benchmarks/verify_toroidal_delta_crit_symbolic.py diff --git a/benchmarks/benchmark_delta_crit_radial_label.jl b/benchmarks/benchmark_delta_crit_radial_label.jl new file mode 100644 index 000000000..1361d8363 --- /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.alpha_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/benchmarks/verify_toroidal_delta_crit_symbolic.py b/benchmarks/verify_toroidal_delta_crit_symbolic.py new file mode 100644 index 000000000..eef0cda8a --- /dev/null +++ b/benchmarks/verify_toroidal_delta_crit_symbolic.py @@ -0,0 +1,109 @@ +"""Symbolic verification of the r_s-referenced Connor et al. 2015 Eq. 59 toroidal critical-Δ factor. + +Checks, with SymPy (run: uv run --with sympy python3 verify_toroidal_delta_crit_symbolic.py): + 1. Λ ≡ ψ'χ'' − χ'ψ'' = ψ'² (ι/2π)' (Connor's two definitions agree) + 2. the GPEC expressions for ψ_t', (ι/2π)', Λ, α in terms of (chi1, v1, q, q1) are the + chain-rule images of the V-derivative definitions + 3. the code formula k_ref·v1·(α²Λ²/(⟨B²⟩ v1²⟨|∇ψ_N|²⟩))^{1/4} equals + V_s (α²Λ²/(⟨B²⟩⟨|∇V|²⟩))^{1/4} · r_s(dV/dr)/V_s (V_s cancels) + 4. circular large-aspect-ratio limit: Eq. 59's geometric factor → ½√(n s r/R), and the + r_s-referenced factor → √(n s r/R), for an arbitrary q(r) + 5. prefactor chain: Eq. 61 · r_s ≡ the rfitzp code formula −√2π^{3/2}D_R/W_d with Connor Eq. 65 + 6. scaling: the r_s-referenced factor is invariant under B → λB and lengths → μ·lengths; the + Fortran STRIDE form (one power of ψ_t' in Λ) scales as (λ/μ)^{-1/2} +""" +import sympy as sp + +ok = True +def check(name, expr): + global ok + res = sp.simplify(expr) + passed = res == 0 + ok &= passed + print(f"[{'PASS' if passed else 'FAIL'}] {name}" + ("" if passed else f" residual = {res}")) + +# ---------------------------------------------------------------- 1. Λ identity +V = sp.symbols('V', positive=True) +psi = sp.Function('psi')(V) # toroidal flux ψ(V) +chi = sp.Function('chi')(V) # poloidal flux χ(V) +iota2pi = sp.diff(chi, V) / sp.diff(psi, V) +Lam_def = sp.diff(psi, V) * sp.diff(chi, V, 2) - sp.diff(chi, V) * sp.diff(psi, V, 2) +check("1. psi'chi'' - chi'psi'' == psi'^2 (iota/2pi)'", Lam_def - sp.diff(psi, V)**2 * sp.diff(iota2pi, V)) + +# ---------------------------------------------------------------- 2. GPEC chain rule +x = sp.symbols('psi_N', positive=True) # normalized poloidal flux +chi1, n = sp.symbols('chi1 n', positive=True) # chi1 = dχ/dψ_N = 2π psio ; toroidal mode number +psit = sp.Function('psi_t')(x) # toroidal flux ψ_t(ψ_N) +Vf = sp.Function('V')(x) # volume V(ψ_N) +v1 = sp.diff(Vf, x) # dV/dψ_N +q = sp.diff(psit, x) / chi1 # q = dψ_t/dχ +q1 = sp.diff(q, x) # dq/dψ_N +# V-derivatives via d/dV = (1/v1) d/dψ_N +psit_V = sp.diff(psit, x) / v1 +iota2pi_V = sp.diff(1 / q, x) / v1 +Lam_V = psit_V**2 * iota2pi_V +alpha_V = 2 * sp.pi * n / (chi1 / v1) # α = 2πn/χ' +# code expressions (LayerInputs.jl toroidal_dgeo) +psit1_code = q * chi1 / v1 +Lam_code = psit1_code**2 * (-q1 / (q**2 * v1)) +alpha_code = 2 * sp.pi * n * v1 / chi1 +check("2a. psi_t' code == chain rule", psit1_code - psit_V) +check("2b. Lambda code == psi_t'^2 (iota/2pi)'", Lam_code - Lam_V) +check("2c. alpha code == 2 pi n / chi'", alpha_code - alpha_V) + +# ---------------------------------------------------------------- 3. reference conversion, V_s cancels +Vs, rs, dadpsi, B2, G, v1s, a2, L2 = sp.symbols('V_s r_s da_dpsi B2 G v1 alpha2 Lambda2', positive=True) +eq59_geo = Vs * (a2 * L2 / (B2 * (v1s**2 * G)))**sp.Rational(1, 4) # ⟨|∇V|²⟩ = v1²⟨|∇ψ_N|²⟩ +dVdr = v1s / dadpsi # dV/dr = (dV/dψ_N)/(da/dψ_N) +conversion = rs * dVdr / Vs # Y=(V−V_s)/V_s → x̂=(r−r_s)/r_s +k_ref = rs / dadpsi +dgeo_code = k_ref * v1s * (a2 * L2 / (B2 * v1s**2 * G))**sp.Rational(1, 4) +check("3. code dgeo == Eq.59 factor x r_s (dV/dr)/V_s (V_s cancels)", dgeo_code - eq59_geo * conversion) + +# ---------------------------------------------------------------- 4. circular LAR limit, arbitrary q(r) +r, R, B = sp.symbols('r R B', positive=True) +qr = sp.Function('q')(r) +V_r = 2 * sp.pi**2 * R * r**2 # volume +psit_r = sp.pi * r**2 * B # toroidal flux (uniform B) +dchi_dr = 2 * sp.pi * r * B / qr # dχ/dr = dψ_t/dr / q (q = dψ_t/dχ) +dVdr_r = sp.diff(V_r, r) +d_dV = lambda f: sp.diff(f, r) / dVdr_r +psit_V_r = d_dV(psit_r) +chi_V_r = dchi_dr / dVdr_r +Lam_r = psit_V_r**2 * d_dV(chi_V_r / psit_V_r) # ψ'^2 (χ'/ψ')' = ψ'^2 (ι/2π)' +alpha_r = 2 * sp.pi * n / chi_V_r +B2_r = B**2 # ⟨B²⟩ at leading order +gradV2_r = dVdr_r**2 # ⟨|∇V|²⟩ = (dV/dr)² +s_r = r * sp.diff(qr, r) / qr # r-based shear +eq59_r = V_r * (alpha_r**2 * Lam_r**2 / (B2_r * gradV2_r))**sp.Rational(1, 4) +lar_half = sp.Rational(1, 2) * sp.sqrt(n * s_r * r / R) +# q' may be negative: compare squares of both sides (both sides positive for s>0) and the sign +check("4a. Eq.59 geometric factor (Y ref) -> 1/2 sqrt(n s r/R)", sp.powsimp(eq59_r**4 - lar_half**4, force=True)) +dgeo_r = eq59_r * r * dVdr_r / V_r +check("4b. r_s-referenced factor -> sqrt(n s r/R)", sp.powsimp(dgeo_r**4 - (n * s_r * r / R)**2, force=True)) +check("4c. conversion r_s (dV/dr)/V_s == 2 on a circle", r * dVdr_r / V_r - 2) + +# ---------------------------------------------------------------- 5. prefactor chain: Eq.61 · r_s == rfitzp +chipar, chiperp, DR, s = sp.symbols('chi_par chi_perp D_R s', positive=True) +eq61 = sp.pi**sp.Rational(3, 2) / 2 * (chipar / chiperp)**sp.Rational(1, 4) * sp.sqrt(n * s / (R * r)) * (-DR) +Wd = sp.sqrt(8) * (chiperp / chipar)**sp.Rational(1, 4) / sp.sqrt(r / R * s * n) # Connor Eq. 65, W_d/r_s +rfitzp = -sp.sqrt(2) * sp.pi**sp.Rational(3, 2) * DR / Wd # LayerParameters.jl :rfitzp +check("5. Eq.61 * r_s == rfitzp code formula", eq61 * r - rfitzp) +# and the toroidal branch with dgeo -> sqrt(n s r/R) equals rfitzp +toroidal_lar = sp.Rational(1, 2) * (-DR) * sp.pi**sp.Rational(3, 2) * (chipar / chiperp)**sp.Rational(1, 4) * sp.sqrt(n * s * r / R) +check("5b. toroidal branch at LAR == rfitzp", toroidal_lar - rfitzp) + +# ---------------------------------------------------------------- 6. scaling +# q is a shape function of r/R, so it is invariant under a uniform length scaling. +lam, mu = sp.symbols('lambda mu', positive=True) +qshape = sp.Function('q')(r / R) +dgeo_shape = dgeo_r.subs(qr, qshape).doit() +fortran_shape = (dgeo_r / sp.sqrt(psit_V_r)).subs(qr, qshape).doit() +scale = {B: lam * B, r: mu * r, R: mu * R} +ratio_code = sp.simplify((dgeo_shape.subs(scale, simultaneous=True) / dgeo_shape)**4) +ratio_fortran = sp.simplify((fortran_shape.subs(scale, simultaneous=True) / fortran_shape)**2) +check("6a. r_s-referenced factor invariant under B->lambda B, lengths->mu lengths", ratio_code - 1) +check("6b. Fortran form (psi_t'^1 in Lambda) scales as (lambda/mu)^(-1/2)", ratio_fortran - mu / lam) + +print("\nALL PASS" if ok else "\nSOME CHECKS FAILED") +raise SystemExit(0 if ok else 1) diff --git a/src/ForceFreeStates/Surfaces/ResistEval.jl b/src/ForceFreeStates/Surfaces/ResistEval.jl index 04953b6c0..7719c4f5b 100644 --- a/src/ForceFreeStates/Surfaces/ResistEval.jl +++ b/src/ForceFreeStates/Surfaces/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 ddbde66c1..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, radial_label, 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, radial_label, 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 252e81ae9..0ab8885db 100644 --- a/src/InnerLayer/SLAYER/LayerInputs.jl +++ b/src/InnerLayer/SLAYER/LayerInputs.jl @@ -154,6 +154,40 @@ function radial_label(equil; rs_method::Symbol=:midplane, theta::Real=0.0) 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 + """ build_slayer_inputs(equil, sings, profiles; …) -> Vector{SLAYERParameters} @@ -199,13 +233,14 @@ 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 @@ -332,27 +367,6 @@ 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. - dgeo_val_k = if dgeo_val === nothing - 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." - ) - ) - 0.0 - else - _eval(dgeo_val, psi) - end - # Reference-length conversion inputs for the outer Δ': K = r_s·(dψ_N/dr)|_s and # α = √(−D_I) (Glasser-Greene-Johnson 1975 Eq. 48), with α clamped to 0 on Mercier-unstable # surfaces (the factor turns complex there) and K = 1 whenever da/dψ is not a usable @@ -364,6 +378,30 @@ function build_slayer_inputs(equil, sings, profiles::KineticProfiles; "Δ' unconverted (k_ref = 1) at this surface.", maxlog=3) 1.0 end + + # 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 + 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 with " * + "dgeo_val=nothing requires `sing.restype` populated " * + "by ForceFreeStates.resist_eval_all!. " * + "Surface k=$k has restype=nothing." + ) + ) + else + 0.0 + end + else + _eval(dgeo_val, psi) + end + alpha_k = if rg === nothing @warn("build_slayer_inputs: sing.restype not populated; using the " * "slab Mercier exponent α = 1/2 for the Δ' reference-length " * diff --git a/src/InnerLayer/SLAYER/LayerParameters.jl b/src/InnerLayer/SLAYER/LayerParameters.jl index ae875dea4..705ff50b3 100644 --- a/src/InnerLayer/SLAYER/LayerParameters.jl +++ b/src/InnerLayer/SLAYER/LayerParameters.jl @@ -41,7 +41,7 @@ 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 | @@ -130,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, @@ -212,7 +212,9 @@ 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, diff --git a/src/InnerLayer/SLAYER/SLAYER.jl b/src/InnerLayer/SLAYER/SLAYER.jl index 033f71eca..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, radial_label, 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/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 b0fc8d2e4..de8f4906e 100644 --- a/src/Tearing/Runner/HDF5Output.jl +++ b/src/Tearing/Runner/HDF5Output.jl @@ -90,7 +90,7 @@ 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",)), diff --git a/test/runtests_slayer_inputs.jl b/test/runtests_slayer_inputs.jl index fce001036..8b431884a 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. @@ -205,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_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 From 9b9efae8972895f1a113d27de7ecfb97cd7453a4 Mon Sep 17 00:00:00 2001 From: d-burg Date: Mon, 14 Sep 2026 16:10:50 -0400 Subject: [PATCH 02/16] =?UTF-8?q?InnerLayer.SLAYER=20-=20DOCS=20-=20Record?= =?UTF-8?q?=20the=20linear=20reference-conversion=20decision=20for=20the?= =?UTF-8?q?=20toroidal=20critical-=CE=94?= MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Co-Authored-By: Claude Fable 5.1 --- src/InnerLayer/SLAYER/LayerInputs.jl | 7 +++++++ 1 file changed, 7 insertions(+) diff --git a/src/InnerLayer/SLAYER/LayerInputs.jl b/src/InnerLayer/SLAYER/LayerInputs.jl index 0ab8885db..9d0d35d5f 100644 --- a/src/InnerLayer/SLAYER/LayerInputs.jl +++ b/src/InnerLayer/SLAYER/LayerInputs.jl @@ -176,6 +176,13 @@ normalized poloidal flux, and `'` = d/dV: 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`. + +The reference conversion is deliberately linear (Frobenius exponent ½, the +paper's H = 0 ordering): the critical-Δ is a layer-side quantity like the slab +Δ̂(Q) and W_d, whereas the outer Δ' is converted with `K^(2α)` because it +genuinely carries the Mercier exponent α. Converting the offset with `c^(2α)` +instead would multiply it by `c^(2α−1)`, 5–13 % at the DIII-D-like 3/2–4/1 +surfaces, without support from the derivation. """ function toroidal_dgeo(; chi1::Real, v1::Real, q::Real, q1::Real, n::Integer, avg_bsq::Real, avg_dpsisq::Real, k_ref::Real) From f13fea683282cf47a38a6cb4a94b54033f4179cf Mon Sep 17 00:00:00 2001 From: d-burg Date: Wed, 23 Sep 2026 16:27:48 -0400 Subject: [PATCH 03/16] =?UTF-8?q?Benchmarks=20-=20MINOR=20-=20Move=20the?= =?UTF-8?q?=20toroidal=20critical-=CE=94=20study=20scripts=20to=20CTM-proc?= =?UTF-8?q?essing?= MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit They are one-off physics studies, archived with their figures under CTM-processing/julia_deltaprime_harness/figures/delta_crit_toroidal/. Co-Authored-By: Claude Opus 5.5 --- .../benchmark_delta_crit_radial_label.jl | 168 ---------------- benchmarks/benchmark_toroidal_delta_crit.jl | 183 ------------------ .../verify_toroidal_delta_crit_symbolic.py | 109 ----------- 3 files changed, 460 deletions(-) delete mode 100644 benchmarks/benchmark_delta_crit_radial_label.jl delete mode 100644 benchmarks/benchmark_toroidal_delta_crit.jl delete mode 100644 benchmarks/verify_toroidal_delta_crit_symbolic.py diff --git a/benchmarks/benchmark_delta_crit_radial_label.jl b/benchmarks/benchmark_delta_crit_radial_label.jl deleted file mode 100644 index 1361d8363..000000000 --- a/benchmarks/benchmark_delta_crit_radial_label.jl +++ /dev/null @@ -1,168 +0,0 @@ -# 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.alpha_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 deleted file mode 100644 index 2b5880aa7..000000000 --- a/benchmarks/benchmark_toroidal_delta_crit.jl +++ /dev/null @@ -1,183 +0,0 @@ -# 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/benchmarks/verify_toroidal_delta_crit_symbolic.py b/benchmarks/verify_toroidal_delta_crit_symbolic.py deleted file mode 100644 index eef0cda8a..000000000 --- a/benchmarks/verify_toroidal_delta_crit_symbolic.py +++ /dev/null @@ -1,109 +0,0 @@ -"""Symbolic verification of the r_s-referenced Connor et al. 2015 Eq. 59 toroidal critical-Δ factor. - -Checks, with SymPy (run: uv run --with sympy python3 verify_toroidal_delta_crit_symbolic.py): - 1. Λ ≡ ψ'χ'' − χ'ψ'' = ψ'² (ι/2π)' (Connor's two definitions agree) - 2. the GPEC expressions for ψ_t', (ι/2π)', Λ, α in terms of (chi1, v1, q, q1) are the - chain-rule images of the V-derivative definitions - 3. the code formula k_ref·v1·(α²Λ²/(⟨B²⟩ v1²⟨|∇ψ_N|²⟩))^{1/4} equals - V_s (α²Λ²/(⟨B²⟩⟨|∇V|²⟩))^{1/4} · r_s(dV/dr)/V_s (V_s cancels) - 4. circular large-aspect-ratio limit: Eq. 59's geometric factor → ½√(n s r/R), and the - r_s-referenced factor → √(n s r/R), for an arbitrary q(r) - 5. prefactor chain: Eq. 61 · r_s ≡ the rfitzp code formula −√2π^{3/2}D_R/W_d with Connor Eq. 65 - 6. scaling: the r_s-referenced factor is invariant under B → λB and lengths → μ·lengths; the - Fortran STRIDE form (one power of ψ_t' in Λ) scales as (λ/μ)^{-1/2} -""" -import sympy as sp - -ok = True -def check(name, expr): - global ok - res = sp.simplify(expr) - passed = res == 0 - ok &= passed - print(f"[{'PASS' if passed else 'FAIL'}] {name}" + ("" if passed else f" residual = {res}")) - -# ---------------------------------------------------------------- 1. Λ identity -V = sp.symbols('V', positive=True) -psi = sp.Function('psi')(V) # toroidal flux ψ(V) -chi = sp.Function('chi')(V) # poloidal flux χ(V) -iota2pi = sp.diff(chi, V) / sp.diff(psi, V) -Lam_def = sp.diff(psi, V) * sp.diff(chi, V, 2) - sp.diff(chi, V) * sp.diff(psi, V, 2) -check("1. psi'chi'' - chi'psi'' == psi'^2 (iota/2pi)'", Lam_def - sp.diff(psi, V)**2 * sp.diff(iota2pi, V)) - -# ---------------------------------------------------------------- 2. GPEC chain rule -x = sp.symbols('psi_N', positive=True) # normalized poloidal flux -chi1, n = sp.symbols('chi1 n', positive=True) # chi1 = dχ/dψ_N = 2π psio ; toroidal mode number -psit = sp.Function('psi_t')(x) # toroidal flux ψ_t(ψ_N) -Vf = sp.Function('V')(x) # volume V(ψ_N) -v1 = sp.diff(Vf, x) # dV/dψ_N -q = sp.diff(psit, x) / chi1 # q = dψ_t/dχ -q1 = sp.diff(q, x) # dq/dψ_N -# V-derivatives via d/dV = (1/v1) d/dψ_N -psit_V = sp.diff(psit, x) / v1 -iota2pi_V = sp.diff(1 / q, x) / v1 -Lam_V = psit_V**2 * iota2pi_V -alpha_V = 2 * sp.pi * n / (chi1 / v1) # α = 2πn/χ' -# code expressions (LayerInputs.jl toroidal_dgeo) -psit1_code = q * chi1 / v1 -Lam_code = psit1_code**2 * (-q1 / (q**2 * v1)) -alpha_code = 2 * sp.pi * n * v1 / chi1 -check("2a. psi_t' code == chain rule", psit1_code - psit_V) -check("2b. Lambda code == psi_t'^2 (iota/2pi)'", Lam_code - Lam_V) -check("2c. alpha code == 2 pi n / chi'", alpha_code - alpha_V) - -# ---------------------------------------------------------------- 3. reference conversion, V_s cancels -Vs, rs, dadpsi, B2, G, v1s, a2, L2 = sp.symbols('V_s r_s da_dpsi B2 G v1 alpha2 Lambda2', positive=True) -eq59_geo = Vs * (a2 * L2 / (B2 * (v1s**2 * G)))**sp.Rational(1, 4) # ⟨|∇V|²⟩ = v1²⟨|∇ψ_N|²⟩ -dVdr = v1s / dadpsi # dV/dr = (dV/dψ_N)/(da/dψ_N) -conversion = rs * dVdr / Vs # Y=(V−V_s)/V_s → x̂=(r−r_s)/r_s -k_ref = rs / dadpsi -dgeo_code = k_ref * v1s * (a2 * L2 / (B2 * v1s**2 * G))**sp.Rational(1, 4) -check("3. code dgeo == Eq.59 factor x r_s (dV/dr)/V_s (V_s cancels)", dgeo_code - eq59_geo * conversion) - -# ---------------------------------------------------------------- 4. circular LAR limit, arbitrary q(r) -r, R, B = sp.symbols('r R B', positive=True) -qr = sp.Function('q')(r) -V_r = 2 * sp.pi**2 * R * r**2 # volume -psit_r = sp.pi * r**2 * B # toroidal flux (uniform B) -dchi_dr = 2 * sp.pi * r * B / qr # dχ/dr = dψ_t/dr / q (q = dψ_t/dχ) -dVdr_r = sp.diff(V_r, r) -d_dV = lambda f: sp.diff(f, r) / dVdr_r -psit_V_r = d_dV(psit_r) -chi_V_r = dchi_dr / dVdr_r -Lam_r = psit_V_r**2 * d_dV(chi_V_r / psit_V_r) # ψ'^2 (χ'/ψ')' = ψ'^2 (ι/2π)' -alpha_r = 2 * sp.pi * n / chi_V_r -B2_r = B**2 # ⟨B²⟩ at leading order -gradV2_r = dVdr_r**2 # ⟨|∇V|²⟩ = (dV/dr)² -s_r = r * sp.diff(qr, r) / qr # r-based shear -eq59_r = V_r * (alpha_r**2 * Lam_r**2 / (B2_r * gradV2_r))**sp.Rational(1, 4) -lar_half = sp.Rational(1, 2) * sp.sqrt(n * s_r * r / R) -# q' may be negative: compare squares of both sides (both sides positive for s>0) and the sign -check("4a. Eq.59 geometric factor (Y ref) -> 1/2 sqrt(n s r/R)", sp.powsimp(eq59_r**4 - lar_half**4, force=True)) -dgeo_r = eq59_r * r * dVdr_r / V_r -check("4b. r_s-referenced factor -> sqrt(n s r/R)", sp.powsimp(dgeo_r**4 - (n * s_r * r / R)**2, force=True)) -check("4c. conversion r_s (dV/dr)/V_s == 2 on a circle", r * dVdr_r / V_r - 2) - -# ---------------------------------------------------------------- 5. prefactor chain: Eq.61 · r_s == rfitzp -chipar, chiperp, DR, s = sp.symbols('chi_par chi_perp D_R s', positive=True) -eq61 = sp.pi**sp.Rational(3, 2) / 2 * (chipar / chiperp)**sp.Rational(1, 4) * sp.sqrt(n * s / (R * r)) * (-DR) -Wd = sp.sqrt(8) * (chiperp / chipar)**sp.Rational(1, 4) / sp.sqrt(r / R * s * n) # Connor Eq. 65, W_d/r_s -rfitzp = -sp.sqrt(2) * sp.pi**sp.Rational(3, 2) * DR / Wd # LayerParameters.jl :rfitzp -check("5. Eq.61 * r_s == rfitzp code formula", eq61 * r - rfitzp) -# and the toroidal branch with dgeo -> sqrt(n s r/R) equals rfitzp -toroidal_lar = sp.Rational(1, 2) * (-DR) * sp.pi**sp.Rational(3, 2) * (chipar / chiperp)**sp.Rational(1, 4) * sp.sqrt(n * s * r / R) -check("5b. toroidal branch at LAR == rfitzp", toroidal_lar - rfitzp) - -# ---------------------------------------------------------------- 6. scaling -# q is a shape function of r/R, so it is invariant under a uniform length scaling. -lam, mu = sp.symbols('lambda mu', positive=True) -qshape = sp.Function('q')(r / R) -dgeo_shape = dgeo_r.subs(qr, qshape).doit() -fortran_shape = (dgeo_r / sp.sqrt(psit_V_r)).subs(qr, qshape).doit() -scale = {B: lam * B, r: mu * r, R: mu * R} -ratio_code = sp.simplify((dgeo_shape.subs(scale, simultaneous=True) / dgeo_shape)**4) -ratio_fortran = sp.simplify((fortran_shape.subs(scale, simultaneous=True) / fortran_shape)**2) -check("6a. r_s-referenced factor invariant under B->lambda B, lengths->mu lengths", ratio_code - 1) -check("6b. Fortran form (psi_t'^1 in Lambda) scales as (lambda/mu)^(-1/2)", ratio_fortran - mu / lam) - -print("\nALL PASS" if ok else "\nSOME CHECKS FAILED") -raise SystemExit(0 if ok else 1) From 7b38745c778a9e66424bc624ba9574df5956cb09 Mon Sep 17 00:00:00 2001 From: d-burg Date: Wed, 23 Sep 2026 16:27:48 -0400 Subject: [PATCH 04/16] =?UTF-8?q?InnerLayer.SLAYER=20-=20MINOR=20-=20Corre?= =?UTF-8?q?ct=20the=20critical-=CE=94=20citation,=20document=20the=20dgeo?= =?UTF-8?q?=5Fval=20reference,=20and=20tighten=20its=20guards?= MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit - Cite Connor, Ham, Hastie & Liu 2015 (PPCF 57 065001) everywhere the source was credited to "Connor-Hastie-Helander". - A prescribed dgeo_val must be in the r_s reference (Eq. 59's value times k_ref·v1/V_s). - toroidal_dgeo requires dV/dψ_N > 0 and a non-zero chi1; the q guard could not fire. - Name the Mercier exponent α_M so it no longer collides with Connor's α, and drop the equilibrium-specific percentages from the docstring. - The k_ref fallback warning names the critical-Δ too; ResistGeometry docs name ψ_N and the eight integrands. Co-Authored-By: Claude Opus 5.5 --- src/ForceFreeStates/Surfaces/ResistEval.jl | 93 +++++++++++----------- src/InnerLayer/SLAYER/LayerInputs.jl | 33 ++++---- src/Tearing/Runner/Control.jl | 2 +- src/Tearing/Runner/HDF5Output.jl | 12 ++- 4 files changed, 74 insertions(+), 66 deletions(-) diff --git a/src/ForceFreeStates/Surfaces/ResistEval.jl b/src/ForceFreeStates/Surfaces/ResistEval.jl index 7719c4f5b..d653357ea 100644 --- a/src/ForceFreeStates/Surfaces/ResistEval.jl +++ b/src/ForceFreeStates/Surfaces/ResistEval.jl @@ -33,24 +33,24 @@ Per-singular-surface Glasser-Greene-Johnson geometric coefficients and supporting flux-surface averages. -| field | meaning | -|-------------|------------------------------------------------------| -| `E`, `F` | Glasser interchange parameters (enter `D_I = E+F+H-¼`) | -| `G` | Coupling coefficient (curvature × pressure gradient) | -| `H` | Pfirsch-Schlüter coefficient | -| `K` | Glasser parameter | -| `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 | -| `R_major` | flux-surface-averaged major radius ⟨R⟩ [m] | -| `eps_local` | (R_max − R_min)/2 / R_major — local inverse aspect ratio | -| `p_local` | Plasma pressure at this surface [Pa] | -| `p1_local` | dp/dψ at this surface | -| `v1_local` | dV/dψ at this surface | +| field | meaning | +|:--------------------- |:-------------------------------------------------------- | +| `E`, `F` | Glasser interchange parameters (enter `D_I = E+F+H-¼`) | +| `G` | Coupling coefficient (curvature × pressure gradient) | +| `H` | Pfirsch-Schlüter coefficient | +| `K` | Glasser parameter | +| `M` | Mass factor | +| `avg_bsq_over_dpsisq` | ⟨B²/ | +| `avg_bsq` | ⟨B²⟩ — needed for τ_R | +| `avg_dpsisq` | ⟨ | +| `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 | +| `R_major` | flux-surface-averaged major radius ⟨R⟩ [m] | +| `eps_local` | (R_max − R_min)/2 / R_major — local inverse aspect ratio | +| `p_local` | Plasma pressure at this surface [Pa] | +| `p1_local` | dp/dψ at this surface | +| `v1_local` | dV/dψ at this surface | `H` here is identical to the `H` reported by `mercier_scan!` and stored in `LocalStability/h` — the GGJ routine recomputes it for convenience. @@ -86,9 +86,9 @@ end resist_geometry(equil, psifac, q1; gamma=5/3) -> ResistGeometry Port of Fortran RDCON `resist_eval` restricted to the -pure-equilibrium geometric coefficients. Integrates the 6 theta integrands -at the given flux surface and combines them into E, F, G, H, K, M via the -standard GGJ formulas. +pure-equilibrium geometric coefficients. Integrates the 8 theta integrands +at the given flux surface and combines the 6 GGJ ones into E, F, G, H, K, M via +the standard GGJ formulas. # Arguments @@ -101,6 +101,7 @@ standard GGJ formulas. - `gamma` — adiabatic index (default 5/3) !!! note "Contract" + `psifac` must be a genuine interior rational surface (`0 < ψ < 1`) with nonzero `q1`, `p1 = dp/dψ`, and `p`. The GGJ combination divides by these and by `|∇ψ|²` (which → 0 at the axis), so calling on the magnetic axis, @@ -109,47 +110,47 @@ standard GGJ formulas. rationals. """ function resist_geometry(equil::Equilibrium.PlasmaEquilibrium, - psifac::Real, q1::Real; gamma::Real=5/3) + psifac::Real, q1::Real; gamma::Real=5 / 3) profiles = equil.profiles - twopi = 2π - chi1 = twopi * equil.psio - psi_f = Float64(psifac) + twopi = 2π + chi1 = twopi * equil.psio + psi_f = Float64(psifac) # Surface-profile quantities (evaluate via the existing splines) twopif = profiles.F_spline(psi_f) - p = profiles.P_spline(psi_f) - p1 = profiles.P_deriv(psi_f) - v1 = profiles.dVdpsi_spline(psi_f) - v2 = profiles.dVdpsi_deriv(psi_f) - q = profiles.q_spline(psi_f) + p = profiles.P_spline(psi_f) + p1 = profiles.P_deriv(psi_f) + v1 = profiles.dVdpsi_spline(psi_f) + v2 = profiles.dVdpsi_deriv(psi_f) + q = profiles.q_spline(psi_f) # 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, 8) - B_max = -Inf - B_min = Inf - R_max = -Inf - R_min = Inf + ff = zeros(Float64, ntheta, 8) + B_max = -Inf + B_min = Inf + R_max = -Inf + R_min = Inf for itheta in 1:ntheta theta = equil.rzphi_ys[itheta] - f1 = equil.rzphi_rsquared((psi_f, theta)) - f2 = equil.rzphi_offset((psi_f, theta)) + f1 = equil.rzphi_rsquared((psi_f, theta)) + f2 = equil.rzphi_offset((psi_f, theta)) jac = equil.rzphi_jac((psi_f, theta)) fy1 = FastInterpolations.deriv_view(equil.rzphi_rsquared, (0, 1))((psi_f, theta)) - fy2 = FastInterpolations.deriv_view(equil.rzphi_offset, (0, 1))((psi_f, theta)) - fy3 = FastInterpolations.deriv_view(equil.rzphi_nu, (0, 1))((psi_f, theta)) + fy2 = FastInterpolations.deriv_view(equil.rzphi_offset, (0, 1))((psi_f, theta)) + fy3 = FastInterpolations.deriv_view(equil.rzphi_nu, (0, 1))((psi_f, theta)) rfac = sqrt(f1) - eta = twopi * (theta + f2) - r = equil.ro + rfac * cos(eta) + eta = twopi * (theta + f2) + r = equil.ro + rfac * cos(eta) v21 = fy1 / (2 * rfac * jac) v22 = (1 + fy2) * twopi * rfac / jac v23 = fy3 * r / jac v33 = twopi * r / jac - bsq = chi1^2 * (v21^2 + v22^2 + (v23 + q*v33)^2) + bsq = chi1^2 * (v21^2 + v22^2 + (v23 + q * v33)^2) dpsisq = (twopi * r)^2 * (v21^2 + v22^2) B_here = sqrt(bsq) @@ -183,7 +184,7 @@ function resist_geometry(equil::Equilibrium.PlasmaEquilibrium, (twopif * q1 * chi1 / avg[5] - v2) F_coef = (p1 * v1 / (q1 * chi1^2))^2 * (avg[1] * avg[3] + (twopif / chi1)^2 * - (avg[1] * avg[4] - avg[2]^2)) + (avg[1] * avg[4] - avg[2]^2)) H_coef = twopif * p1 * v1 / (q1 * chi1^3) * (avg[2] - avg[1] / avg[5]) M_coef = avg[1] * (avg[6] + (twopif / chi1)^2 * (avg[3] - 1.0 / avg[5])) @@ -195,7 +196,7 @@ function resist_geometry(equil::Equilibrium.PlasmaEquilibrium, E_coef, F_coef, G_coef, H_coef, K_coef, M_coef, avg[1], avg[5], avg[8], avg_B, B_max, B_min, f_trap, R_major, eps_local, - p, p1, v1, + p, p1, v1 ) end @@ -207,8 +208,8 @@ Populate `sing.restype` for every `SingType` in `intr.sing` using filled. """ function resist_eval_all!(intr::ForceFreeStatesInternal, - equil::Equilibrium.PlasmaEquilibrium; - gamma::Real=5/3) + equil::Equilibrium.PlasmaEquilibrium; + gamma::Real=5 / 3) for sing in intr.sing sing.restype === nothing || continue sing.restype = resist_geometry(equil, sing.psifac, sing.q1; gamma=gamma) diff --git a/src/InnerLayer/SLAYER/LayerInputs.jl b/src/InnerLayer/SLAYER/LayerInputs.jl index 9d0d35d5f..0f9e75b92 100644 --- a/src/InnerLayer/SLAYER/LayerInputs.jl +++ b/src/InnerLayer/SLAYER/LayerInputs.jl @@ -179,15 +179,13 @@ is recovered after dividing by `r_s`. The reference conversion is deliberately linear (Frobenius exponent ½, the paper's H = 0 ordering): the critical-Δ is a layer-side quantity like the slab -Δ̂(Q) and W_d, whereas the outer Δ' is converted with `K^(2α)` because it -genuinely carries the Mercier exponent α. Converting the offset with `c^(2α)` -instead would multiply it by `c^(2α−1)`, 5–13 % at the DIII-D-like 3/2–4/1 -surfaces, without support from the derivation. +Δ̂(Q) and W_d, whereas the outer Δ' carries the Mercier exponent `α_M = √(−D_I)` +and is converted with `K^(2α_M)`. """ 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")) + v1 > 0 || throw(ArgumentError("toroidal_dgeo: dV/dψ_N must be positive, got $v1")) + chi1 != 0 || throw(ArgumentError("toroidal_dgeo: chi1 must be non-zero")) alpha = 2π * n * v1 / chi1 psit1 = q * chi1 / v1 lambda = psit1^2 * (-q1 / (q^2 * v1)) @@ -231,9 +229,9 @@ profiles, without an intermediate file round-trip. (`:lar`, `:rfitzp`, `:toroidal`). When `nothing` (default), Julia derives it per-surface from the equilibrium as `dr_val_k = D_R(ψ_k) = E_k + F_k + H_k²`, - consistent with Connor-Hastie-Helander 2015 (PPCF 57 065001) Eq. 59 - which uses `(−D_R)` in the χ_‖-matching critical-Δ. Pass a scalar / - vector / callable to override. + consistent with Connor, Ham, Hastie & Liu 2015 (PPCF 57 065001) Eq. 59, + which uses `(−D_R)` in the χ_‖-matching critical-Δ. Pass a scalar or a + callable of `psi` to override. **NOTE**: the χ_‖-matching critical-Δ requires the resistive interchange index `D_R = E + F + H²` (Glasser-Greene-Johnson 1975), @@ -246,8 +244,10 @@ profiles, without an intermediate file round-trip. 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. + one. Pass a scalar or a callable of `psi` to prescribe it; a prescribed + value must already be in the `r_s` reference, i.e. Eq. 59's own value + times `k_ref·v1/V_s` (≈ 2 at large aspect ratio). 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 @@ -355,7 +355,7 @@ function build_slayer_inputs(equil, sings, profiles::KineticProfiles; # dr_val: per-surface resistive interchange index D_R = E + F + H² # (Glasser-Greene-Johnson 1975). Used by `_solve_dc_tmp` to compute - # the χ_‖-matching critical-Δ via Connor-Hastie-Helander 2015 Eq. 59, + # the χ_‖-matching critical-Δ via Connor, Ham, Hastie & Liu 2015 Eq. 59, # which has `(−D_R)` as a multiplier. NOT the Mercier index # D_I = E + F + H − 1/4 (see this function's docstring); we use the # physically correct D_R here. @@ -382,7 +382,7 @@ function build_slayer_inputs(equil, sings, profiles::KineticProfiles; rs / da_dpsi else @warn("build_slayer_inputs: da/dψ = $da_dpsi at ψ = $psi is not usable; leaving " * - "Δ' unconverted (k_ref = 1) at this surface.", maxlog=3) + "Δ' and the toroidal critical-Δ unconverted (k_ref = 1) at this surface.", maxlog = 3) 1.0 end @@ -410,9 +410,10 @@ function build_slayer_inputs(equil, sings, profiles::KineticProfiles; end alpha_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) + @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)) diff --git a/src/Tearing/Runner/Control.jl b/src/Tearing/Runner/Control.jl index e95c69a62..06e997a73 100644 --- a/src/Tearing/Runner/Control.jl +++ b/src/Tearing/Runner/Control.jl @@ -23,7 +23,7 @@ constructor. (multi-surface determinant) - `dc_type` -- critical-Δ offset selector, one of `:none`, `:lar`, `:rfitzp`, `:toroidal` (χ_‖-matching critical-Δ formulas, - Connor-Hastie-Helander 2015) + Connor, Ham, Hastie & Liu 2015) - `msing_max` -- number of surfaces to include in the coupled determinant (default 3; capped at `length(sings)` at runtime) diff --git a/src/Tearing/Runner/HDF5Output.jl b/src/Tearing/Runner/HDF5Output.jl index de8f4906e..224fb1c64 100644 --- a/src/Tearing/Runner/HDF5Output.jl +++ b/src/Tearing/Runner/HDF5Output.jl @@ -93,7 +93,7 @@ const TEARING_H5_ANNOTATIONS = [ "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_offset" => (; long_name="critical-Δ offset from χ_∥/χ_⊥ matching (Connor, Ham, Hastie & Liu 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/alpha_mercier" => (; long_name="Mercier Frobenius exponent α = √(−D_I) at each surface (Glasser-Greene-Johnson 1975 Eq. 48)", dims=("surface",)), @@ -106,8 +106,14 @@ 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, ψ_N-referenced (GGJ path; identical to SingularSurfaces/Delta_prime_matrix)", dims=("surface_row", "surface_col")), - "PerSurface/Delta_prime_matrix_rs" => (; long_name="full complex Δ' matrix as used in the slab-layer matching, converted to the r_s reference length (K^(2α) on the diagonal; the ψ_N-referenced BVP matrix is SingularSurfaces/Delta_prime_matrix)", dims=("surface_row", "surface_col")), + "PerSurface/Delta_prime_matrix" => (; + long_name="full complex Δ' matrix coupling the rational surfaces, ψ_N-referenced (GGJ path; identical to SingularSurfaces/Delta_prime_matrix)", + dims=("surface_row", "surface_col") + ), + "PerSurface/Delta_prime_matrix_rs" => (; + long_name="full complex Δ' matrix as used in the slab-layer matching, converted to the r_s reference length (K^(2α) on the diagonal; 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",)), From e4e0ef2ca02e45786dca0af99af96c0579bc64ed Mon Sep 17 00:00:00 2001 From: d-burg Date: Wed, 23 Sep 2026 16:27:48 -0400 Subject: [PATCH 05/16] =?UTF-8?q?InnerLayer.SLAYER=20-=20TEST=20-=20Pin=20?= =?UTF-8?q?the=20toroidal=20critical-=CE=94=20prefactor=20against=20rfitzp?= =?UTF-8?q?=20and=20the=20stable-D=5FR=20sign?= MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit dc_tmp(:toroidal)/dc_tmp(:rfitzp) must equal dgeo/√(n|s|r_s/R₀) exactly (up to the W_d iteration's convergence), which checks the shared χ-matching prefactor end to end; a stable D_R < 0 surface must get a positive offset. The direct-evaluation and label checks are labelled as the wiring checks they are. Co-Authored-By: Claude Opus 5.5 --- test/runtests_slayer_inputs.jl | 39 +++++++++++++++++++++------------- 1 file changed, 24 insertions(+), 15 deletions(-) diff --git a/test/runtests_slayer_inputs.jl b/test/runtests_slayer_inputs.jl index 8b431884a..2cb170df2 100644 --- a/test/runtests_slayer_inputs.jl +++ b/test/runtests_slayer_inputs.jl @@ -25,7 +25,7 @@ omega_i=fill(5.0e3, length(psi_pts))) # Helper to build a minimal SingType without touching unused fields - _mk_sing(; psi, q, q1, m, n, delta_prime=-10.0+0im) = SingType( + _mk_sing(; psi, q, q1, m, n, delta_prime=-10.0 + 0im) = SingType(; psifac=psi, rho=sqrt(psi), m=[m], n=[n], q=q, q1=q1, delta_prime=ComplexF64[delta_prime], delta_prime_col=zeros(ComplexF64, 0, 0), @@ -72,8 +72,8 @@ 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)] + 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)] # dr_val=0.0 bypasses the build_slayer_inputs requirement that sing.restype be # pre-populated by ForceFreeStates.resist_eval_all! — the test sings here are # minimal stubs without restype, so we supply dr_val explicitly. @@ -126,7 +126,7 @@ # survived. It shows up on a deck whose normalization differs from its field -- the # DIII-D-like EFIT deck has b0exp = 1.0 against a physical ~1.95 T. What is pinned here is # therefore the resolution rule itself, which is deck-independent. - sings = [_mk_sing(psi=0.3, q=2.0, q1=1.5, m=2, n=1)] + sings = [_mk_sing(; psi=0.3, q=2.0, q1=1.5, m=2, n=1)] bt_phys = Float64(equil.profiles.F_spline(0.3)) / (2π * equil.ro) sl_default = build_slayer_inputs(equil, sings, profiles; dr_val=0.0, compute_omega_star=false) @@ -143,14 +143,14 @@ end @testset "build_slayer_inputs: chi_perp/chi_tor as scalars and callables" begin - sings = [_mk_sing(psi=0.5, q=2.4, q1=1.2, m=2, n=1)] + sings = [_mk_sing(; psi=0.5, q=2.4, q1=1.2, m=2, n=1)] # Scalar (dr_val=0.0 bypasses the sing.restype requirement; see comment above) sl_s = build_slayer_inputs(equil, sings, profiles; bt=2.0, chi_perp=2.0, chi_tor=1.5, dr_val=0.0) # Callable with matching value - chi_p(psi) = 2.0 + 0.0*psi - chi_t(psi) = 1.5 + 0.0*psi + chi_p(psi) = 2.0 + 0.0 * psi + chi_t(psi) = 1.5 + 0.0 * psi sl_c = build_slayer_inputs(equil, sings, profiles; bt=2.0, chi_perp=chi_p, chi_tor=chi_t, dr_val=0.0) @test sl_s[1].P_perp ≈ sl_c[1].P_perp @@ -167,7 +167,7 @@ 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)] + 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) @@ -186,7 +186,7 @@ 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)] + sings = [_mk_sing(; psi=0.5, q=2.4, q1=1.2, m=2, n=1)] # dc_type=:none and dr_val=0.0 → dc_tmp = 0 regardless of dr_val sl_none = build_slayer_inputs(equil, sings, profiles; @@ -207,7 +207,7 @@ @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) + 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; @@ -218,8 +218,18 @@ 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 + # An interchange-stable surface (D_R < 0) gets a positive, stabilizing critical-Δ. + p_stable = build_slayer_inputs(equil, [sing], profiles; + bt=2.0, dc_type=:toroidal, dr_val=-0.01)[1] + @test p_stable.dc_tmp > 0 + + # Same χ matching and D_R prefactor as :rfitzp, which carries √(n|s|r_s/R₀) in place of dgeo; + # exact up to the W_d iteration's 1e-10 convergence. + p_rf0 = build_slayer_inputs(equil, [sing], profiles; + bt=2.0, dc_type=:rfitzp, dr_val=0.01)[1] + @test p.dc_tmp / p_rf0.dc_tmp ≈ p.dgeo_val / sqrt(p.n * abs(p.sval_r) * p.rs / p.R0) rtol = 1e-8 - # Matches a direct evaluation of the r_s-referenced Eq. 59 factor. + # Wiring check, not physics: the derivation reaches toroidal_dgeo with the surface's inputs. rg = sing.restype rs = surface_minor_radius(equil, psi_s) k_ref = rs / surface_da_dpsi(equil, psi_s) @@ -238,8 +248,7 @@ 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. + # Wiring check, not physics: the radial label reaches the factor only through k_ref. 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] @@ -259,8 +268,8 @@ # and n = 2 must double Q/tauk = -ω_*, while the ratio iota_e must not move. # Asserted through the returned parameters so the check survives a refactor of # where the factor is applied. - s1 = [_mk_sing(psi=0.3, q=2.0, q1=1.5, m=2, n=1)] - s2 = [_mk_sing(psi=0.3, q=2.0, q1=1.5, m=4, n=2)] + s1 = [_mk_sing(; psi=0.3, q=2.0, q1=1.5, m=2, n=1)] + s2 = [_mk_sing(; psi=0.3, q=2.0, q1=1.5, m=4, n=2)] sl1 = build_slayer_inputs(equil, s1, profiles; bt=2.0, dr_val=0.0) sl2 = build_slayer_inputs(equil, s2, profiles; bt=2.0, dr_val=0.0) From c0dba9f7f655b013cecf695678a148fa6d3bb5c0 Mon Sep 17 00:00:00 2001 From: d-burg Date: Wed, 23 Sep 2026 16:27:48 -0400 Subject: [PATCH 06/16] =?UTF-8?q?Regression=20-=20TEST=20-=20Track=20D=5Fg?= =?UTF-8?q?eo=20and=20add=20a=20toroidal=20critical-=CE=94=20harness=20cas?= =?UTF-8?q?e?= MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Co-Authored-By: Claude Opus 5.5 --- regression-harness/cases/diiid_slayer_n1.toml | 10 +++ .../cases/diiid_slayer_n1_toroidal.toml | 71 +++++++++++++++++++ 2 files changed, 81 insertions(+) create mode 100644 regression-harness/cases/diiid_slayer_n1_toroidal.toml diff --git a/regression-harness/cases/diiid_slayer_n1.toml b/regression-harness/cases/diiid_slayer_n1.toml index 90edf2e09..c53315ba1 100644 --- a/regression-harness/cases/diiid_slayer_n1.toml +++ b/regression-harness/cases/diiid_slayer_n1.toml @@ -109,6 +109,16 @@ label = "SLAYER Q_i" noise_threshold = 1e-12 order = 21 +# Toroidal critical-Δ geometric factor, derived on every surface that carries a ResistGeometry +# whatever dc_type is selected, so it is tracked here at the deck's dc_type = "none". +[quantities.slayer_D_geo] +h5path = "Tearing/PerSurface/D_geo" +type = "real_vector" +extract = "all_real" +label = "SLAYER toroidal critical-Δ factor D_geo" +noise_threshold = 1e-10 +order = 22 + # Tearing eigenvalue (coupled mode → length 1). The headline deliverable: # a real, nonzero growth rate on a realistic equilibrium. Root-extraction # is sensitive to the AMR cell topology and ODE solver, so the Q/ω/γ diff --git a/regression-harness/cases/diiid_slayer_n1_toroidal.toml b/regression-harness/cases/diiid_slayer_n1_toroidal.toml new file mode 100644 index 000000000..6d0c220d4 --- /dev/null +++ b/regression-harness/cases/diiid_slayer_n1_toroidal.toml @@ -0,0 +1,71 @@ +# Regression case: DIII-D-like H-mode n=1 SLAYER with the toroidal critical-Δ offset. +# Same deck as diiid_slayer_n1 with dc_type = "toroidal", so the geometric factor D_geo derived +# from each surface's ResistGeometry feeds the χ_‖-matching offset D_c and, through it, the roots. +# Each [quantities.*] block names an HDF5 path in the run output, how to extract it, +# and the noise floor below which a difference is treated as zero. +[case] +name = "diiid_slayer_n1_toroidal" +description = "DIII-D-like H-mode, n=1, SLAYER with the toroidal critical-Δ offset (dc_type = toroidal)" +example_dir = "examples/DIIID-like_SLAYER_example" + +[overrides] +"SLAYER.dc_type" = "toroidal" + +# Critical-Δ inputs and offset: the resistive interchange index, the geometric factor, and D_c. +[quantities.slayer_dr_val] +h5path = "Tearing/PerSurface/D_R" +type = "real_vector" +extract = "all_real" +label = "SLAYER resistive interchange D_R" +noise_threshold = 1e-10 +order = 10 + +[quantities.slayer_D_geo] +h5path = "Tearing/PerSurface/D_geo" +type = "real_vector" +extract = "all_real" +label = "SLAYER toroidal critical-Δ factor D_geo" +noise_threshold = 1e-10 +order = 11 + +[quantities.slayer_D_c_offset] +h5path = "Tearing/PerSurface/D_c_offset" +type = "real_vector" +extract = "all_real" +label = "SLAYER critical-Δ offset D_c" +noise_threshold = 1e-10 +order = 12 + +# Roots for the inner three rational surfaces, with the same floors as diiid_slayer_n1: the +# outermost surfaces' contour search is unreliable, so they are not tracked. +[quantities.slayer_Q] +h5path = "Tearing/Roots/Q_root" +type = "complex_vector" +extract = "first_3_complex" +label = "SLAYER Q_root [2/1,3/1,4/1]" +noise_threshold = 1e-4 +order = 30 + +[quantities.slayer_gamma_Hz] +h5path = "Tearing/Roots/gamma" +type = "real_vector" +extract = "first_3" +label = "SLAYER γ_Hz [2/1,3/1,4/1]" +noise_threshold = 1e-1 +order = 33 + +[quantities.slayer_no_root] +h5path = "Tearing/Roots/no_root" +type = "real_vector" +extract = "first_3" +label = "SLAYER no_root flags [2/1,3/1,4/1]" +noise_threshold = 0 +order = 34 + +[quantities.runtime] +h5path = "" +type = "runtime" +extract = "value" +label = "Runtime (s)" +noise_threshold = 0.0 +order = 999 From 13b709efbfaf934d8978679b6332a1ce816ee521 Mon Sep 17 00:00:00 2001 From: d-burg Date: Thu, 24 Sep 2026 15:53:12 -0400 Subject: [PATCH 07/16] =?UTF-8?q?InnerLayer.SLAYER=20-=20TEST=20-=20Bound?= =?UTF-8?q?=20the=20toroidal=20critical-=CE=94=20large-aspect-ratio=20resi?= =?UTF-8?q?dual=20by=20local=20=CE=B5=20and=20require=20it=20to=20halve?= MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit The old rtol = 1e-2 sat inside the measured spread of a genuine O(ε) toroidal correction and passed only because the two sampled surfaces fell under it. The residual is now bounded by the local r_s/R₀ (measured coefficient 0.15–0.18 at ε = 0.05) and must halve from ε = 0.05 to 0.025 within 10 % (measured 0.491–0.495), which a wrong O(1) factor cannot satisfy. Co-Authored-By: Claude Opus 5.5 --- test/runtests_tj_analytic.jl | 93 +++++++++++++++++++----------------- 1 file changed, 49 insertions(+), 44 deletions(-) diff --git a/test/runtests_tj_analytic.jl b/test/runtests_tj_analytic.jl index 91d134c8d..5c39e7d0a 100644 --- a/test/runtests_tj_analytic.jl +++ b/test/runtests_tj_analytic.jl @@ -19,12 +19,12 @@ using GeneralizedPerturbedEquilibrium.Equilibrium: TJAnalyticConfig, Equilibrium @testset "TJ-analytic model" begin @testset "tj_analytic_run (inverse) — basic invariants at ε = 0.25" begin # Keep ε, mpsi, mtheta modest so the whole block runs in ~1 s. - tj = TJAnalyticConfig(lar_r0 = 1.0 / 0.25, 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) + tj = TJAnalyticConfig(; lar_r0=1.0 / 0.25, 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) # psio is a physical-scale ψ; regressions in the a→a² normalization @@ -33,12 +33,12 @@ using GeneralizedPerturbedEquilibrium.Equilibrium: TJAnalyticConfig, Equilibrium @test isfinite(pe.psio) # ν root-find pins q₂(x=1) = qa; qmax at psihigh=0.995 lands ~0.04 below. - @test pe.params.q0 ≈ 1.5 rtol = 1e-3 + @test pe.params.q0 ≈ 1.5 rtol = 1e-3 @test pe.params.qmax > 3.5 @test pe.params.qmax < 3.7 # Magnetic axis at R = R0, Z = 0 for the shifted-circle benchmark. - @test pe.ro ≈ 4.0 rtol = 1e-3 + @test pe.ro ≈ 4.0 rtol = 1e-3 @test abs(pe.zo) < 1e-8 end @@ -46,12 +46,12 @@ using GeneralizedPerturbedEquilibrium.Equilibrium: TJAnalyticConfig, Equilibrium # ε = 0.60 sits on the stable side of the ideal-external-kink pole at # ε ≈ 0.665 for this (qc, qa, pc, μ) combination. Pole-approach shape # (δW_t small, Δ' > 0 and growing) is the Option B success criterion. - tj = TJAnalyticConfig(lar_r0 = 1.0 / 0.60, 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_direct", - psilow = 0.01, psihigh = 0.995, - mpsi = 64, mtheta = 128, etol = 1e-7) + tj = TJAnalyticConfig(; lar_r0=1.0 / 0.60, 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_direct", + psilow=0.01, psihigh=0.995, + mpsi=64, mtheta=128, etol=1e-7) pe = setup_equilibrium(eq, tj) @test pe.psio > 0 @@ -59,13 +59,13 @@ using GeneralizedPerturbedEquilibrium.Equilibrium: TJAnalyticConfig, Equilibrium # Direct-GS line integration at ε=0.60 gives qmax between 3.8 and 4.0. # If the εa³·L shape terms in f_R / f_Z regress, qmax jumps above 5. - @test pe.params.q0 ≈ 1.5 rtol = 1e-2 + @test pe.params.q0 ≈ 1.5 rtol = 1e-2 @test pe.params.qmax > 3.75 @test pe.params.qmax < 4.1 # Magnetic axis at R = R0. Shafranov shift of the O-point itself is # zero by construction (H₁(0) = 0). - @test pe.ro ≈ (1.0 / 0.60) rtol = 1e-3 + @test pe.ro ≈ (1.0 / 0.60) rtol = 1e-3 @test abs(pe.zo) < 1e-4 end @@ -73,48 +73,53 @@ using GeneralizedPerturbedEquilibrium.Equilibrium: TJAnalyticConfig, Equilibrium # At the magnetic axis ψ_in should equal psio (axis convention: ψ # positive at axis, zero at LCFS); sampling well outside the LCFS should # give a negative value (the vacuum branch of psi_rz). - tj = TJAnalyticConfig(lar_r0 = 1.0 / 0.25, 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_direct", - psilow = 0.01, psihigh = 0.995, - mpsi = 64, mtheta = 128, etol = 1e-7) + tj = TJAnalyticConfig(; lar_r0=1.0 / 0.25, 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_direct", + psilow=0.01, psihigh=0.995, + mpsi=64, mtheta=128, etol=1e-7) inp = tj_analytic_run_direct(eq, tj) # ψ at the geometric axis matches psio (see DirectRunInput docstring for # the sign convention: psi_in is positive at axis, zero at LCFS). R0 = 1.0 / 0.25 - @test inp.psi_in((R0, 0.0)) ≈ inp.psio rtol = 1e-3 + @test inp.psi_in((R0, 0.0)) ≈ inp.psio rtol = 1e-3 # Well outside the LCFS → negative ψ_in (vacuum branch of the grid). 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. + @testset "toroidal critical-Δ factor → √(n s r_s/R₀) as ε → 0" begin + # Connor et al. 2015 Eq. 59 in the r_s reference must reduce to the rfitzp factor at large + # aspect ratio. The residual is a genuine O(ε) toroidal correction, so it is bounded by the + # local r_s/R₀ (measured coefficient 0.15–0.18) and must halve with ε (measured 0.49–0.50); + # a wrong O(1) factor would not shrink. 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)) + function _lar_equilibrium(eps) + tj = TJAnalyticConfig(; lar_r0=1.0 / eps, 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) + return setup_equilibrium(eq, tj) + end + function _lar_deviation(pe, psi_s, n) 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, + dgeo = toroidal_dgeo(; chi1=2π * pe.psio, 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 + return dgeo / sqrt(n * r_based_shear(rs, q, q1, da) * rs / pe.ro) - 1, rs / pe.ro + end + pe, pe_half = _lar_equilibrium(0.05), _lar_equilibrium(0.025) + for (psi_s, n) in ((0.3, 1), (0.6, 2)) + dev, eps_local = _lar_deviation(pe, psi_s, n) + @test abs(dev) <= eps_local + dev_half, _ = _lar_deviation(pe_half, psi_s, n) + @test dev_half / dev ≈ 0.5 rtol = 0.1 end end @@ -124,12 +129,12 @@ using GeneralizedPerturbedEquilibrium.Equilibrium: TJAnalyticConfig, Equilibrium 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) + 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) From bb92e49a3abee2d74578efd59c645f52cefdaac4 Mon Sep 17 00:00:00 2001 From: d-burg Date: Fri, 25 Sep 2026 11:04:24 -0400 Subject: [PATCH 08/16] ForceFreeStates - DOCS - Repair the ResistGeometry field table and cite GGJ 1975 MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit The bare | in |∇ψ|² split the Markdown table cells, truncating avg_bsq_over_dpsisq and leaving avg_dpsisq undocumented. Use ‖∇ψ_N‖, state the normalized-average convention, list the eighth integrand in the file header, and replace the Fortran routine citation in the edited resist_geometry docstring with Glasser, Greene & Johnson 1975. Co-Authored-By: Claude Opus 5.5 --- src/ForceFreeStates/Surfaces/ResistEval.jl | 19 +++++++++++-------- 1 file changed, 11 insertions(+), 8 deletions(-) diff --git a/src/ForceFreeStates/Surfaces/ResistEval.jl b/src/ForceFreeStates/Surfaces/ResistEval.jl index d653357ea..5ade7d19c 100644 --- a/src/ForceFreeStates/Surfaces/ResistEval.jl +++ b/src/ForceFreeStates/Surfaces/ResistEval.jl @@ -1,9 +1,9 @@ # ResistEval.jl # # Per-singular-surface Glasser-Greene-Johnson geometric coefficients (E, F, -# G, H, K, M) and the two flux-surface averages (⟨B²/|∇ψ|²⟩, ⟨B²⟩) that -# downstream callers need to turn geometry into τ_A / τ_R with kinetic -# profiles. +# G, H, K, M) and the three flux-surface averages (⟨B²/|∇ψ|²⟩, ⟨B²⟩, ⟨|∇ψ|²⟩) +# that downstream callers need to turn geometry into τ_A / τ_R with kinetic +# profiles and into the toroidal critical-Δ geometric factor. # # Port of Fortran RDCON `resist_eval` (geometric part only). # Unlike the Fortran, this routine produces *only* the pure-equilibrium @@ -19,6 +19,8 @@ # 4: 1 / (B² · |∇ψ|²) # 5: B² # 6: |∇ψ|² / B² +# 7: B (see below) +# 8: |∇ψ|² (toroidal critical-Δ, ⟨|∇V|²⟩ = v1²⟨|∇ψ_N|²⟩) # All weighted by `jac / v1` (jacobian / dV/dψ) before integration. # # A seventh integrand, B, is added (beyond the Fortran set) so that ⟨B⟩ is @@ -31,7 +33,8 @@ ResistGeometry Per-singular-surface Glasser-Greene-Johnson geometric coefficients and -supporting flux-surface averages. +supporting flux-surface averages. All averages ⟨·⟩ are normalized flux-surface +averages (⟨1⟩ = 1), and `‖∇ψ_N‖` is the gradient of the normalized poloidal flux. | field | meaning | |:--------------------- |:-------------------------------------------------------- | @@ -40,9 +43,9 @@ supporting flux-surface averages. | `H` | Pfirsch-Schlüter coefficient | | `K` | Glasser parameter | | `M` | Mass factor | -| `avg_bsq_over_dpsisq` | ⟨B²/ | +| `avg_bsq_over_dpsisq` | ⟨B²/‖∇ψ_N‖²⟩ — needed for τ_R | | `avg_bsq` | ⟨B²⟩ — needed for τ_R | -| `avg_dpsisq` | ⟨ | +| `avg_dpsisq` | ⟨‖∇ψ_N‖²⟩ — needed for the toroidal critical-Δ | | `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 | @@ -85,8 +88,8 @@ end """ resist_geometry(equil, psifac, q1; gamma=5/3) -> ResistGeometry -Port of Fortran RDCON `resist_eval` restricted to the -pure-equilibrium geometric coefficients. Integrates the 8 theta integrands +Glasser, Greene & Johnson 1975 (Phys. Fluids 18, 875) resistive-layer coefficients, +restricted to the pure-equilibrium geometric part. Integrates the 8 theta integrands at the given flux surface and combines the 6 GGJ ones into E, F, G, H, K, M via the standard GGJ formulas. From 127e3c619eb7268ee7c28d58a5753f0d45f6d1d9 Mon Sep 17 00:00:00 2001 From: d-burg Date: Fri, 25 Sep 2026 11:04:24 -0400 Subject: [PATCH 09/16] Tearing - DOCS - Scope the zero-override wording to dr_val and annotate the D_geo override An explicit dgeo_val = 0.0 only disables the offset for dc_type = :toroidal; dr_val = 0.0 disables it for every dc_type. The D_geo dataset holds the prescribed value when dgeo_val is overridden. Co-Authored-By: Claude Opus 5.5 --- src/Tearing/Runner/Control.jl | 7 ++++--- src/Tearing/Runner/HDF5Output.jl | 3 ++- 2 files changed, 6 insertions(+), 4 deletions(-) diff --git a/src/Tearing/Runner/Control.jl b/src/Tearing/Runner/Control.jl index 06e997a73..e440806c2 100644 --- a/src/Tearing/Runner/Control.jl +++ b/src/Tearing/Runner/Control.jl @@ -42,9 +42,10 @@ constructor. 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 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) + (consumed only by `dc_type=:toroidal`). Supply a scalar only to override + the auto-derivation. An explicit `dr_val = 0.0` disables the critical-Δ + offset (Δ_crit ≡ 0) for every `dc_type`; `dgeo_val = 0.0` does so only for + `:toroidal` - `theta_sample` -- poloidal angle at which to sample minor radius (default 0.0, outboard midplane) - `resistivity_model` -- η closure setting τ_R = μ₀r_s²/η: `:sauter` diff --git a/src/Tearing/Runner/HDF5Output.jl b/src/Tearing/Runner/HDF5Output.jl index 224fb1c64..6942a0564 100644 --- a/src/Tearing/Runner/HDF5Output.jl +++ b/src/Tearing/Runner/HDF5Output.jl @@ -90,7 +90,8 @@ 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 et al. 2015 Eq. 59 toroidal critical-Δ geometric factor in the r_s reference (0 when no ResistGeometry)", dims=("surface",)), + "PerSurface/D_geo" => + (; long_name="Connor et al. 2015 Eq. 59 toroidal critical-Δ geometric factor in the r_s reference or the dgeo_val override (0 without 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, Ham, Hastie & Liu 2015 Eq. 59)", dims=("surface",)), From d0e8a0261dd94e1102e2bed88eedc01b1644ff99 Mon Sep 17 00:00:00 2001 From: d-burg Date: Fri, 25 Sep 2026 11:04:24 -0400 Subject: [PATCH 10/16] Regression - MINOR - Track D_geo once, in diiid_slayer_n1 D_geo is derived independently of dc_type, so the toroidal variant of the same deck tracked an identical vector twice. Co-Authored-By: Claude Opus 5.5 --- .../cases/diiid_slayer_n1_toroidal.toml | 11 ++--------- 1 file changed, 2 insertions(+), 9 deletions(-) diff --git a/regression-harness/cases/diiid_slayer_n1_toroidal.toml b/regression-harness/cases/diiid_slayer_n1_toroidal.toml index 6d0c220d4..7140b34e0 100644 --- a/regression-harness/cases/diiid_slayer_n1_toroidal.toml +++ b/regression-harness/cases/diiid_slayer_n1_toroidal.toml @@ -11,7 +11,8 @@ example_dir = "examples/DIIID-like_SLAYER_example" [overrides] "SLAYER.dc_type" = "toroidal" -# Critical-Δ inputs and offset: the resistive interchange index, the geometric factor, and D_c. +# Critical-Δ input and offset: the resistive interchange index and D_c. The geometric factor +# D_geo does not depend on dc_type and is tracked once, in diiid_slayer_n1. [quantities.slayer_dr_val] h5path = "Tearing/PerSurface/D_R" type = "real_vector" @@ -20,14 +21,6 @@ label = "SLAYER resistive interchange D_R" noise_threshold = 1e-10 order = 10 -[quantities.slayer_D_geo] -h5path = "Tearing/PerSurface/D_geo" -type = "real_vector" -extract = "all_real" -label = "SLAYER toroidal critical-Δ factor D_geo" -noise_threshold = 1e-10 -order = 11 - [quantities.slayer_D_c_offset] h5path = "Tearing/PerSurface/D_c_offset" type = "real_vector" From 2762880651b2b33b90182600a55d1fb6f5dff937 Mon Sep 17 00:00:00 2001 From: d-burg Date: Sat, 26 Sep 2026 15:07:34 -0400 Subject: [PATCH 11/16] =?UTF-8?q?InnerLayer.SLAYER=20-=20FEATURE!=20-=20Us?= =?UTF-8?q?e=20toroidal=20field-line=20geometry=20in=20the=20toroidal=20cr?= =?UTF-8?q?itical-=CE=94=20parallel-conduction=20closure?= MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit dc_type = "toroidal" took Connor, Ham, Hastie & Liu 2015 (PPCF 57 065001) Eq. 59's toroidal geometry only in its final prefactor; the self-consistent W_d loop and the free-streaming χ∥ inside it still used the cylindrical field-line pitch n·s/R₀ (Fitzpatrick 1995). Both now come from the same Eq. 20 transport balance, χ⊥⟨|∇V|²⟩∂²_x = χ∥(α²Λ²/⟨B²⟩)x²: W_d = √8·(χ⊥/χ∥)^¼ / D_geo χ∥,fs = 2·v_te / (√π·K∥·W_d), K∥ = |αΛ|·r_s(dV/dr)/√⟨B²⟩ (new toroidal_kpar) K∥ is not D_geo²/r_s outside a cylinder (they differ by r_s√⟨|∇ψ_N|²⟩/k_ref, 0.57–0.76 on the DIII-D-like deck), so it is derived separately. Both reduce to the cylindrical forms at large aspect ratio (K∥ at O(ε²), D_geo at O(ε)), and the whole toroidal critical-Δ is now invariant under the choice of radial label up to k_ref. :lar and :rfitzp are algebraically unchanged. Co-Authored-By: Claude Opus 5.5 --- src/InnerLayer/InnerLayer.jl | 4 +- src/InnerLayer/SLAYER/LayerInputs.jl | 55 +++++++++++++--- src/InnerLayer/SLAYER/LayerParameters.jl | 83 ++++++++++++++---------- src/InnerLayer/SLAYER/SLAYER.jl | 2 +- src/Tearing/Runner/Control.jl | 3 +- test/runtests_slayer_inputs.jl | 30 ++++++--- test/runtests_slayer_params.jl | 37 ++++++++++- test/runtests_tj_analytic.jl | 34 +++++++--- 8 files changed, 181 insertions(+), 67 deletions(-) diff --git a/src/InnerLayer/InnerLayer.jl b/src/InnerLayer/InnerLayer.jl index 78d2df4c4..01ad7df7e 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, radial_label, build_slayer_inputs, toroidal_dgeo +import .SLAYER: surface_minor_radius, surface_da_dpsi, radial_label, build_slayer_inputs, toroidal_dgeo, toroidal_kpar 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, radial_label, build_slayer_inputs, toroidal_dgeo +export surface_minor_radius, surface_da_dpsi, radial_label, build_slayer_inputs, toroidal_dgeo, toroidal_kpar end # module InnerLayer diff --git a/src/InnerLayer/SLAYER/LayerInputs.jl b/src/InnerLayer/SLAYER/LayerInputs.jl index 0f9e75b92..6495c0fca 100644 --- a/src/InnerLayer/SLAYER/LayerInputs.jl +++ b/src/InnerLayer/SLAYER/LayerInputs.jl @@ -184,13 +184,40 @@ and is converted with `K^(2α_M)`. """ 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ψ_N must be positive, got $v1")) - chi1 != 0 || throw(ArgumentError("toroidal_dgeo: chi1 must be non-zero")) + alpha, lambda = _connor_alpha_lambda("toroidal_dgeo", chi1, v1, q, q1, n) + grad_v_sq = v1^2 * avg_dpsisq + return k_ref * v1 * (alpha^2 * lambda^2 / (avg_bsq * grad_v_sq))^0.25 +end + +""" + toroidal_kpar(; chi1, v1, q, q1, n, avg_bsq, k_ref) -> Float64 + +Parallel-wavenumber gradient `K∥` [1/m] of the toroidal critical-Δ closure: an island of normalized +width `W` (in the `x̂ = (r−r_s)/r_s` reference) sees `k∥ = K∥·W`. From Connor, Ham, Hastie & Liu 2015 +(PPCF 57 065001) Eq. 20, whose transport balance reads `χ⊥⟨|∇V|²⟩ ∂²_x = χ∥ (α²Λ²/⟨B²⟩) x²` in +`x = V − V_s`, the parallel wavenumber is `k∥ = |αΛ|·x/√⟨B²⟩`; with `x = W·r_s·(dV/dr) = W·k_ref·v1`, + +``` +K∥ = |α·Λ|·k_ref·v1 / √⟨B²⟩, +``` + +with `α`, `Λ` as in [`toroidal_dgeo`](@ref). At large aspect ratio it reduces to Fitzpatrick's +(1995, Phys. Plasmas 2 825, Eq. 132) cylindrical `n·s/R₀`. It is not a function of `D_geo` alone: +`K∥ = (D_geo²/r_s)·r_s√⟨|∇ψ_N|²⟩/k_ref`, and the last factor is 1 only in a cylinder. +""" +function toroidal_kpar(; chi1::Real, v1::Real, q::Real, q1::Real, n::Integer, avg_bsq::Real, k_ref::Real) + alpha, lambda = _connor_alpha_lambda("toroidal_kpar", chi1, v1, q, q1, n) + return abs(alpha * lambda) * k_ref * v1 / sqrt(avg_bsq) +end + +# Connor et al. 2015 α = 2πn/χ' and Λ = ψ_t'²·(ι/2π)' (' = d/dV) from ψ_N-grid quantities. +function _connor_alpha_lambda(caller, chi1, v1, q, q1, n) + v1 > 0 || throw(ArgumentError("$caller: dV/dψ_N must be positive, got $v1")) + chi1 != 0 || throw(ArgumentError("$caller: chi1 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 + return alpha, lambda end """ @@ -247,7 +274,10 @@ profiles, without an intermediate file round-trip. one. Pass a scalar or a callable of `psi` to prescribe it; a prescribed value must already be in the `r_s` reference, i.e. Eq. 59's own value times `k_ref·v1/V_s` (≈ 2 at large aspect ratio). Only `dc_type=:toroidal` - consumes it. + consumes it. The same dc_type's χ∥ closure also uses the toroidal + parallel-wavenumber gradient [`toroidal_kpar`](@ref), which is always derived + from the surface's `ResistGeometry`, also when `dgeo_val` is prescribed; + without one it falls back to the cylindrical `n·|s|/R₀` with a warning. - `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 @@ -386,9 +416,17 @@ function build_slayer_inputs(equil, sings, profiles::KineticProfiles; 1.0 end - # 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 and kpar_val: Connor et al. 2015 Eq. 59 geometric factor in the r_s reference and + # the Eq. 20 parallel-wavenumber gradient (see `toroidal_dgeo`, `toroidal_kpar`), derived + # whenever the surface carries a ResistGeometry; only dc_type=:toroidal consumes them. + kpar_val_k = if rg === nothing + dc_type === :toroidal && @warn( + "build_slayer_inputs: no ResistGeometry, so the toroidal critical-Δ χ∥ closure uses " * + "the cylindrical parallel wavenumber n·|s|/R₀.", maxlog = 1) + nothing + else + toroidal_kpar(; chi1=chi1, v1=rg.v1_local, q=q, q1=q1, n=n_res, avg_bsq=rg.avg_bsq, k_ref=k_ref_k) + end dgeo_val_k = if dgeo_val === nothing if rg !== nothing toroidal_dgeo(; chi1=chi1, v1=rg.v1_local, q=q, q1=q1, n=n_res, @@ -429,6 +467,7 @@ function build_slayer_inputs(equil, sings, profiles::KineticProfiles; m=m_res, n=n_res, dr_val=dr_val_k, dgeo_val=dgeo_val_k, + kpar_val=kpar_val_k, dc_type=dc_type, ising=k, resistivity_model=resistivity_model, f_trap=f_trap_kw, diff --git a/src/InnerLayer/SLAYER/LayerParameters.jl b/src/InnerLayer/SLAYER/LayerParameters.jl index 705ff50b3..2cce6db28 100644 --- a/src/InnerLayer/SLAYER/LayerParameters.jl +++ b/src/InnerLayer/SLAYER/LayerParameters.jl @@ -20,33 +20,34 @@ Fitzpatrick two-fluid drift-MHD SLAYER inner-layer model (Fitzpatrick de-normalization. The parametrization uses `P_perp`, `P_tor`, and `D_norm` (not the older `pr`/`pe`/`ds` set). -| field | meaning | -|:---------- |:----------------------------------------------------------------- | -| `ising` | Singular-surface index (traceability only) | -| `m`, `n` | Poloidal / toroidal mode numbers at this surface | -| `tau` | T_i / T_e | -| `lu` | Lundquist number S = τ_R / τ_H | -| `c_beta` | Compressibility √(β_local / (1 + β_local)) | -| `D_norm` | (d_β/r_s) · S^(1/3) · √ι_e (Fitzpatrick normalized scale) | -| `P_perp` | Perpendicular Prandtl number τ_R / τ_⊥ | -| `P_tor` | Toroidal-direction Prandtl number τ_R / τ_‖tor | -| `Q_e` | Normalized electron diamagnetic: −tauk · ω_*e | -| `Q_i` | Normalized ion diamagnetic: −tauk · ω_*i | -| `iota_e` | Q_e / (Q_e − Q_i) | -| `tauk` | Q-conversion factor S^(1/3) · τ_H [s] — multiplies ω to get Q | -| `tau_r` | Resistive diffusion time [s] | -| `delta_n` | Δ-normalization factor S^(1/3) / r_s [m⁻¹] | -| `rs` | Minor radius at this surface [m] | -| `R0` | Major radius [m] | -| `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 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) | +| field | meaning | +|:--------------- |:------------------------------------------------------------------------------------------------------------------------------------------------------- | +| `ising` | Singular-surface index (traceability only) | +| `m`, `n` | Poloidal / toroidal mode numbers at this surface | +| `tau` | T_i / T_e | +| `lu` | Lundquist number S = τ_R / τ_H | +| `c_beta` | Compressibility √(β_local / (1 + β_local)) | +| `D_norm` | (d_β/r_s) · S^(1/3) · √ι_e (Fitzpatrick normalized scale) | +| `P_perp` | Perpendicular Prandtl number τ_R / τ_⊥ | +| `P_tor` | Toroidal-direction Prandtl number τ_R / τ_‖tor | +| `Q_e` | Normalized electron diamagnetic: −tauk · ω_*e | +| `Q_i` | Normalized ion diamagnetic: −tauk · ω_*i | +| `iota_e` | Q_e / (Q_e − Q_i) | +| `tauk` | Q-conversion factor S^(1/3) · τ_H [s] — multiplies ω to get Q | +| `tau_r` | Resistive diffusion time [s] | +| `delta_n` | Δ-normalization factor S^(1/3) / r_s [m⁻¹] | +| `rs` | Minor radius at this surface [m] | +| `R0` | Major radius [m] | +| `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 et al. 2015 Eq. 59 toroidal critical-Δ geometric factor in the r_s reference (see `toroidal_dgeo`) | +| `kpar_val` | Parallel-wavenumber gradient K∥ [1/m] of the `:toroidal` χ∥ closure, k∥ = K∥·W_d (see `toroidal_kpar`); the cylindrical n·abs(s)/R₀ when not derived | +| `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 (for `:toroidal`, W_d and the free-streaming χ∥ use `dgeo_val` and `kpar_val`) | +| `dc_type` | Selector for `dc_tmp` formula | +| `k_ref` | Reference-length ratio K = r_s · (dψ_N/dr) at this surface (1 = no Δ' conversion) | | `alpha_mercier` | Mercier Frobenius exponent α = √(−D_I) (Glasser-Greene-Johnson 1975 Eq. 48) governing the Δ' reference-length conversion (1/2 = slab/cylindrical value) | `k_ref` and `alpha_mercier` feed the ψ_N → r_s reference-length conversion of @@ -88,6 +89,7 @@ Base.@kwdef struct SLAYERParameters <: InnerLayerParameters sval_r::Float64 dr_val::Float64 = 0.0 dgeo_val::Float64 = 0.0 + kpar_val::Float64 = 0.0 eta::Float64 d_beta::Float64 @@ -134,23 +136,32 @@ end 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, + kpar_val::Union{Real,Nothing}=nothing, max_iter::Integer=100, tol::Real=1e-10) dc_type in ALLOWED_DC_TYPES || throw(ArgumentError("SLAYERParameters: unknown dc_type=$dc_type. " * "Allowed: $(ALLOWED_DC_TYPES)")) - (dc_type === :none || dr_val == 0.0) && return 0.0 + (dc_type === :none || dr_val == 0.0 || (dc_type === :toroidal && dgeo_val == 0.0)) && return 0.0 vte = sqrt(2.0 * t_e * E_CHG / M_E) chi_par_smfp = (1.581 * tau_ee * vte^2) / (1.0 + 0.2535 * zeff) + # W_d balance χ∥k∥² = χ⊥k⊥² (Fitzpatrick 1995 Sec. VII): W_d = √8·(χ⊥/χ∥)^¼/g_w, free-streaming k∥ = K∥·W_d. + # Cylindrical g_w = √(n|s|r_s/R₀), K∥ = n|s|/R₀; :toroidal uses Connor et al. 2015 Eqs. 20/59 (toroidal_dgeo/kpar). + kpar_cyl = n_tor * abs(sval_r) / R0 + g_w, kpar = if dc_type === :toroidal + (dgeo_val, kpar_val === nothing ? kpar_cyl : kpar_val) + else + (sqrt((rs / R0) * abs(sval_r) * n_tor), kpar_cyl) + end + Wd = 0.1 converged = false for _ in 1:max_iter - chi_par_lmfp = (2.0 * R0 * vte) / (sqrt(π) * n_tor * abs(sval_r) * Wd) + chi_par_lmfp = (2.0 * vte) / (sqrt(π) * kpar * Wd) chi_par = (chi_par_smfp * chi_par_lmfp) / (chi_par_smfp + chi_par_lmfp) - Wd_new = sqrt(8.0) * (chi_perp / chi_par)^0.25 * - (1.0 / sqrt((rs / R0) * abs(sval_r) * n_tor)) + Wd_new = sqrt(8.0) * (chi_perp / chi_par)^0.25 / g_w if abs(Wd_new - Wd) / max(abs(Wd), 1e-30) < tol Wd = Wd_new converged = true @@ -160,7 +171,7 @@ function _solve_dc_tmp(; dc_type::Symbol, dr_val::Real, dgeo_val::Real, end converged || error("SLAYERParameters: Wd iteration failed to converge") - chi_par_lmfp = (2.0 * R0 * vte) / (sqrt(π) * n_tor * abs(sval_r) * Wd) + chi_par_lmfp = (2.0 * vte) / (sqrt(π) * kpar * Wd) chi_par = (chi_par_smfp * chi_par_lmfp) / (chi_par_smfp + chi_par_lmfp) if dc_type === :lar @@ -181,7 +192,7 @@ end qval, sval_r, bt, rs, R0, mu_i, zeff, chi_perp, chi_tor, m, n, - dr_val=0.0, dgeo_val=0.0, + dr_val=0.0, dgeo_val=0.0, kpar_val=nothing, dc_type=:none, ising=0, resistivity_model=SauterNeoModel(), f_trap=nothing, nu_e_star=nothing, @@ -215,6 +226,8 @@ parametrization (P_perp/P_tor/D_norm; the older magnetic/electron Prandtl - `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`) + - `kpar_val` -- parallel-wavenumber gradient K∥ [1/m] of the `:toroidal` χ∥ + closure (`toroidal_kpar`); `nothing` (default) uses the cylindrical `n·|s|/R₀` - `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, @@ -268,6 +281,7 @@ function slayer_parameters(; chi_perp::Real, chi_tor::Real, m::Integer, n::Integer, dr_val::Real=0.0, dgeo_val::Real=0.0, + kpar_val::Union{Real,Nothing}=nothing, dc_type::Symbol=:none, ising::Integer=0, resistivity_model::NeoResistivityModel=SauterNeoModel(), f_trap::Union{Real,Nothing}=nothing, @@ -382,7 +396,7 @@ function slayer_parameters(; dc_tmp = _solve_dc_tmp(; dc_type=dc_type, dr_val=dr_val, dgeo_val=dgeo_val, chi_perp=chi_perp, t_e=t_e, zeff=zeff, tau_ee=tau_ee, rs=rs, R0=R0, sval_r=sval_r, - n_tor=n) + n_tor=n, kpar_val=kpar_val) return SLAYERParameters(; ising=ising, m=m, n=n, @@ -392,6 +406,7 @@ function slayer_parameters(; tauk=tauk, tau_r=tau_r, delta_n=delta_n, rs=rs, R0=R0, bt=bt, sval_r=sval_r, dr_val=dr_val, dgeo_val=dgeo_val, + kpar_val=kpar_val === nothing ? n * abs(sval_r) / R0 : kpar_val, eta=eta, d_beta=d_beta, dc_tmp=dc_tmp, dc_type=dc_type, k_ref=k_ref, alpha_mercier=alpha_mercier diff --git a/src/InnerLayer/SLAYER/SLAYER.jl b/src/InnerLayer/SLAYER/SLAYER.jl index 201e1c3c2..5586666fa 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, radial_label, build_slayer_inputs, toroidal_dgeo +export surface_minor_radius, surface_da_dpsi, radial_label, build_slayer_inputs, toroidal_dgeo, toroidal_kpar export NeoResistivityModel, SpitzerModel, SpitzerHarmModel, SauterNeoModel, RedlNeoModel end # module SLAYER diff --git a/src/Tearing/Runner/Control.jl b/src/Tearing/Runner/Control.jl index e440806c2..4285a5576 100644 --- a/src/Tearing/Runner/Control.jl +++ b/src/Tearing/Runner/Control.jl @@ -45,7 +45,8 @@ constructor. (consumed only by `dc_type=:toroidal`). Supply a scalar only to override the auto-derivation. An explicit `dr_val = 0.0` disables the critical-Δ offset (Δ_crit ≡ 0) for every `dc_type`; `dgeo_val = 0.0` does so only for - `:toroidal` + `:toroidal`. `:toroidal` also takes the χ∥ closure's parallel wavenumber + from the equilibrium (`toroidal_kpar`), which has no override - `theta_sample` -- poloidal angle at which to sample minor radius (default 0.0, outboard midplane) - `resistivity_model` -- η closure setting τ_R = μ₀r_s²/η: `:sauter` diff --git a/test/runtests_slayer_inputs.jl b/test/runtests_slayer_inputs.jl index 2cb170df2..ba071df61 100644 --- a/test/runtests_slayer_inputs.jl +++ b/test/runtests_slayer_inputs.jl @@ -223,38 +223,48 @@ bt=2.0, dc_type=:toroidal, dr_val=-0.01)[1] @test p_stable.dc_tmp > 0 - # Same χ matching and D_R prefactor as :rfitzp, which carries √(n|s|r_s/R₀) in place of dgeo; - # exact up to the W_d iteration's 1e-10 convergence. - p_rf0 = build_slayer_inputs(equil, [sing], profiles; - bt=2.0, dc_type=:rfitzp, dr_val=0.01)[1] - @test p.dc_tmp / p_rf0.dc_tmp ≈ p.dgeo_val / sqrt(p.n * abs(p.sval_r) * p.rs / p.R0) rtol = 1e-8 - - # Wiring check, not physics: the derivation reaches toroidal_dgeo with the surface's inputs. + # Wiring check, not physics: the derivation reaches toroidal_dgeo and toroidal_kpar with the surface's inputs. 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) + kpar_ref = toroidal_kpar(; chi1=2π * equil.psio, v1=rg.v1_local, q=q_s, q1=q1_s, n=1, avg_bsq=rg.avg_bsq, k_ref=k_ref) @test p.dgeo_val ≈ dgeo_ref rtol = 1e-12 + @test p.kpar_val ≈ kpar_ref rtol = 1e-12 @test p.k_ref ≈ k_ref rtol = 1e-12 - # An explicit dgeo_val still overrides the derivation. + # K∥ is not D_geo²/r_s outside a cylinder: the metric factor r_s·√⟨|∇ψ_N|²⟩/k_ref enters. + @test p.kpar_val ≈ p.dgeo_val^2 * sqrt(rg.avg_dpsisq) / k_ref rtol = 1e-10 + + # An explicit dgeo_val still overrides the derivation; K∥ stays derived. 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 + @test p_fix.kpar_val ≈ kpar_ref rtol = 1e-12 - # The factor is derived for every dc_type once a ResistGeometry is present. + # The factors are 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 + @test p_rf.kpar_val ≈ kpar_ref rtol = 1e-12 - # Wiring check, not physics: the radial label reaches the factor only through k_ref. + # The radial label reaches the geometry only through k_ref, and the whole toroidal critical-Δ, closure + # included, scales with k_ref like the converted outer Δ' (exact up to the W_d iteration's 1e-10). 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 + @test p_lab.kpar_val / p_lab.k_ref ≈ p.kpar_val / p.k_ref rtol = 1e-10 + @test p_lab.dc_tmp / p_lab.k_ref ≈ p.dc_tmp / p.k_ref rtol = 1e-8 end + + # Without a ResistGeometry a prescribed dgeo_val runs with the cylindrical K∥, and says so. + sing_bare = _mk_sing(; psi=psi_s, q=q_s, q1=q1_s, m=2, n=1) + p_bare = (@test_logs (:warn, r"cylindrical") match_mode = :any build_slayer_inputs(equil, [sing_bare], profiles; + bt=2.0, dc_type=:toroidal, dr_val=0.01, dgeo_val=0.3))[1] + @test p_bare.kpar_val ≈ p_bare.n * abs(p_bare.sval_r) / p_bare.R0 rtol = 1e-14 end @testset "build_slayer_inputs: empty sings returns empty vector" begin diff --git a/test/runtests_slayer_params.jl b/test/runtests_slayer_params.jl index 1f665fb1d..a9d925faf 100644 --- a/test/runtests_slayer_params.jl +++ b/test/runtests_slayer_params.jl @@ -148,6 +148,41 @@ @test_throws ArgumentError slayer_parameters(; _drifts(0.0, 5.0e3)...) # no electron drift: iota_e = 0 end + @testset "Test 1e: toroidal χ∥ closure algebra" begin + # The :toroidal W_d balance takes its field-line geometry from (g_w, K∥) = (dgeo_val, kpar_val): + # W_d = √8·(χ⊥/χ∥)^¼/g_w and the free-streaming χ∥ = 2v_te/(√π·K∥·W_d). + using GeneralizedPerturbedEquilibrium.InnerLayer.SLAYER: _solve_dc_tmp + common = (; dr_val=0.01, chi_perp=1.0, t_e=1000.0, zeff=1.0, rs=0.5, R0=1.7, sval_r=1.0, n_tor=1) + vte = sqrt(2.0 * common.t_e * E_CHG / M_E) + g_cyl = sqrt(common.rs / common.R0 * common.sval_r * common.n_tor) + k_cyl = common.n_tor * common.sval_r / common.R0 + dc(dc_type, g, k; tau_ee=1e-6) = _solve_dc_tmp(; common..., dc_type=dc_type, dgeo_val=g, kpar_val=k, tau_ee=tau_ee) + + # Cylindrical geometry reproduces :lar (which is per metre, hence the r_s) and :rfitzp. + @test dc(:toroidal, g_cyl, k_cyl) ≈ common.rs * dc(:lar, 0.0, nothing) rtol = 1e-12 + @test dc(:toroidal, g_cyl, k_cyl) ≈ dc(:rfitzp, 0.0, nothing) rtol = 1e-8 + # kpar_val = nothing falls back to the cylindrical K∥. + @test dc(:toroidal, g_cyl, nothing) == dc(:toroidal, g_cyl, k_cyl) + + # Collisional limit: χ∥ is the Spitzer-Härm value, so K∥ drops out. + tau_c = 1e-20 + chi_smfp = 1.581 * tau_c * vte^2 / (1.0 + 0.2535 * common.zeff) + for (g, k) in ((0.3, 0.4), (0.3, 4.0), (0.6, 0.4)) + @test dc(:toroidal, g, k; tau_ee=tau_c) ≈ 0.5 * π^1.5 * (-common.dr_val) * (chi_smfp / common.chi_perp)^0.25 * g rtol = 1e-8 + end + + # Free-streaming limit: the fixed point solves in closed form, + # χ∥^(3/4) = 2·v_te·g_w / (√(8π)·K∥·χ⊥^(1/4)), so Δ_crit ∝ g_w^(4/3)·K∥^(-1/3). + tau_f = 1e20 + for (g, k) in ((0.3, 0.4), (0.3, 4.0), (0.6, 0.4)) + chi_fs = (2.0 * vte * g / (sqrt(8π) * k * common.chi_perp^0.25))^(4 / 3) + @test dc(:toroidal, g, k; tau_ee=tau_f) ≈ 0.5 * π^1.5 * (-common.dr_val) * (chi_fs / common.chi_perp)^0.25 * g rtol = 1e-8 + end + + # dgeo_val = 0 keeps its meaning of no toroidal offset. + @test dc(:toroidal, 0.0, k_cyl) == 0.0 + end + @testset "Test 2: r-based shear conversion" begin # Direct application of r_s · (dq/dψ) / (q · da/dψ). @test r_based_shear(0.5, 2.0, 4.0, 0.5) ≈ 2.0 @@ -191,7 +226,7 @@ # Every normalized layer quantity is bit-identical: abs(-1.0) === 1.0, # so the whole downstream chain reproduces exactly. for f in (:tau, :lu, :c_beta, :D_norm, :P_perp, :P_tor, :Q_e, :Q_i, - :iota_e, :tauk, :tau_r, :delta_n, :eta, :d_beta, :dc_tmp) + :iota_e, :tauk, :tau_r, :delta_n, :eta, :d_beta, :dc_tmp) @test getfield(neg, f) == getfield(pos, f) end diff --git a/test/runtests_tj_analytic.jl b/test/runtests_tj_analytic.jl index 5c39e7d0a..274fc7248 100644 --- a/test/runtests_tj_analytic.jl +++ b/test/runtests_tj_analytic.jl @@ -97,7 +97,7 @@ using GeneralizedPerturbedEquilibrium.Equilibrium: TJAnalyticConfig, Equilibrium # local r_s/R₀ (measured coefficient 0.15–0.18) and must halve with ε (measured 0.49–0.50); # a wrong O(1) factor would not shrink. using GeneralizedPerturbedEquilibrium.ForceFreeStates: resist_geometry - using GeneralizedPerturbedEquilibrium.InnerLayer: toroidal_dgeo, r_based_shear, + using GeneralizedPerturbedEquilibrium.InnerLayer: toroidal_dgeo, toroidal_kpar, r_based_shear, surface_minor_radius, surface_da_dpsi function _lar_equilibrium(eps) tj = TJAnalyticConfig(; lar_r0=1.0 / eps, lar_a=1.0, qc=1.5, qa=3.6, pc=0.001, mu=2.0, B0=12.0, ma=64, mtau=64) @@ -112,14 +112,20 @@ using GeneralizedPerturbedEquilibrium.Equilibrium: TJAnalyticConfig, Equilibrium da = surface_da_dpsi(pe, psi_s) dgeo = toroidal_dgeo(; chi1=2π * pe.psio, v1=rg.v1_local, q=q, q1=q1, n=n, avg_bsq=rg.avg_bsq, avg_dpsisq=rg.avg_dpsisq, k_ref=rs / da) - return dgeo / sqrt(n * r_based_shear(rs, q, q1, da) * rs / pe.ro) - 1, rs / pe.ro + kpar = toroidal_kpar(; chi1=2π * pe.psio, v1=rg.v1_local, q=q, q1=q1, n=n, avg_bsq=rg.avg_bsq, k_ref=rs / da) + s = r_based_shear(rs, q, q1, da) + return dgeo / sqrt(n * s * rs / pe.ro) - 1, kpar / (n * s / pe.ro) - 1, rs / pe.ro end pe, pe_half = _lar_equilibrium(0.05), _lar_equilibrium(0.025) for (psi_s, n) in ((0.3, 1), (0.6, 2)) - dev, eps_local = _lar_deviation(pe, psi_s, n) + dev, dev_k, eps_local = _lar_deviation(pe, psi_s, n) + dev_half, dev_k_half, _ = _lar_deviation(pe_half, psi_s, n) @test abs(dev) <= eps_local - dev_half, _ = _lar_deviation(pe_half, psi_s, n) @test dev_half / dev ≈ 0.5 rtol = 0.1 + # The χ∥ closure's K∥ reduces to Fitzpatrick's cylindrical n·s/R₀ at second order (measured + # coefficient 0.37–0.47 of ε²): it has no |∇ψ| metric factor, which carries D_geo's O(ε) term. + @test abs(dev_k) <= eps_local^2 + @test dev_k_half / dev_k ≈ 0.25 rtol = 0.1 end end @@ -127,7 +133,7 @@ using GeneralizedPerturbedEquilibrium.Equilibrium: TJAnalyticConfig, Equilibrium # 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 + using GeneralizedPerturbedEquilibrium.InnerLayer: toroidal_dgeo, toroidal_kpar, 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, @@ -139,11 +145,19 @@ using GeneralizedPerturbedEquilibrium.Equilibrium: TJAnalyticConfig, Equilibrium 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)) + q = pe.profiles.q_spline(psi) + k_ref = rs / surface_da_dpsi(pe, psi) + dgeo = toroidal_dgeo(; chi1=2π * pe.psio, v1=rg.v1_local, q=q, q1=q1, n=1, + avg_bsq=rg.avg_bsq, avg_dpsisq=rg.avg_dpsisq, k_ref=k_ref) + # K∥ is per metre, so K∥·a is the dimensionless combination. + kpar_a = a * toroidal_kpar(; chi1=2π * pe.psio, v1=rg.v1_local, q=q, q1=q1, n=1, avg_bsq=rg.avg_bsq, k_ref=k_ref) + return dgeo, kpar_a + end + d_ref, k_ref_a = _dgeo(1.0, 12.0, 0.5) + for (a, B0) in ((1.0, 6.0), (2.0, 12.0)) + d, k_a = _dgeo(a, B0, 0.5) + @test d ≈ d_ref rtol = 1e-6 + @test k_a ≈ k_ref_a rtol = 1e-6 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 From ff50706610b43ec2474bc409e8bc043fe374889f Mon Sep 17 00:00:00 2001 From: d-burg Date: Sun, 4 Oct 2026 02:21:25 -0400 Subject: [PATCH 12/16] =?UTF-8?q?InnerLayer.SLAYER=20-=20BUGFIX=20-=20Erro?= =?UTF-8?q?r=20on=20an=20unusable=20reference=20length=20under=20the=20tor?= =?UTF-8?q?oidal=20critical-=CE=94=20model?= MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit With k_ref = 1 the toroidal critical Δ stayed in the ψ_N reference while Δ' was reported per r_s. Co-Authored-By: Claude Opus 5.5 --- src/InnerLayer/SLAYER/LayerInputs.jl | 27 ++++++++++++++++++--------- test/runtests_slayer_inputs.jl | 9 +++++++++ 2 files changed, 27 insertions(+), 9 deletions(-) diff --git a/src/InnerLayer/SLAYER/LayerInputs.jl b/src/InnerLayer/SLAYER/LayerInputs.jl index 6495c0fca..c0f106af7 100644 --- a/src/InnerLayer/SLAYER/LayerInputs.jl +++ b/src/InnerLayer/SLAYER/LayerInputs.jl @@ -154,6 +154,22 @@ function radial_label(equil; rs_method::Symbol=:midplane, theta::Real=0.0) return (_rs_at, _da_dpsi_at) end +# Reference-length factor K = r_s·(dψ_N/dr)|_s. When da/dψ is not a usable positive number (zero, negative +# or non-finite would corrupt the Δ' diagonal), Δ' is left unconverted with K = 1; under dc_type=:toroidal +# that is an error, because the critical Δ would be compared with Δ' in a different reference. +function _reference_length_factor(rs::Real, da_dpsi::Real, psi::Real, dc_type::Symbol) + isfinite(da_dpsi) && da_dpsi > 0.0 && return rs / da_dpsi + dc_type === :toroidal && throw( + ArgumentError( + "build_slayer_inputs: da/dψ = $da_dpsi at ψ = $psi is not usable, so the toroidal critical-Δ " * + "cannot be converted to the r_s reference at this surface." + ) + ) + @warn("build_slayer_inputs: da/dψ = $da_dpsi at ψ = $psi is not usable; leaving " * + "Δ' unconverted (k_ref = 1) at this surface.", maxlog = 3) + return 1.0 +end + """ toroidal_dgeo(; chi1, v1, q, q1, n, avg_bsq, avg_dpsisq, k_ref) -> Float64 @@ -406,15 +422,8 @@ function build_slayer_inputs(equil, sings, profiles::KineticProfiles; # Reference-length conversion inputs for the outer Δ': K = r_s·(dψ_N/dr)|_s and # α = √(−D_I) (Glasser-Greene-Johnson 1975 Eq. 48), with α clamped to 0 on Mercier-unstable - # surfaces (the factor turns complex there) and K = 1 whenever da/dψ is not a usable - # positive number (zero/non-finite/negative would corrupt the Δ' diagonal). - k_ref_k = if isfinite(da_dpsi) && da_dpsi > 0.0 - rs / da_dpsi - else - @warn("build_slayer_inputs: da/dψ = $da_dpsi at ψ = $psi is not usable; leaving " * - "Δ' and the toroidal critical-Δ unconverted (k_ref = 1) at this surface.", maxlog = 3) - 1.0 - end + # surfaces (the factor turns complex there). + k_ref_k = _reference_length_factor(rs, da_dpsi, psi, dc_type) # dgeo_val and kpar_val: Connor et al. 2015 Eq. 59 geometric factor in the r_s reference and # the Eq. 20 parallel-wavenumber gradient (see `toroidal_dgeo`, `toroidal_kpar`), derived diff --git a/test/runtests_slayer_inputs.jl b/test/runtests_slayer_inputs.jl index ba071df61..a092bd418 100644 --- a/test/runtests_slayer_inputs.jl +++ b/test/runtests_slayer_inputs.jl @@ -267,6 +267,15 @@ @test p_bare.kpar_val ≈ p_bare.n * abs(p_bare.sval_r) / p_bare.R0 rtol = 1e-14 end + @testset "reference-length factor: an unusable da/dψ is an error only under dc_type=:toroidal" begin + ref_factor = InnerLayer.SLAYER._reference_length_factor + @test ref_factor(0.5, 0.25, 0.4, :toroidal) == 2.0 + for da_dpsi in (0.0, -0.25, NaN, Inf) + @test_throws r"cannot be converted to the r_s reference" ref_factor(0.5, da_dpsi, 0.4, :toroidal) + @test (@test_logs (:warn, r"leaving Δ' unconverted") ref_factor(0.5, da_dpsi, 0.4, :rfitzp)) == 1.0 + 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} From e3274310dc96e39c7263a2d6cf778cea5b89b07f Mon Sep 17 00:00:00 2001 From: d-burg Date: Sun, 4 Oct 2026 02:21:25 -0400 Subject: [PATCH 13/16] =?UTF-8?q?Repo=20-=20DOCS=20-=20Cite=20Connor=20et?= =?UTF-8?q?=20al.=202015=20for=20the=20toroidal=20critical-=CE=94?= MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Co-Authored-By: Claude Opus 5.5 --- docs/development/references.md | 6 ++++++ 1 file changed, 6 insertions(+) diff --git a/docs/development/references.md b/docs/development/references.md index e19b92e32..abf7fd223 100644 --- a/docs/development/references.md +++ b/docs/development/references.md @@ -108,3 +108,9 @@ The KineticForces module (formerly PENTRC) implements neoclassical toroidal visc - Published: Physical Review Letters **102**, 065002 (2009) - Link: https://doi.org/10.1103/PhysRevLett.102.065002 - Describes: Trapped-particle nonambipolar transport theory underpinning the NTV calculation + +- **Connor et al. (2015)**: "The role of thermal conduction in tearing mode theory" + - Location: not in `docs/resources/`; preprint at https://arxiv.org/abs/1410.7240 + - Published: Plasma Physics and Controlled Fusion **57**, 065001 (2015) + - Link: https://doi.org/10.1088/0741-3335/57/6/065001 + - Describes: Toroidal critical Δ′ from finite parallel thermal conduction at a rational surface. Eq. 59 gives the geometric factor (`toroidal_dgeo`), Eq. 20 the parallel-wavenumber gradient (`toroidal_kpar`), and Eqs. 61 and 65 the critical width `W_d` used by `dc_type = :toroidal` in `InnerLayer.SLAYER` From ed203bd46287f7b14416aac138c3fc1654f948a3 Mon Sep 17 00:00:00 2001 From: d-burg Date: Sun, 4 Oct 2026 02:21:25 -0400 Subject: [PATCH 14/16] Regression - MINOR - Declare the tolerance class of each tracked quantity Co-Authored-By: Claude Opus 5.5 --- regression-harness/cases/diiid_slayer_n1.toml | 1 + regression-harness/cases/diiid_slayer_n1_toroidal.toml | 6 ++++++ 2 files changed, 7 insertions(+) diff --git a/regression-harness/cases/diiid_slayer_n1.toml b/regression-harness/cases/diiid_slayer_n1.toml index c53315ba1..c43964ea4 100644 --- a/regression-harness/cases/diiid_slayer_n1.toml +++ b/regression-harness/cases/diiid_slayer_n1.toml @@ -118,6 +118,7 @@ extract = "all_real" label = "SLAYER toroidal critical-Δ factor D_geo" noise_threshold = 1e-10 order = 22 +class = "physics_converged" # Tearing eigenvalue (coupled mode → length 1). The headline deliverable: # a real, nonzero growth rate on a realistic equilibrium. Root-extraction diff --git a/regression-harness/cases/diiid_slayer_n1_toroidal.toml b/regression-harness/cases/diiid_slayer_n1_toroidal.toml index 7140b34e0..36f978321 100644 --- a/regression-harness/cases/diiid_slayer_n1_toroidal.toml +++ b/regression-harness/cases/diiid_slayer_n1_toroidal.toml @@ -20,6 +20,7 @@ extract = "all_real" label = "SLAYER resistive interchange D_R" noise_threshold = 1e-10 order = 10 +class = "physics_converged" [quantities.slayer_D_c_offset] h5path = "Tearing/PerSurface/D_c_offset" @@ -28,6 +29,7 @@ extract = "all_real" label = "SLAYER critical-Δ offset D_c" noise_threshold = 1e-10 order = 12 +class = "physics_converged" # Roots for the inner three rational surfaces, with the same floors as diiid_slayer_n1: the # outermost surfaces' contour search is unreliable, so they are not tracked. @@ -38,6 +40,7 @@ extract = "first_3_complex" label = "SLAYER Q_root [2/1,3/1,4/1]" noise_threshold = 1e-4 order = 30 +class = "physics_converged" [quantities.slayer_gamma_Hz] h5path = "Tearing/Roots/gamma" @@ -46,6 +49,7 @@ extract = "first_3" label = "SLAYER γ_Hz [2/1,3/1,4/1]" noise_threshold = 1e-1 order = 33 +class = "physics_converged" [quantities.slayer_no_root] h5path = "Tearing/Roots/no_root" @@ -54,6 +58,7 @@ extract = "first_3" label = "SLAYER no_root flags [2/1,3/1,4/1]" noise_threshold = 0 order = 34 +class = "physics_converged" [quantities.runtime] h5path = "" @@ -62,3 +67,4 @@ extract = "value" label = "Runtime (s)" noise_threshold = 0.0 order = 999 +class = "diagnostic" From ade1ba93d5e752b47e5eaed1ad63eb3d21c03902 Mon Sep 17 00:00:00 2001 From: d-burg Date: Tue, 6 Oct 2026 00:15:05 -0400 Subject: [PATCH 15/16] =?UTF-8?q?InnerLayer.SLAYER=20-=20BUGFIX=20-=20Go?= =?UTF-8?q?=20back=20to=20warning=20on=20an=20unusable=20reference=20lengt?= =?UTF-8?q?h=20for=20every=20critical-=CE=94=20model?= MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Erroring under the toroidal model only was inconsistent: with k_ref = 1 the derived toroidal critical Δ stays consistent with Δ', and the mismatch against the layer response exists under every model. This restores the behaviour the branch had before ff5070661. Co-Authored-By: Claude Opus 5.5 --- src/InnerLayer/SLAYER/LayerInputs.jl | 27 +++++++++------------------ test/runtests_slayer_inputs.jl | 9 --------- 2 files changed, 9 insertions(+), 27 deletions(-) diff --git a/src/InnerLayer/SLAYER/LayerInputs.jl b/src/InnerLayer/SLAYER/LayerInputs.jl index c0f106af7..6495c0fca 100644 --- a/src/InnerLayer/SLAYER/LayerInputs.jl +++ b/src/InnerLayer/SLAYER/LayerInputs.jl @@ -154,22 +154,6 @@ function radial_label(equil; rs_method::Symbol=:midplane, theta::Real=0.0) return (_rs_at, _da_dpsi_at) end -# Reference-length factor K = r_s·(dψ_N/dr)|_s. When da/dψ is not a usable positive number (zero, negative -# or non-finite would corrupt the Δ' diagonal), Δ' is left unconverted with K = 1; under dc_type=:toroidal -# that is an error, because the critical Δ would be compared with Δ' in a different reference. -function _reference_length_factor(rs::Real, da_dpsi::Real, psi::Real, dc_type::Symbol) - isfinite(da_dpsi) && da_dpsi > 0.0 && return rs / da_dpsi - dc_type === :toroidal && throw( - ArgumentError( - "build_slayer_inputs: da/dψ = $da_dpsi at ψ = $psi is not usable, so the toroidal critical-Δ " * - "cannot be converted to the r_s reference at this surface." - ) - ) - @warn("build_slayer_inputs: da/dψ = $da_dpsi at ψ = $psi is not usable; leaving " * - "Δ' unconverted (k_ref = 1) at this surface.", maxlog = 3) - return 1.0 -end - """ toroidal_dgeo(; chi1, v1, q, q1, n, avg_bsq, avg_dpsisq, k_ref) -> Float64 @@ -422,8 +406,15 @@ function build_slayer_inputs(equil, sings, profiles::KineticProfiles; # Reference-length conversion inputs for the outer Δ': K = r_s·(dψ_N/dr)|_s and # α = √(−D_I) (Glasser-Greene-Johnson 1975 Eq. 48), with α clamped to 0 on Mercier-unstable - # surfaces (the factor turns complex there). - k_ref_k = _reference_length_factor(rs, da_dpsi, psi, dc_type) + # surfaces (the factor turns complex there) and K = 1 whenever da/dψ is not a usable + # positive number (zero/non-finite/negative would corrupt the Δ' diagonal). + k_ref_k = if isfinite(da_dpsi) && da_dpsi > 0.0 + rs / da_dpsi + else + @warn("build_slayer_inputs: da/dψ = $da_dpsi at ψ = $psi is not usable; leaving " * + "Δ' and the toroidal critical-Δ unconverted (k_ref = 1) at this surface.", maxlog = 3) + 1.0 + end # dgeo_val and kpar_val: Connor et al. 2015 Eq. 59 geometric factor in the r_s reference and # the Eq. 20 parallel-wavenumber gradient (see `toroidal_dgeo`, `toroidal_kpar`), derived diff --git a/test/runtests_slayer_inputs.jl b/test/runtests_slayer_inputs.jl index a092bd418..ba071df61 100644 --- a/test/runtests_slayer_inputs.jl +++ b/test/runtests_slayer_inputs.jl @@ -267,15 +267,6 @@ @test p_bare.kpar_val ≈ p_bare.n * abs(p_bare.sval_r) / p_bare.R0 rtol = 1e-14 end - @testset "reference-length factor: an unusable da/dψ is an error only under dc_type=:toroidal" begin - ref_factor = InnerLayer.SLAYER._reference_length_factor - @test ref_factor(0.5, 0.25, 0.4, :toroidal) == 2.0 - for da_dpsi in (0.0, -0.25, NaN, Inf) - @test_throws r"cannot be converted to the r_s reference" ref_factor(0.5, da_dpsi, 0.4, :toroidal) - @test (@test_logs (:warn, r"leaving Δ' unconverted") ref_factor(0.5, da_dpsi, 0.4, :rfitzp)) == 1.0 - 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} From e84e1afa77ad11684b65959bcfcfe5df7118b797 Mon Sep 17 00:00:00 2001 From: d-burg Date: Tue, 6 Oct 2026 00:15:05 -0400 Subject: [PATCH 16/16] Repo - DOCS - Attribute the Connor et al. 2015 equations as the code cites them Co-Authored-By: Claude Opus 5.5 --- docs/development/references.md | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/docs/development/references.md b/docs/development/references.md index abf7fd223..a815ca6f8 100644 --- a/docs/development/references.md +++ b/docs/development/references.md @@ -113,4 +113,4 @@ The KineticForces module (formerly PENTRC) implements neoclassical toroidal visc - Location: not in `docs/resources/`; preprint at https://arxiv.org/abs/1410.7240 - Published: Plasma Physics and Controlled Fusion **57**, 065001 (2015) - Link: https://doi.org/10.1088/0741-3335/57/6/065001 - - Describes: Toroidal critical Δ′ from finite parallel thermal conduction at a rational surface. Eq. 59 gives the geometric factor (`toroidal_dgeo`), Eq. 20 the parallel-wavenumber gradient (`toroidal_kpar`), and Eqs. 61 and 65 the critical width `W_d` used by `dc_type = :toroidal` in `InnerLayer.SLAYER` + - Describes: Toroidal critical Δ′ from finite parallel thermal conduction at a rational surface. Eq. 59 gives the geometric factor (`toroidal_dgeo`) and Eq. 20 the parallel-wavenumber gradient (`toroidal_kpar`) used by `dc_type = :toroidal` in `InnerLayer.SLAYER`; Eq. 61 is the large-aspect-ratio form of Eq. 59, and Eq. 65 gives the corresponding critical island width of Fitzpatrick (1995)