diff --git a/docs/development/references.md b/docs/development/references.md index e19b92e32..a815ca6f8 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`) 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) diff --git a/regression-harness/cases/diiid_slayer_n1.toml b/regression-harness/cases/diiid_slayer_n1.toml index 90edf2e09..c43964ea4 100644 --- a/regression-harness/cases/diiid_slayer_n1.toml +++ b/regression-harness/cases/diiid_slayer_n1.toml @@ -109,6 +109,17 @@ 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 +class = "physics_converged" + # 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..36f978321 --- /dev/null +++ b/regression-harness/cases/diiid_slayer_n1_toroidal.toml @@ -0,0 +1,70 @@ +# 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-Δ 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" +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" +type = "real_vector" +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. +[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 +class = "physics_converged" + +[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 +class = "physics_converged" + +[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 +class = "physics_converged" + +[quantities.runtime] +h5path = "" +type = "runtime" +extract = "value" +label = "Runtime (s)" +noise_threshold = 0.0 +order = 999 +class = "diagnostic" diff --git a/src/ForceFreeStates/Surfaces/ResistEval.jl b/src/ForceFreeStates/Surfaces/ResistEval.jl index 09bd6dd80..74224c796 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,25 +33,27 @@ ResistGeometry 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_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 | +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 | +|:--------------------- |:-------------------------------------------------------- | +| `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²/‖∇ψ_N‖²⟩ — needed for τ_R | +| `avg_bsq` | ⟨B²⟩ — needed for τ_R | +| `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 | +| `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. @@ -69,6 +73,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 @@ -83,10 +88,10 @@ 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. +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. # Arguments @@ -99,6 +104,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, @@ -107,47 +113,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) - - # Build the 6 GGJ θ-integrands plus a 7th (B) for the neoclassical - # resistivity f_t calculation, and accumulate running extrema of + 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, 7) - 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) @@ -163,6 +169,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 # Snap the repeated endpoint exactly equal to the start @@ -182,7 +189,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])) @@ -192,9 +199,9 @@ 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, + p, p1, v1 ) end @@ -206,8 +213,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/InnerLayer.jl b/src/InnerLayer/InnerLayer.jl index ddbde66c1..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 +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 +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 252e81ae9..6495c0fca 100644 --- a/src/InnerLayer/SLAYER/LayerInputs.jl +++ b/src/InnerLayer/SLAYER/LayerInputs.jl @@ -154,6 +154,72 @@ 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`. + +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 Δ' 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) + 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)) + return alpha, lambda +end + """ build_slayer_inputs(equil, sings, profiles; …) -> Vector{SLAYERParameters} @@ -190,22 +256,28 @@ 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), 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 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. 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 @@ -313,7 +385,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. @@ -332,27 +404,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 @@ -361,13 +412,46 @@ 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 + + # 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, + 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 " * - "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)) @@ -383,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 ae875dea4..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-Hastie-Helander 2015 Eq. 59 geometric factor (0 unless supplied) | -| `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 @@ -130,27 +132,36 @@ 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, + 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, @@ -212,7 +223,11 @@ 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`) + - `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, @@ -266,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, @@ -380,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, @@ -390,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 033f71eca..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 +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 1a12ff7e6..4285a5576 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) @@ -41,9 +41,12 @@ 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 - scalar only to override the auto-derivation; an explicit `0.0` disables the - critical-Δ offset (Δ_crit ≡ 0) + 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 `dr_val = 0.0` disables the critical-Δ + offset (Δ_crit ≡ 0) for every `dc_type`; `dgeo_val = 0.0` does so only for + `: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/src/Tearing/Runner/HDF5Output.jl b/src/Tearing/Runner/HDF5Output.jl index b0fc8d2e4..6942a0564 100644 --- a/src/Tearing/Runner/HDF5Output.jl +++ b/src/Tearing/Runner/HDF5Output.jl @@ -90,10 +90,11 @@ 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 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-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 +107,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",)), diff --git a/test/runtests_slayer_inputs.jl b/test/runtests_slayer_inputs.jl index fce001036..ba071df61 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. @@ -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; @@ -205,6 +205,68 @@ @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 + # 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 + + # 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 + + # 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 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 + + # 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 sl = build_slayer_inputs(equil, SingType[], profiles; bt=2.0) @test sl isa Vector{SLAYERParameters} @@ -216,8 +278,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) 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 5bbcb25d2..274fc7248 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,21 +73,91 @@ 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₀) 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, 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) + 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) + 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) + 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, 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 + @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 + + @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, 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, + 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) + 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 + end end