diff --git a/docs/src/stability.md b/docs/src/stability.md index 91f953fc1..814a5207a 100644 --- a/docs/src/stability.md +++ b/docs/src/stability.md @@ -379,7 +379,7 @@ do not resolve the H-mode pedestal transport barrier, where power-balance χ fal | Dataset | Contents | |---|---| -| `psi`, `m`, `n`, `r_s` | scanned surfaces and their minor radius | +| `psi`, `m`, `n`, `r_s` | scanned surfaces and their minor radius √(2ψ_t/B₀) | | `delta_s_abs`, `width_delta_s` | Riccati layer thickness, in metres and as a Δψ width | | `width_visco`, `width_dr` | visco-resistive and diffusive-resistive (Eq. 100) Δψ widths; `width_dr` sets the cut | | `extrapolated` | 1 where the surface lies beyond the equilibrium grid | diff --git a/src/InnerLayer/SLAYER/LayerInputs.jl b/src/InnerLayer/SLAYER/LayerInputs.jl index 97a6c28db..0e72705b7 100644 --- a/src/InnerLayer/SLAYER/LayerInputs.jl +++ b/src/InnerLayer/SLAYER/LayerInputs.jl @@ -7,10 +7,9 @@ # round-trip. # # Geometry extraction: -# - Minor radius at the outboard midplane (θ = 0) via -# `equil.rzphi_rsquared((ψ, 0.0))`. -# - `da/dψ` from the interpolant's own analytic ψ-derivative. -# - r-based magnetic shear via `r_based_shear(rs, q, q1, da/dψ)` (defined +# - r_s and dr/dψ from one radial label (`radial_label`; by default the +# toroidal-flux label r = √(2ψ_t/B₀)), both analytic. +# - r-based magnetic shear via `r_based_shear(rs, q, q1, dr/dψ)` (defined # in LayerParameters.jl). using ..Utilities: KineticProfiles @@ -54,7 +53,7 @@ end end """ - radial_label(equil; rs_method=:midplane, theta=0.0) -> (r_at, dr_dpsi_at) + radial_label(equil; rs_method=:flux, 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, @@ -91,7 +90,7 @@ 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, q_beyond=nothing) +function radial_label(equil; rs_method::Symbol=:flux, theta::Real=0.0, q_beyond=nothing) theta_f = Float64(theta) _flux_r, _flux_dr = if rs_method === :flux @@ -294,14 +293,14 @@ profiles, without an intermediate file round-trip. without one it falls back to the cylindrical `n·|s|/R₀` with a warning. - `dc_type` -- `:none` (default), `:lar`, `:fitzpatrick`, or `:toroidal`. - `rs_method` -- radial label defining `r_s` for the whole layer stack: - `:midplane` (default), `:halfwidth`, `:fsa`, `:volume`, or `:flux`. See + `:flux` (default), `:midplane`, `:halfwidth`, `:fsa`, or `:volume`. 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). + - `theta` -- poloidal angle of the minor-radius chord, used only by + `rs_method=:midplane` (default `0.0`, outboard midplane). - `resistivity_model` -- `SauterNeoModel()` (default), `RedlNeoModel()`, `SpitzerModel()`, or `SpitzerHarmModel()` (legacy Fitzpatrick σ_∥). @@ -317,7 +316,7 @@ profiles, without an intermediate file round-trip. function build_slayer_inputs(equil, sings, profiles::KineticProfiles; bt=nothing, R0=nothing, - rs_method::Symbol=:midplane, + rs_method::Symbol=:flux, q_beyond=nothing, mu_i::Real=2.0, zeff=1.0, diff --git a/src/Tearing/LayerOverlap.jl b/src/Tearing/LayerOverlap.jl index 4fc96d21e..4d48ecf0e 100644 --- a/src/Tearing/LayerOverlap.jl +++ b/src/Tearing/LayerOverlap.jl @@ -151,9 +151,9 @@ Riccati solve; `delta_s` is solved per surface as a reported comparison only. ψ = 1 on the real q. Geometry needs no clamp of its own -- `radial_label` derivatives are analytic and valid under extrapolation - `extrapolate` -- search past `psihigh_safe` on the separatrix q-law (default `true`) - - `theta` -- poloidal angle for the minor-radius chord (default `0.0`) + - `theta` -- poloidal angle for the minor-radius chord, used only by `rs_method=:midplane` (default `0.0`) - `rs_method` -- radial label the widths and their ψ-conversion Jacobian are taken in - (default `:midplane`; see `radial_label`). The caller wiring this scan into the domain + (default `:flux`; see `radial_label`). The caller wiring this scan into the domain truncation passes `:flux`, the label Fitzpatrick Eq. (100) is anchored to - `bt`, `mu_i`, `zeff`, `chi_perp_e`, `chi_perp_i`, `chi_tor`, `tau_E`, `P_perp_model`, `resistivity_model`, `lnLambda_form` -- passed through to `build_slayer_inputs`. The defaults are placeholders: a caller @@ -172,7 +172,7 @@ function resistive_layer_overlap(equil, profiles::KineticProfiles; psi_gap::Real=1e-9, extrapolate::Bool=true, theta::Real=0.0, - rs_method::Symbol=:midplane, + rs_method::Symbol=:flux, bt=nothing, mu_i::Real=2.0, zeff=1.0, diff --git a/src/Tearing/Runner/Control.jl b/src/Tearing/Runner/Control.jl index 0c49e0ccb..0e0d6368a 100644 --- a/src/Tearing/Runner/Control.jl +++ b/src/Tearing/Runner/Control.jl @@ -104,8 +104,10 @@ constructor. offset (Δ_crit ≡ 0) for every `D_c_type`; `D_geo = 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) + - `theta_sample` -- poloidal angle of the minor-radius chord of the + `:midplane` radial label, used only by API callers that pass + `rs_method = :midplane`; the default toroidal-flux label does not depend + on angle, so a nonzero value warns (default 0.0) - `resistivity_model` -- η closure setting τ_R = μ₀r_s²/η: `:sauter` (default, neoclassical F_33), `:redl`, `:spitzer` (Sauter 18a, no trapped correction), or `:spitzer_harm` (legacy Fitzpatrick σ_∥) @@ -242,13 +244,19 @@ const _VALID_LNLAMBDA_FORMS = (:nrl, :sauter, :wesson) const _VALID_ROOT_CLASSIFICATIONS = (:threshold, :winding) function validate(ctrl::SLAYERControl) + iszero(ctrl.theta_sample) || @warn "SLAYERControl: theta_sample = $(ctrl.theta_sample) has no effect: the layer quantities " * + "use the toroidal-flux radial label, which does not depend on poloidal angle." maxlog = 1 ctrl.inner_model in _VALID_INNER_MODELS || ctrl.inner_model in _UNSUPPORTED_INNER_MODELS || throw(ArgumentError("SLAYERControl: inner_model=$(ctrl.inner_model) " * "not in $(_VALID_INNER_MODELS)")) # A disabled section is never solved, so it may keep a GGJ model while serving only its profiles. !ctrl.enabled || ctrl.inner_model in _VALID_INNER_MODELS || - throw(ArgumentError("SLAYERControl: inner_model=$(ctrl.inner_model) is not currently supported because GGJ tearing " * - "growth rates are not implemented; use inner_model=slayer_fitzpatrick")) + throw( + ArgumentError( + "SLAYERControl: inner_model=$(ctrl.inner_model) is not currently supported because GGJ tearing " * + "growth rates are not implemented; use inner_model=slayer_fitzpatrick" + ) + ) ctrl.scan_mode in _VALID_SCAN_MODES || throw(ArgumentError("SLAYERControl: scan_mode=$(ctrl.scan_mode) " * "not in $(_VALID_SCAN_MODES)")) diff --git a/src/Tearing/Runner/HDF5Output.jl b/src/Tearing/Runner/HDF5Output.jl index 9d81d3a6c..6cf61d92a 100644 --- a/src/Tearing/Runner/HDF5Output.jl +++ b/src/Tearing/Runner/HDF5Output.jl @@ -114,7 +114,7 @@ const TEARING_H5_ANNOTATIONS = [ (; long_name="real E×B Doppler offset ΔRe(Q) = −τ_k·n·Ω_E applied to each surface's inner-layer Q in the coupled determinant (0 in uncoupled mode)", dims=("surface",)), "PerSurface/tau_R" => (; long_name="resistive diffusion time τ_R = μ₀r_s²/η per surface", units="s", dims=("surface",)), "PerSurface/Delta_prime_norm" => (; long_name="Δ'-normalization factor S^(1/3)/r_s per surface", units="1/m", dims=("surface",)), - "PerSurface/rs" => (; long_name="minor radius of each rational surface", units="m", dims=("surface",)), + "PerSurface/rs" => (; long_name="minor radius √(2ψ_t/B₀) of each rational surface (toroidal-flux label)", units="m", dims=("surface",)), "PerSurface/R0" => (; long_name="major radius", units="m", dims=("surface",)), "PerSurface/bt" => (; long_name="toroidal field", units="T", dims=("surface",)), "PerSurface/sval_r" => (; long_name="r-based magnetic shear r_s·(dq/dr)/q (Fitzpatrick convention)", dims=("surface",)), diff --git a/test/runtests_layer_overlap_fitzpatrick.jl b/test/runtests_layer_overlap_fitzpatrick.jl index b185ae5d9..de6d8bb01 100644 --- a/test/runtests_layer_overlap_fitzpatrick.jl +++ b/test/runtests_layer_overlap_fitzpatrick.jl @@ -226,14 +226,14 @@ end @test binding([s1]) === s1 # The scan must see the run's own plasma: every [SLAYER] setting the layer width depends on - # reaches it, and changing one changes the widths. + # reaches it, and changing one changes the widths. theta_sample is passed through but warns: the flux label ignores it. trunc_dir = joinpath(dirname(@__DIR__), "examples", "DIIID-like_truncation_example") trunc_deck = TOML.parsefile(joinpath(trunc_dir, "gpec.toml")) slayer_section = trunc_deck["SLAYER"] kin = Equilibrium.read_kinetic_file(joinpath(trunc_dir, trunc_deck["Equilibrium"]["kinetic_file"])) slayer_inputs(over) = Dict{String,Any}("SLAYER" => merge(slayer_section, over)) overlap_inputs(over; kinetic=kin) = GeneralizedPerturbedEquilibrium._overlap_slayer_inputs(slayer_inputs(over), kinetic) - oi = overlap_inputs(Dict("ion_mass" => 3.0, "bt" => 1.7, "theta_sample" => 0.3, "lnLambda_form" => "sauter")) + oi = @test_logs (:warn, r"theta_sample") match_mode = :any overlap_inputs(Dict("ion_mass" => 3.0, "bt" => 1.7, "theta_sample" => 0.3, "lnLambda_form" => "sauter")) @test (oi.mu_i, oi.bt, oi.theta, oi.lnLambda_form) == (3.0, 1.7, 0.3, :sauter) # A disabled section feeds the scan whatever its inner model, and missing kinetic profiles only skip the scan; diff --git a/test/runtests_slayer_inputs.jl b/test/runtests_slayer_inputs.jl index 3561d064b..a43d98c53 100644 --- a/test/runtests_slayer_inputs.jl +++ b/test/runtests_slayer_inputs.jl @@ -94,9 +94,10 @@ @test sl[1].R0 ≈ equil.ro @test sl[1].bt == 2.0 - # Minor radius and r-based shear recovered from the equilibrium - rs1 = surface_minor_radius(equil, 0.3) - da1 = surface_da_dpsi(equil, 0.3) + # Minor radius and r-based shear recovered from the equilibrium in the default, toroidal-flux label + r_at, dr_at = radial_label(equil; rs_method=:flux) + rs1 = r_at(0.3) + da1 = dr_at(0.3) @test sl[1].rs ≈ rs1 @test sl[1].sval_r ≈ rs1 * 1.5 / (2.0 * da1)