From 9ea6d15943efb0b6c6774d0660bf60a7c841f840 Mon Sep 17 00:00:00 2001 From: d-burg Date: Mon, 17 Aug 2026 16:42:57 -0400 Subject: [PATCH 1/3] =?UTF-8?q?TEARING=20-=20BUG=20FIX=20-=20Convert=20?= =?UTF-8?q?=CE=94'=20to=20the=20r=5Fs=20reference=20length=20before=20slab?= =?UTF-8?q?-layer=20matching?= MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit The SLAYER dispersion relation paired a ψ_N-referenced BVP Δ' with an r_s-referenced layer side (Δ(Q), S=τ_R/τ_H on r_s, and the critical-Δ, whose Ŵ_d is W/r_s). Apply the Frobenius reference-length transform Δ̂_ij = K_i^(1/2+μ_i)·Δ'_ij·K_j^(μ_j−1/2), K = r_s·(dψ_N/dr)|_s, μ = √(−D_I), at the matching point. Verified parameter-free against the TJ circular benchmarks (median residual +1% over 23 points; absolute 2/1 agreement 17.5%→2.5% (β) and 17.6%→1.1% (ε)). BVP Δ' outputs are unchanged; only the SLAYER matching (γ, Δ_eff) moves. GGJ untouched. Co-Authored-By: Claude Opus 5 --- src/InnerLayer/SLAYER/LayerInputs.jl | 18 +++++++++- src/InnerLayer/SLAYER/LayerParameters.jl | 26 ++++++++++++-- src/Tearing/Dispersion/SurfaceCoupling.jl | 12 ++++--- src/Tearing/Runner/HDF5Output.jl | 10 +++++- src/Tearing/Runner/Result.jl | 4 ++- src/Tearing/Runner/run_slayer.jl | 41 +++++++++++++++++++++++ test/runtests_slayer_runner.jl | 35 +++++++++++++++++++ 7 files changed, 137 insertions(+), 9 deletions(-) diff --git a/src/InnerLayer/SLAYER/LayerInputs.jl b/src/InnerLayer/SLAYER/LayerInputs.jl index 96177903d..8a408dd29 100644 --- a/src/InnerLayer/SLAYER/LayerInputs.jl +++ b/src/InnerLayer/SLAYER/LayerInputs.jl @@ -289,6 +289,20 @@ function build_slayer_inputs(equil, sings, profiles::KineticProfiles; _eval(dgeo_val, psi) end + # Reference-length conversion inputs for the outer Δ': K = r_s·(dψ_N/dr)|_s + # and μ = √(−D_I) with D_I = E + F + H − 1/4 (Glasser-Greene-Johnson 1975). + # A Mercier-unstable surface (D_I ≥ 0) has no real exponent; μ → 0 there, + # which makes the conversion factor K^(2μ) → 1 continuously (Δ' left raw). + k_ref_k = da_dpsi == 0.0 ? 1.0 : rs / da_dpsi + mu_k = if rg === nothing + @warn("build_slayer_inputs: sing.restype not populated; using the " * + "slab Mercier exponent μ = 1/2 for the Δ' reference-length " * + "conversion at all such surfaces.", maxlog=1) + 0.5 + else + sqrt(max(-(rg.E + rg.F + rg.H - 0.25), 0.0)) + end + out[k] = slayer_parameters(; n_e=prof.n_e, t_e=prof.T_e, t_i=prof.T_i, omega=prof.omega, omega_e=ω_e_use, omega_i=ω_i_use, @@ -304,7 +318,9 @@ function build_slayer_inputs(equil, sings, profiles::KineticProfiles; f_trap=f_trap_kw, nu_e_star=nu_e_star_kw, R_major_eff=R_major_eff, - lnLambda_form=lnLambda_form + lnLambda_form=lnLambda_form, + k_ref=k_ref_k, + mu_mercier=mu_k ) end return out diff --git a/src/InnerLayer/SLAYER/LayerParameters.jl b/src/InnerLayer/SLAYER/LayerParameters.jl index 5dbf9163f..28976e6d3 100644 --- a/src/InnerLayer/SLAYER/LayerParameters.jl +++ b/src/InnerLayer/SLAYER/LayerParameters.jl @@ -46,6 +46,15 @@ de-normalization. The parametrization uses `P_perp`, `P_tor`, and | `d_beta` | Beta-weighted ion length scale c_β · d_i [m] | | `dc_tmp` | Critical-Δ offset from chi_parallel matching | | `dc_type` | Selector for `dc_tmp` formula | +| `k_ref` | Reference-length ratio K = r_s · (dψ_N/dr) at this surface (1 = no Δ' conversion) | +| `mu_mercier` | Mercier exponent μ = √(−D_I) governing the Δ' reference-length conversion (1/2 = slab/cylindrical value) | + +`k_ref` and `mu_mercier` feed the ψ_N → r_s reference-length conversion of +the outer Δ' matrix (see `delta_prime_to_rs_reference` in the Tearing +runner): the slab layer, its `dc_tmp` critical-Δ, and the `S^(1/3)` Δ(Q) +scale are all referenced to unit `x̂ = (r−r_s)/r_s`, while the outer BVP Δ' +is referenced to unit Δψ_N. Hand-built parameters default to `k_ref = 1`, +which makes the conversion the identity. The complex normalized growth rate `Q = ω + iγ` is **not** stored here; it is passed as a separate argument to `solve_inner`. @@ -85,6 +94,10 @@ Base.@kwdef struct SLAYERParameters <: InnerLayerParameters # Critical-Δ offset dc_tmp::Float64 = 0.0 dc_type::Symbol = :none + + # Reference-length conversion inputs for the outer Δ' (ψ_N → r_s-based x̂) + k_ref::Float64 = 1.0 + mu_mercier::Float64 = 0.5 end # Allowed dc_type values for the critical-Δ offset. `:none` is the default @@ -202,6 +215,12 @@ parametrization (P_perp/P_tor/D_norm; the older magnetic/electron Prandtl - `dr_val`, `dgeo_val` -- inputs for the critical-Δ formula - `dc_type` -- one of `:none`, `:lar`, `:rfitzp`, `:toroidal` - `ising` -- singular-surface index for traceability + - `k_ref` -- reference-length ratio K = r_s·(dψ_N/dr) at the surface, + used by the Tearing runner to convert the ψ_N-referenced outer Δ' to + the r_s-referenced convention this layer works in (default `1.0`, + i.e. no conversion; `build_slayer_inputs` fills it from the equilibrium) + - `mu_mercier` -- Mercier exponent μ = √(−D_I) for the same conversion + (default `0.5`, the slab/cylindrical value at D_I = −1/4) # Resistivity kwargs @@ -252,7 +271,9 @@ function slayer_parameters(; f_trap::Union{Real,Nothing}=nothing, nu_e_star::Union{Real,Nothing}=nothing, R_major_eff::Union{Real,Nothing}=nothing, - lnLambda_form::Symbol=:nrl) + lnLambda_form::Symbol=:nrl, + k_ref::Real=1.0, + mu_mercier::Real=0.5) # Coulomb logarithm shared by the resistivity closure and τ_ee. lnLamb = coulomb_log_e(n_e, t_e; form=lnLambda_form) @@ -358,6 +379,7 @@ function slayer_parameters(; rs=rs, R0=R0, bt=bt, sval_r=sval_r, dr_val=dr_val, dgeo_val=dgeo_val, eta=eta, d_beta=d_beta, - dc_tmp=dc_tmp, dc_type=dc_type + dc_tmp=dc_tmp, dc_type=dc_type, + k_ref=k_ref, mu_mercier=mu_mercier ) end diff --git a/src/Tearing/Dispersion/SurfaceCoupling.jl b/src/Tearing/Dispersion/SurfaceCoupling.jl index 271162c1a..0d57207ae 100644 --- a/src/Tearing/Dispersion/SurfaceCoupling.jl +++ b/src/Tearing/Dispersion/SurfaceCoupling.jl @@ -53,10 +53,14 @@ end surface_coupling(model::SLAYERModel, params::SLAYERParameters, dp_diag::Number; dc::Real=0.0) -> SurfaceCoupling -SLAYER convenience constructor. `scale` is set to `params.lu^(1/3)` so that -the dimensionless Δ from `riccati_f` is mapped to outer ψ-units before -subtraction from the Δ' diagonal. `tauk` is taken from `params.tauk` for use -by `MultiSurfaceCoupling` Q rescaling. +SLAYER convenience constructor. `scale` is set to `params.lu^(1/3)`, which +maps the dimensionless inner-layer Δ from `riccati_f` to the r_s-referenced +outer convention (S^(1/3) with S = τ_R/τ_H built on r_s; the Lundquist +number, `dc`, and the layer all share the `x̂ = (r−r_s)/r_s` reference +length). `dp_diag` must already be in that same r_s reference — the Tearing +runner converts the ψ_N-referenced BVP Δ' via `delta_prime_to_rs_reference` +before building couplings. `tauk` is taken from `params.tauk` for use by +`MultiSurfaceCoupling` Q rescaling. """ function surface_coupling(model::SLAYERModel, params::SLAYERParameters, dp_diag::Number; dc::Real=0.0) diff --git a/src/Tearing/Runner/HDF5Output.jl b/src/Tearing/Runner/HDF5Output.jl index ca4287147..b0f9110f3 100644 --- a/src/Tearing/Runner/HDF5Output.jl +++ b/src/Tearing/Runner/HDF5Output.jl @@ -95,6 +95,9 @@ const TEARING_H5_ANNOTATIONS = [ "PerSurface/d_beta" => (; long_name="β-weighted ion drift scale d_β", units="m", dims=("surface",)), "PerSurface/D_c_offset" => (; long_name="critical-Δ offset from χ_∥/χ_⊥ matching (Connor-Hastie-Helander 2015 Eq. 59)", dims=("surface",)), "PerSurface/D_c_type" => (; long_name="per-surface D_c prescription label", dims=("surface",)), + "PerSurface/k_ref" => (; long_name="reference-length ratio K = r_s·(dψ_N/dr) at each surface", dims=("surface",)), + "PerSurface/mu_mercier" => (; long_name="Mercier exponent μ = √(−D_I) at each surface", dims=("surface",)), + "PerSurface/delta_prime_conversion" => (; long_name="ψ_N → r_s reference-length factor K^(2μ) applied to each Δ' diagonal", dims=("surface",)), "PerSurface/E" => (; long_name="Glasser-Greene-Johnson coefficient E per surface", dims=("surface",)), "PerSurface/F" => (; long_name="Glasser-Greene-Johnson coefficient F per surface", dims=("surface",)), "PerSurface/G" => (; long_name="Glasser-Greene-Johnson coefficient G per surface", dims=("surface",)), @@ -103,7 +106,7 @@ const TEARING_H5_ANNOTATIONS = [ "PerSurface/M" => (; long_name="Glasser-Greene-Johnson coefficient M per surface", dims=("surface",)), "PerSurface/tau_A" => (; long_name="Alfvén time τ_A per surface (GGJ layer parameters)", units="s", dims=("surface",)), "PerSurface/dVdpsi" => (; long_name="dV/dψ_N at each surface", units="m^3", dims=("surface",)), - "PerSurface/Delta_prime_matrix" => (; long_name="full complex Δ' matrix coupling the rational surfaces", dims=("surface_row", "surface_col")), + "PerSurface/Delta_prime_matrix" => (; long_name="full complex Δ' matrix coupling the rational surfaces, as used in the matching (SLAYER path: r_s-referenced via K^(2μ); the ψ_N-referenced BVP matrix is SingularSurfaces/Delta_prime_matrix)", dims=("surface_row", "surface_col")), "Roots/Q_root" => (; long_name="complex dispersion-root normalized frequency Q (NaN = no root)", dims=("surface",)), "Roots/omega" => (; long_name="mode rotation angular frequency ω = Re(Q)/τ_k of each root", units="rad/s", dims=("surface",)), @@ -184,6 +187,11 @@ function _write_per_surface!(g, params::AbstractVector{SLAYERParameters}, ps["D_geo"] = Float64[p.dgeo_val for p in params] # Store dc_type per-surface as string array ps["D_c_type"] = String[String(p.dc_type) for p in params] + # Reference-length conversion applied to Δ' (ψ_N → r_s): K, μ, and the + # diagonal factor K^(2μ) actually multiplying each Δ'_kk. + ps["k_ref"] = Float64[p.k_ref for p in params] + ps["mu_mercier"] = Float64[p.mu_mercier for p in params] + ps["delta_prime_conversion"] = Float64[p.k_ref^(2 * p.mu_mercier) for p in params] ps["Delta_prime_matrix"] = dp_matrix return nothing diff --git a/src/Tearing/Runner/Result.jl b/src/Tearing/Runner/Result.jl index 0cffc21c9..d81508989 100644 --- a/src/Tearing/Runner/Result.jl +++ b/src/Tearing/Runner/Result.jl @@ -20,7 +20,9 @@ downstream inspection and HDF5 output. factor q of each analyzed surface, aligned with `params`. Empty when the analysis was built from bare parameters (`run_slayer_from_inputs` without the surface list), in which case the HDF5 writer skips them. - - `dp_matrix` -- outer-region Δ' matrix used in the analysis + - `dp_matrix` -- outer-region Δ' matrix used in the analysis. + SLAYER path: r_s-referenced (the ψ_N BVP matrix transformed by + `delta_prime_to_rs_reference`); GGJ path: the ψ_N matrix unchanged - `Q_root` -- tearing eigenvalue(s) in normalized Q * length `nsurfaces` in `:uncoupled` mode * length `1` in `:coupled` mode (global eigenvalue normalized by diff --git a/src/Tearing/Runner/run_slayer.jl b/src/Tearing/Runner/run_slayer.jl index 2a29c63a8..650761b3c 100644 --- a/src/Tearing/Runner/run_slayer.jl +++ b/src/Tearing/Runner/run_slayer.jl @@ -151,6 +151,39 @@ function _build_surface_coupling(model::GGJModel, params::GGJParameters, return surface_coupling(model, params, dp_diag) end +# --------------------------------------------------------------------- +# Reference-length conversion of the outer Δ' for the slab layer +# --------------------------------------------------------------------- +""" + delta_prime_to_rs_reference(dp_matrix, params) -> Matrix{ComplexF64} + +Convert the outer-region Δ' matrix from its ψ_N reference length (the +STRIDE/BVP convention: Frobenius coefficients normalized per unit Δψ_N) to +the r_s-based `x̂ = (r − r_s)/r_s` reference the slab layer works in +(Fitzpatrick 2023 convention — the same reference used by `dc_tmp` and by +the `S^(1/3)` Δ(Q) scale, both built on r_s). + +Near surface `k` the tearing-parity outer solution is +`A_L·|x|^(1/2−μ) + A_S·|x|^(1/2+μ)` with `μ = √(−D_I)` (Glasser, Wang & +Park 2016, Eq. 26). Rescaling the radial variable `x_ψ = K·x̂` with +`K = r_s·(dψ_N/dr)|_s` maps the coefficients as `Â_L = A_L·K^(1/2−μ)` and +`Â_S = A_S·K^(1/2+μ)`, so the response matrix (small coefficient at surface +`i` per unit large coefficient at surface `j`) transforms as + + Δ̂_ij = K_i^(1/2+μ_i) · Δ'_ij · K_j^(μ_j−1/2) + +whose diagonal is `K^(2μ)·Δ'_kk`; at D_I = −1/4 (μ = 1/2) this reduces to +the textbook `Δ̂ = r_s·Δ'_phys`. `K` and `μ` are carried per surface in +`SLAYERParameters.k_ref` / `.mu_mercier`; hand-built parameters default to +`k_ref = 1`, making the conversion the identity. +""" +function delta_prime_to_rs_reference(dp_matrix::AbstractMatrix, + params::AbstractVector) + dl = [p.k_ref^(0.5 + p.mu_mercier) for p in params] + dr = [p.k_ref^(p.mu_mercier - 0.5) for p in params] + return Diagonal(dl) * Matrix{ComplexF64}(dp_matrix) * Diagonal(dr) +end + # --------------------------------------------------------------------- # Core analysis entry point that takes pre-built parameters. # --------------------------------------------------------------------- @@ -194,6 +227,14 @@ function run_slayer_from_inputs(params::AbstractVector{<:InnerLayerParameters}, "(build_slayer_inputs for SLAYER, build_ggj_inputs for GGJ).") ) + # Slab-layer path: convert Δ' from its ψ_N reference length to the + # r_s-based convention shared by the layer Δ(Q) and the critical-Δ (see + # `delta_prime_to_rs_reference`). GGJ is genuinely toroidal/ψ-based + # (its `rescale_delta` handles inner→outer units natively) — no conversion. + if !_is_ggj(model) + dp = delta_prime_to_rs_reference(dp, params) + end + # The coupled determinant uses the reduced m×m (tearing-only) form, which # drops the interchange channel. For GGJ that channel carries the Glasser # interchange stabilization, so coupled-GGJ results omit real physics. diff --git a/test/runtests_slayer_runner.jl b/test/runtests_slayer_runner.jl index 0eacc18ed..fea52dc46 100644 --- a/test/runtests_slayer_runner.jl +++ b/test/runtests_slayer_runner.jl @@ -153,6 +153,41 @@ @test_throws ArgumentError run_slayer_from_inputs(params, bad_dp, c) end + @testset "delta_prime_to_rs_reference: ψ_N → r_s conversion" begin + # Two surfaces with distinct K and μ; the transform is + # Δ̂_ij = K_i^(1/2+μ_i) · Δ'_ij · K_j^(μ_j−1/2). + K1, mu1 = 0.9, 0.54 + K2, mu2 = 1.2, 0.60 + p1 = SLAYERParameters(; tau=1.0, lu=1e7, c_beta=0.1, D_norm=2.0, + P_perp=20.0, P_tor=10.0, Q_e=-1.0, Q_i=0.5, iota_e=2 / 3, + tauk=1e-4, tau_r=1.0, delta_n=1.0, rs=0.4, R0=1.7, bt=2.0, + sval_r=1.0, eta=2.5e-8, d_beta=4e-3, m=2, n=1, ising=1, + k_ref=K1, mu_mercier=mu1) + p2 = SLAYERParameters(; tau=1.0, lu=1e7, c_beta=0.1, D_norm=2.0, + P_perp=20.0, P_tor=10.0, Q_e=-1.0, Q_i=0.5, iota_e=2 / 3, + tauk=1e-4, tau_r=1.0, delta_n=1.0, rs=0.5, R0=1.7, bt=2.0, + sval_r=1.0, eta=2.5e-8, d_beta=4e-3, m=3, n=1, ising=2, + k_ref=K2, mu_mercier=mu2) + dp = ComplexF64[10.0+1im 2.0-0.5im; 3.0+0im 1.5+2im] + out = Runner.delta_prime_to_rs_reference(dp, [p1, p2]) + # Diagonal carries the scalar K^(2μ) + @test out[1, 1] ≈ K1^(2mu1) * dp[1, 1] + @test out[2, 2] ≈ K2^(2mu2) * dp[2, 2] + # Off-diagonals carry the split row/column factors + @test out[1, 2] ≈ K1^(0.5 + mu1) * K2^(mu2 - 0.5) * dp[1, 2] + @test out[2, 1] ≈ K2^(0.5 + mu2) * K1^(mu1 - 0.5) * dp[2, 1] + # At the slab point μ = 1/2 the diagonal factor is exactly K + pslab = SLAYERParameters(; tau=1.0, lu=1e7, c_beta=0.1, D_norm=2.0, + P_perp=20.0, P_tor=10.0, Q_e=-1.0, Q_i=0.5, iota_e=2 / 3, + tauk=1e-4, tau_r=1.0, delta_n=1.0, rs=0.4, R0=1.7, bt=2.0, + sval_r=1.0, eta=2.5e-8, d_beta=4e-3, m=2, n=1, ising=1, + k_ref=0.8, mu_mercier=0.5) + dp1 = ComplexF64[5.0+0im;;] + @test Runner.delta_prime_to_rs_reference(dp1, [pslab])[1, 1] ≈ 0.8 * 5.0 + # Default parameters (k_ref = 1) give the identity regardless of μ + @test Runner.delta_prime_to_rs_reference(dp, [_mk_params(), _mk_params()]) ≈ dp + end + @testset "run_slayer_from_inputs: coupled mode finds known root" begin # Build a 2-surface problem with a known coupled root by construction. p1 = _mk_params(; rs=0.5, lu=1.0e7, tauk=1.0e-4, Q_e=-1.0, Q_i=0.5, From f8e454cce22cc58be995c826d0435a86ad41e712 Mon Sep 17 00:00:00 2001 From: d-burg Date: Thu, 20 Aug 2026 13:44:24 -0400 Subject: [PATCH 2/3] InnerLayer.SLAYER - NEW FEATURE - Add halfwidth and volume radial-label options to build_slayer_inputs MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit The cylindrical layer model needs a single minor radius; on a shaped torus that label is ambiguous. Add rs_method options :halfwidth (midplane half-chord, the shift-free stand-in for the circular-theory flux label; reproduces it to ~1% on circular benchmarks) and :volume (cylinder-equivalent √(V/2π²R₀), the Rutherford-literature convention) alongside the existing :midplane default and :fsa. Every label feeds S, the r-based shear, W_d, and k_ref together, so each choice is self-consistent by construction. Programmatic API only — not exposed via TOML; the default and all TOML-driven results are unchanged. Co-Authored-By: Claude Opus 5 --- src/InnerLayer/SLAYER/LayerInputs.jl | 29 +++++++++++++++++++++++++--- test/runtests_slayer_inputs.jl | 18 +++++++++++++++++ 2 files changed, 44 insertions(+), 3 deletions(-) diff --git a/src/InnerLayer/SLAYER/LayerInputs.jl b/src/InnerLayer/SLAYER/LayerInputs.jl index 8a408dd29..091dcb49a 100644 --- a/src/InnerLayer/SLAYER/LayerInputs.jl +++ b/src/InnerLayer/SLAYER/LayerInputs.jl @@ -16,7 +16,7 @@ using ..Utilities: KineticProfiles using ...Utilities.NeoclassicalResistivity: NeoResistivityModel, SpitzerModel, coulomb_log_e, nu_star_e -using FastInterpolations: DerivOp +using FastInterpolations: DerivOp, integrate """ surface_minor_radius(equil, psi; theta=0.0) -> Float64 @@ -114,6 +114,20 @@ profiles, without an intermediate file round-trip. a prescribed value. (For `dc_type=:rfitzp` and `:lar`, dgeo_val is not consulted.) - `dc_type` -- `:none` (default), `:lar`, `:rfitzp`, or `:toroidal`. + - `rs_method` -- radial label defining `r_s` for the whole layer stack + (S, r-based shear, W_d, and the Δ' reference-length factor `k_ref` all + follow it together, so every choice is self-consistent). `:midplane` + (default): outboard-midplane chord from the magnetic axis — natural for + comparison with midplane diagnostics. `:halfwidth`: midplane half-chord + `(R_out − R_in)/2` — shift-free, the closest stand-in for the circular + flux label the Fitzpatrick layer formulas are derived in (reproduces the + cylindrical-theory label to ~1% on circular benchmark equilibria). + `:fsa`: θ-mean surface radius. `:volume`: cylinder-equivalent label + `√(V(ψ)/(2π²R₀))`, the Rutherford-literature convention. On shaped + equilibria the labels agree at low-q surfaces (~±3% in growth rate at + q=2) and diverge strongly near the edge, where the slab-layer matching + is label-ambiguous regardless of choice. Not exposed via TOML — + programmatic use only. - `theta` -- poloidal angle at which to measure minor radius (default `0.0`, outboard midplane). - `resistivity_model` -- `SauterNeoModel()` (default), `RedlNeoModel()`, @@ -158,7 +172,9 @@ function build_slayer_inputs(equil, sings, profiles::KineticProfiles; # Minor-radius extractor: `:midplane` = outboard-midplane chord # (original behavior); `:fsa` = θ-mean of √rzphi_rsquared, the - # flux-surface-averaged minor radius. + # flux-surface-averaged minor radius; `:halfwidth` = midplane half-chord + # (R_out − R_in)/2, the shift-free circular flux label; `:volume` = the + # cylinder-equivalent label r_V = √(V(ψ)/(2π²R₀)). _rs_at(ψ) = if rs_method === :fsa integrand(θ) = sqrt(equil.rzphi_rsquared((Float64(ψ), Float64(θ)))) @@ -168,11 +184,18 @@ function build_slayer_inputs(equil, sings, profiles::KineticProfiles; s += integrand((k - 0.5) / N) end s / N + elseif rs_method === :halfwidth + 0.5 * (surface_minor_radius(equil, ψ; theta=0.0) + + surface_minor_radius(equil, ψ; theta=0.5)) + elseif rs_method === :volume + lo = 1e-4 + V = integrate(equil.profiles.dVdpsi_spline, lo, Float64(ψ)) + sqrt(max(V, 0.0) / (2π^2 * equil.ro)) else surface_minor_radius(equil, ψ; theta=theta) end _da_dpsi_at(ψ) = - if rs_method === :fsa + if rs_method !== :midplane # central finite difference on _rs_at h = 1e-5 lo = ψ - h diff --git a/test/runtests_slayer_inputs.jl b/test/runtests_slayer_inputs.jl index 55d6124e9..159406380 100644 --- a/test/runtests_slayer_inputs.jl +++ b/test/runtests_slayer_inputs.jl @@ -158,6 +158,24 @@ @test sl_var[1].P_perp ≈ sl_s[1].P_perp * 6.0 / 2.0 rtol = 1e-10 end + @testset "build_slayer_inputs: rs_method radial labels are self-consistent" begin + sings = [_mk_sing(psi=0.5, q=2.4, q1=1.2, m=2, n=1)] + got = Dict{Symbol,Any}() + for rsm in (:midplane, :halfwidth, :fsa, :volume) + sl = build_slayer_inputs(equil, sings, profiles; bt=2.0, dr_val=0.0, rs_method=rsm) + got[rsm] = sl[1] + @test isfinite(sl[1].rs) && sl[1].rs > 0 + @test isfinite(sl[1].k_ref) && sl[1].k_ref > 0 + @test isfinite(sl[1].sval_r) + end + # The outboard-shifted axis compresses the outboard chord, so the shift-free + # half-chord is at least the outboard axis-to-edge distance. + @test got[:halfwidth].rs >= got[:midplane].rs + # Labels genuinely differ (each self-consistent set has its own rs, S, k_ref) + @test got[:fsa].rs != got[:midplane].rs + @test got[:volume].lu != got[:midplane].lu + 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)] From baf674659134d650187c6bf729ab3a70fae4f324 Mon Sep 17 00:00:00 2001 From: d-burg Date: Thu, 20 Aug 2026 19:11:35 -0400 Subject: [PATCH 3/3] InnerLayer.SLAYER - FEATURE - Add the Fitzpatrick flux radial label and unify the label plumbing MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Factor the rs_method radial-label options out of build_slayer_inputs into one exported radial_label(equil; rs_method, theta) returning (r, dr/dpsi) closures, and add a fifth label :flux — Fitzpatrick's toroidal-flux surface label (Nucl. Fusion 2025, Eq. 30), r = sqrt(2 psi_t / B0) with the g = F/(B0 R0) correction carried. All five labels now use analytic psi-derivatives from the interpolants themselves; the finite-difference stencils (which silently clamped within 1e-4 of the flux boundaries) are removed. The :fsa branch gains the max(r^2, 0) guard needed on extrapolated surfaces. k_ref and the K^(2mu) Delta-prime conversion consume the same closures, so the label choice drives S, the r-based shear, W_d and the Delta-prime reference length together for every option. Default :midplane behavior is unchanged up to the analytic-derivative refinement. Shared infrastructure for the resistive layer-overlap psihigh cap branch, which evaluates its criterion in :flux. Co-Authored-By: Claude Fable 5 --- src/InnerLayer/InnerLayer.jl | 4 +- src/InnerLayer/SLAYER/LayerInputs.jl | 200 ++++++++++++++++----------- src/InnerLayer/SLAYER/SLAYER.jl | 2 +- test/runtests_slayer_inputs.jl | 23 ++- 4 files changed, 140 insertions(+), 89 deletions(-) diff --git a/src/InnerLayer/InnerLayer.jl b/src/InnerLayer/InnerLayer.jl index 43d1fc214..ddbde66c1 100644 --- a/src/InnerLayer/InnerLayer.jl +++ b/src/InnerLayer/InnerLayer.jl @@ -25,7 +25,7 @@ import .GGJ: delta_convergence, solution_profile, asymptotic_profile, q4_surface import .SLAYER: SLAYERModel, SLAYERParameters, slayer_parameters, r_based_shear import .SLAYER: riccati_del_s, slayer_layer_thickness, LayerWidths -import .SLAYER: surface_minor_radius, surface_da_dpsi, build_slayer_inputs +import .SLAYER: surface_minor_radius, surface_da_dpsi, radial_label, build_slayer_inputs export InnerLayerModel, InnerLayerParameters, InnerLayerResponse, solve_inner, solve_inner_profile export GGJ, GGJModel, GGJParameters @@ -37,6 +37,6 @@ export delta_convergence, solution_profile, asymptotic_profile, q4_surface_bench export SLAYER, SLAYERModel, SLAYERParameters, slayer_parameters, r_based_shear export riccati_del_s, slayer_layer_thickness, LayerWidths -export surface_minor_radius, surface_da_dpsi, build_slayer_inputs +export surface_minor_radius, surface_da_dpsi, radial_label, build_slayer_inputs end # module InnerLayer diff --git a/src/InnerLayer/SLAYER/LayerInputs.jl b/src/InnerLayer/SLAYER/LayerInputs.jl index 091dcb49a..c649aeb01 100644 --- a/src/InnerLayer/SLAYER/LayerInputs.jl +++ b/src/InnerLayer/SLAYER/LayerInputs.jl @@ -9,14 +9,14 @@ # Geometry extraction: # - Minor radius at the outboard midplane (θ = 0) via # `equil.rzphi_rsquared((ψ, 0.0))`. -# - `da/dψ` via central finite difference on the same bicubic. +# - `da/dψ` from the interpolant's own analytic ψ-derivative. # - r-based magnetic shear via `r_based_shear(rs, q, q1, da/dψ)` (defined # in LayerParameters.jl). using ..Utilities: KineticProfiles using ...Utilities.NeoclassicalResistivity: NeoResistivityModel, SpitzerModel, coulomb_log_e, nu_star_e -using FastInterpolations: DerivOp, integrate +using FastInterpolations: DerivOp, integrate, cubic_interp, cumulative_integrate, ExtendExtrap """ surface_minor_radius(equil, psi; theta=0.0) -> Float64 @@ -32,33 +32,123 @@ function surface_minor_radius(equil, psi::Real; theta::Real=0.0) end """ - surface_da_dpsi(equil, psi; theta=0.0, h=1e-5) -> Float64 + surface_da_dpsi(equil, psi; theta=0.0) -> Float64 -Central finite-difference approximation of `d(minor radius)/dψ` at `psi`. -Falls back to one-sided differences near the flux-coordinate boundaries -(0 or 1). +Analytic ψ-derivative of the minor radius at `psi` and poloidal angle +`theta`, taken from the `rzphi_rsquared` interpolant's own ψ-derivative as +`da/dψ = (∂r²/∂ψ) / (2a)`. Valid wherever the interpolant is, including +under extrapolation past the ψ grid. Diverges at the magnetic axis, where +`a ~ √ψ`; callers evaluating near `ψ = 0` must check `isfinite`. """ -function surface_da_dpsi(equil, psi::Real; theta::Real=0.0, h::Real=1e-5) - psi_f = Float64(psi) - # Clamp to safe sampling range within (0, 1) - eps_edge = 10 * h - lo = psi_f - h - hi = psi_f + h - if lo < eps_edge - # one-sided forward - a0 = surface_minor_radius(equil, max(psi_f, eps_edge); theta=theta) - a1 = surface_minor_radius(equil, max(psi_f, eps_edge) + h; theta=theta) - return (a1 - a0) / h - elseif hi > 1.0 - eps_edge - # one-sided backward - a0 = surface_minor_radius(equil, min(psi_f, 1.0 - eps_edge) - h; theta=theta) - a1 = surface_minor_radius(equil, min(psi_f, 1.0 - eps_edge); theta=theta) - return (a1 - a0) / h +function surface_da_dpsi(equil, psi::Real; theta::Real=0.0) + return _da_dpsi_at_theta(equil, Float64(psi), Float64(theta)) +end + +# d(√r²)/dψ at one (ψ, θ) from the interpolant's own ψ-derivative. Shared by all radial-label +# conventions so none carries its own stencil. +@inline function _da_dpsi_at_theta(equil, psi::Float64, theta::Float64) + r_sq = equil.rzphi_rsquared((psi, theta)) + a = sqrt(max(r_sq, 0.0)) + a > 0 || return Inf # magnetic axis: a ~ √ψ, so da/dψ genuinely diverges + return equil.rzphi_rsquared((psi, theta); deriv=DerivOp(1, 0)) / (2a) +end + +""" + radial_label(equil; rs_method=:midplane, theta=0.0) -> (r_at, dr_dpsi_at) + +Build the pair of closures `r_at(ψ)` and `dr_dpsi_at(ψ)` defining one radial +label for the layer stack. Both closures must come from the same label, +because the r-based shear `(r/q)dq/dr`, `τ_R = μ₀r²/η`, `τ_E = r²/χ`, the +`d_β/r` normalization, the Δ' reference-length factor `k_ref = r_s/(da/dψ)`, +and any metre-to-ψ width conversion all have to live in one coordinate; a +mismatched `r` and `dr/dψ` silently corrupts every one of them. All +derivatives are analytic — no label carries a finite-difference stencil. + +# Labels + + - `:midplane` -- outboard-midplane chord from the magnetic axis at `theta` + (historical default), natural for comparison with midplane diagnostics. + - `:halfwidth` -- midplane half-chord, the mean of the outboard and inboard + chords at `θ = 0` and `θ = 0.5`; shift-free. Coincides with the flux label + on circular equilibria but is its own convention on shaped ones. + - `:fsa` -- θ-mean surface radius, a 128-point midpoint mean of the local + minor radius; the closest geometric approximation to `:flux` at interior + surfaces of shaped equilibria. + - `:volume` -- cylinder-equivalent label `√(V(ψ)/(2π²R₀))`, the + Rutherford-literature convention. + - `:flux` -- toroidal-flux label. Fitzpatrick, Nucl. Fusion (2025), + Eq. 30: `dψ_p/dr = B₀ r g/q` integrates to `ψ_t = B₀r²/2`, so + `r = √(2ψ_t/B₀)` with `ψ_t = psio·∫₀^ψ (q/g) dψ′` and `g = F/(B₀R₀)`. + Defined from flux alone, it carries no circular-cross-section assumption, + and its derivative `dr/dψ ∝ q` grows toward a separatrix where the + geometric labels' `da/dψ` collapses. + +On shaped equilibria the labels agree at low-q surfaces and diverge strongly +near the edge, where the slab-layer matching is label-ambiguous regardless of +choice. The label is selected programmatically; it is not exposed via TOML. +""" +function radial_label(equil; rs_method::Symbol=:midplane, theta::Real=0.0) + theta_f = Float64(theta) + + _flux_r, _flux_dr = if rs_method === :flux + b0f = Float64(equil.params.b0) + R0f = Float64(equil.ro) + psiof = Float64(equil.psio) + xs_f = collect(Float64, equil.profiles.xs) + # g = F/(B0 R0) departs from 1 by ~3% on a DIII-D-like deck and ~0.8% on a circular + # one, so carry it rather than assuming g = 1: r² = 2∫(q/g)dψ_p/B0. + _g_at(x) = Float64(equil.profiles.F_spline(x)) / (2π * b0f * R0f) + qg = [Float64(equil.profiles.q_spline(x)) / _g_at(x) for x in xs_f] + Phi = collect(Float64, cumulative_integrate(cubic_interp(xs_f, qg))) + r_knots = sqrt.(max.(2 .* psiof .* Phi ./ b0f, 0.0)) + rspl = cubic_interp(xs_f, r_knots; extrap=ExtendExtrap()) + (ψ -> Float64(rspl(Float64(ψ))), + ψ -> psiof * Float64(equil.profiles.q_spline(Float64(ψ))) / + (b0f * _g_at(Float64(ψ)) * max(Float64(rspl(Float64(ψ))), eps()))) else - a_plus = surface_minor_radius(equil, psi_f + h; theta=theta) - a_minus = surface_minor_radius(equil, psi_f - h; theta=theta) - return (a_plus - a_minus) / (2h) + (nothing, nothing) end + + _a_at(ψ, θ) = sqrt(max(equil.rzphi_rsquared((Float64(ψ), Float64(θ))), 0.0)) + + _rs_at(ψ) = + if rs_method === :fsa + N = 128 + s = 0.0 + @inbounds for k in 1:N + s += _a_at(ψ, (k - 0.5) / N) + end + s / N + elseif rs_method === :halfwidth + 0.5 * (_a_at(ψ, 0.0) + _a_at(ψ, 0.5)) + elseif rs_method === :volume + V = integrate(equil.profiles.dVdpsi_spline, 1e-4, Float64(ψ)) + sqrt(max(V, 0.0) / (2π^2 * equil.ro)) + elseif rs_method === :flux + _flux_r(ψ) + else + surface_minor_radius(equil, ψ; theta=theta_f) + end + + _da_dpsi_at(ψ) = + if rs_method === :fsa + N = 128 + s = 0.0 + @inbounds for k in 1:N + s += _da_dpsi_at_theta(equil, Float64(ψ), (k - 0.5) / N) + end + s / N + elseif rs_method === :halfwidth + 0.5 * (_da_dpsi_at_theta(equil, Float64(ψ), 0.0) + _da_dpsi_at_theta(equil, Float64(ψ), 0.5)) + elseif rs_method === :volume + Float64(equil.profiles.dVdpsi_spline(ψ)) / (4π^2 * equil.ro * max(_rs_at(ψ), eps())) + elseif rs_method === :flux + _flux_dr(ψ) + else + _da_dpsi_at_theta(equil, Float64(ψ), theta_f) + end + + return (_rs_at, _da_dpsi_at) end """ @@ -114,19 +204,11 @@ profiles, without an intermediate file round-trip. a prescribed value. (For `dc_type=:rfitzp` and `:lar`, dgeo_val is not consulted.) - `dc_type` -- `:none` (default), `:lar`, `:rfitzp`, or `:toroidal`. - - `rs_method` -- radial label defining `r_s` for the whole layer stack - (S, r-based shear, W_d, and the Δ' reference-length factor `k_ref` all - follow it together, so every choice is self-consistent). `:midplane` - (default): outboard-midplane chord from the magnetic axis — natural for - comparison with midplane diagnostics. `:halfwidth`: midplane half-chord - `(R_out − R_in)/2` — shift-free, the closest stand-in for the circular - flux label the Fitzpatrick layer formulas are derived in (reproduces the - cylindrical-theory label to ~1% on circular benchmark equilibria). - `:fsa`: θ-mean surface radius. `:volume`: cylinder-equivalent label - `√(V(ψ)/(2π²R₀))`, the Rutherford-literature convention. On shaped - equilibria the labels agree at low-q surfaces (~±3% in growth rate at - q=2) and diverge strongly near the edge, where the slab-layer matching - is label-ambiguous regardless of choice. Not exposed via TOML — + - `rs_method` -- radial label defining `r_s` for the whole layer stack: + `:midplane` (default), `:halfwidth`, `:fsa`, `:volume`, or `:flux`. See + [`radial_label`](@ref) for the definitions. S, the r-based shear, W_d, + and the Δ' reference-length factor `k_ref` all follow the choice + together, so every option is self-consistent. Not exposed via TOML — programmatic use only. - `theta` -- poloidal angle at which to measure minor radius (default `0.0`, outboard midplane). @@ -170,47 +252,7 @@ function build_slayer_inputs(equil, sings, profiles::KineticProfiles; Float64(bt(ψ)) end - # Minor-radius extractor: `:midplane` = outboard-midplane chord - # (original behavior); `:fsa` = θ-mean of √rzphi_rsquared, the - # flux-surface-averaged minor radius; `:halfwidth` = midplane half-chord - # (R_out − R_in)/2, the shift-free circular flux label; `:volume` = the - # cylinder-equivalent label r_V = √(V(ψ)/(2π²R₀)). - _rs_at(ψ) = - if rs_method === :fsa - integrand(θ) = sqrt(equil.rzphi_rsquared((Float64(ψ), Float64(θ)))) - N = 128 - s = 0.0 - @inbounds for k in 1:N - s += integrand((k - 0.5) / N) - end - s / N - elseif rs_method === :halfwidth - 0.5 * (surface_minor_radius(equil, ψ; theta=0.0) + - surface_minor_radius(equil, ψ; theta=0.5)) - elseif rs_method === :volume - lo = 1e-4 - V = integrate(equil.profiles.dVdpsi_spline, lo, Float64(ψ)) - sqrt(max(V, 0.0) / (2π^2 * equil.ro)) - else - surface_minor_radius(equil, ψ; theta=theta) - end - _da_dpsi_at(ψ) = - if rs_method !== :midplane - # central finite difference on _rs_at - h = 1e-5 - lo = ψ - h - hi = ψ + h - eps_edge = 10h - if lo < eps_edge - (_rs_at(max(ψ, eps_edge) + h) - _rs_at(max(ψ, eps_edge))) / h - elseif hi > 1.0 - eps_edge - (_rs_at(min(ψ, 1.0 - eps_edge)) - _rs_at(min(ψ, 1.0 - eps_edge) - h)) / h - else - (_rs_at(ψ + h) - _rs_at(ψ - h)) / (2h) - end - else - surface_da_dpsi(equil, ψ; theta=theta) - end + _rs_at, _da_dpsi_at = radial_label(equil; rs_method=rs_method, theta=theta) # Per-surface ω_*e, ω_*i (diamagnetic frequencies) from spline # derivatives. When `compute_omega_star=true` we override any ω_*e/ω_*i diff --git a/src/InnerLayer/SLAYER/SLAYER.jl b/src/InnerLayer/SLAYER/SLAYER.jl index ece897c15..033f71eca 100644 --- a/src/InnerLayer/SLAYER/SLAYER.jl +++ b/src/InnerLayer/SLAYER/SLAYER.jl @@ -56,7 +56,7 @@ include("LayerInputs.jl") export SLAYERModel, SLAYERParameters, slayer_parameters export r_based_shear export riccati_del_s, slayer_layer_thickness, LayerWidths -export surface_minor_radius, surface_da_dpsi, build_slayer_inputs +export surface_minor_radius, surface_da_dpsi, radial_label, build_slayer_inputs export NeoResistivityModel, SpitzerModel, SpitzerHarmModel, SauterNeoModel, RedlNeoModel end # module SLAYER diff --git a/test/runtests_slayer_inputs.jl b/test/runtests_slayer_inputs.jl index 159406380..a5a6da414 100644 --- a/test/runtests_slayer_inputs.jl +++ b/test/runtests_slayer_inputs.jl @@ -42,8 +42,8 @@ @test r1 > 0 end - @testset "surface_da_dpsi: FD agrees with numerical derivative" begin - # Reference via a tighter FD + @testset "surface_da_dpsi: agrees with FD reference" begin + # The analytic ψ-derivative must reproduce a tight FD of the minor radius. for psi in (0.1, 0.4, 0.7) h_ref = 1e-4 r_p = surface_minor_radius(equil, psi + h_ref) @@ -53,16 +53,24 @@ end end - @testset "surface_da_dpsi: one-sided near boundaries" begin - # Near ψ=0 and ψ=1, the function falls back to one-sided FD and - # should still produce a finite positive number (minor radius is - # still increasing). + @testset "surface_da_dpsi: finite near the boundaries" begin + # The analytic form needs no clamping: near ψ=0 and ψ=1 it still returns a + # finite positive number (large near the axis, where a ~ √ψ). d_near_axis = surface_da_dpsi(equil, 1e-6) d_near_edge = surface_da_dpsi(equil, 1.0 - 1e-6) @test isfinite(d_near_axis) && d_near_axis > 0 @test isfinite(d_near_edge) && d_near_edge > 0 end + @testset "radial_label: analytic derivatives match FD" begin + for rsm in (:midplane, :halfwidth, :fsa, :volume, :flux) + r_at, dr_at = radial_label(equil; rs_method=rsm) + h = 1e-5 + ref = (r_at(0.5 + h) - r_at(0.5 - h)) / (2h) + @test dr_at(0.5) ≈ ref rtol = 1e-3 + end + end + @testset "build_slayer_inputs: returns correct per-surface data" begin sings = [_mk_sing(psi=0.3, q=2.0, q1=1.5, m=2, n=1), _mk_sing(psi=0.6, q=3.0, q1=2.5, m=3, n=1)] @@ -161,7 +169,7 @@ @testset "build_slayer_inputs: rs_method radial labels are self-consistent" begin sings = [_mk_sing(psi=0.5, q=2.4, q1=1.2, m=2, n=1)] got = Dict{Symbol,Any}() - for rsm in (:midplane, :halfwidth, :fsa, :volume) + for rsm in (:midplane, :halfwidth, :fsa, :volume, :flux) sl = build_slayer_inputs(equil, sings, profiles; bt=2.0, dr_val=0.0, rs_method=rsm) got[rsm] = sl[1] @test isfinite(sl[1].rs) && sl[1].rs > 0 @@ -174,6 +182,7 @@ # Labels genuinely differ (each self-consistent set has its own rs, S, k_ref) @test got[:fsa].rs != got[:midplane].rs @test got[:volume].lu != got[:midplane].lu + @test got[:flux].rs != got[:midplane].rs end @testset "build_slayer_inputs: dc_type propagates and dr_val activates offset" begin