diff --git a/regression-harness/cases/diiid_slayer_n1_coupled.toml b/regression-harness/cases/diiid_slayer_n1_coupled.toml new file mode 100644 index 000000000..f3f095698 --- /dev/null +++ b/regression-harness/cases/diiid_slayer_n1_coupled.toml @@ -0,0 +1,97 @@ +# Regression case: DIII-D-like H-mode n=1 SLAYER coupled-determinant growth rate. +# Reuses the n=1 SLAYER deck via [overrides] in coupled mode, the only harness coverage of the +# multi-surface determinant: its inter-surface Q normalization and the per-surface E×B Doppler +# shift taken from the kinetic file's omega_E. The Re(Q) window is widened because that shift +# moves the coupled root to the lab frame, outside the deck's [-2, 2]. +# 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_coupled" +description = "DIII-D-like H-mode equilibrium, n=1, SLAYER coupled multi-surface determinant with kinetic-file E×B rotation (reuses the n=1 deck via overrides)" +example_dir = "examples/DIIID-like_SLAYER_example" + +# Run the shared SLAYER deck through the coupled determinant. +[overrides] +"SLAYER.coupling_mode" = "coupled" +"SLAYER.Q_re_range" = [-4.0, 6.0] + +# Per-surface inputs to the determinant +[quantities.slayer_tauk] +h5path = "Tearing/PerSurface/tau_k" +type = "real_vector" +extract = "all_real" +label = "SLAYER tauk" +noise_threshold = 1e-12 +order = 10 +class = "physics_converged" + +# Applied E×B Doppler offset ΔRe(Q) = −τ_k·n·Ω_E per surface +[quantities.slayer_q_shift] +h5path = "Tearing/PerSurface/q_shift" +type = "real_vector" +extract = "all_real" +label = "SLAYER q_shift" +noise_threshold = 1e-10 +order = 11 +class = "physics_converged" + +# Coupled eigenvalue (one root for the whole determinant). The coupled path has no root +# polishing, so the thresholds match the uncoupled case's absolute, modestly loose floors. +# Pin Q, ω and γ with a measured tolerance, not the class default: re-extracting one set of scan +# samples under different triangulations gave γ = 190.646 and 191.888 s⁻¹ (0.65 %). ω and Q +# need the same measurement. +[quantities.slayer_Q_coupled] +h5path = "Tearing/Roots/Q_root" +type = "complex_vector" +extract = "all_complex" +label = "SLAYER coupled Q_root" +noise_threshold = 1e-4 +order = 30 +class = "physics_converged" + +[quantities.slayer_omega_coupled] +h5path = "Tearing/Roots/omega" +type = "real_vector" +extract = "all_real" +label = "SLAYER coupled ω" +noise_threshold = 1.0 +order = 32 +class = "physics_converged" + +[quantities.slayer_gamma_coupled] +h5path = "Tearing/Roots/gamma" +type = "real_vector" +extract = "all_real" +label = "SLAYER coupled γ" +noise_threshold = 1e-1 +order = 33 +class = "physics_converged" + +# no_root flag (1 = extraction failed); guards that the coupled root stays inside the window. +[quantities.slayer_no_root_coupled] +h5path = "Tearing/Roots/no_root" +type = "real_vector" +extract = "all_real" +label = "SLAYER coupled no_root flag" +noise_threshold = 0 +order = 34 +class = "physics_converged" + +# Settings (catches accidental config drift) +[quantities.slayer_enabled] +h5path = "Tearing/enabled" +type = "int_scalar" +extract = "value" +label = "SLAYER enabled flag" +noise_threshold = 0 +order = 90 +class = "topological" + +[quantities.runtime] +h5path = "" +type = "runtime" +extract = "value" +label = "Runtime (s)" +noise_threshold = 0.0 +order = 999 +class = "diagnostic" diff --git a/src/InnerLayer/SLAYER/LayerInputs.jl b/src/InnerLayer/SLAYER/LayerInputs.jl index 252e81ae9..ea3c89189 100644 --- a/src/InnerLayer/SLAYER/LayerInputs.jl +++ b/src/InnerLayer/SLAYER/LayerInputs.jl @@ -238,7 +238,6 @@ function build_slayer_inputs(equil, sings, profiles::KineticProfiles; dgeo_val=nothing, dc_type::Symbol=:none, theta::Real=0.0, - compute_omega_star::Bool=true, resistivity_model::NeoResistivityModel=SauterNeoModel(), lnLambda_form::Symbol=:nrl) R0_use = R0 === nothing ? equil.ro : Float64(R0) @@ -291,9 +290,8 @@ function build_slayer_inputs(equil, sings, profiles::KineticProfiles; n_res = sing.n[1] prof = profiles(psi) - # Take ω_*e, ω_*i from the spline derivatives, or from `profiles` when the caller - # supplies them directly. `run_slayer` supplies zeros, so the latter is a library path. - ω_e_use, ω_i_use = compute_omega_star ? _omega_star_at(psi, n_res) : (prof.omega_e, prof.omega_i) + # ω_*e, ω_*i from the density and temperature spline derivatives, at this surface's n. + ω_e_use, ω_i_use = _omega_star_at(psi, n_res) # Pull geometric trapped-fraction inputs from ResistGeometry when # available (populated by ForceFreeStates.resist_eval_all!); else @@ -375,7 +373,7 @@ function build_slayer_inputs(equil, sings, profiles::KineticProfiles; 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, + omega_e=ω_e_use, omega_i=ω_i_use, qval=q, sval_r=sval_r, bt=_bt_at(psi), rs=rs, R0=R0_use, mu_i=mu_i, zeff=zeff, chi_perp=_eval(chi_perp, psi), diff --git a/src/InnerLayer/SLAYER/LayerParameters.jl b/src/InnerLayer/SLAYER/LayerParameters.jl index ae875dea4..a18dbc81a 100644 --- a/src/InnerLayer/SLAYER/LayerParameters.jl +++ b/src/InnerLayer/SLAYER/LayerParameters.jl @@ -177,7 +177,7 @@ function _solve_dc_tmp(; dc_type::Symbol, dr_val::Real, dgeo_val::Real, end """ - slayer_parameters(; n_e, t_e, t_i, omega, omega_e, omega_i, + slayer_parameters(; n_e, t_e, t_i, omega_e, omega_i, qval, sval_r, bt, rs, R0, mu_i, zeff, chi_perp, chi_tor, m, n, @@ -199,7 +199,6 @@ parametrization (P_perp/P_tor/D_norm; the older magnetic/electron Prandtl - `n_e` -- electron density [m⁻³] - `t_e` -- electron temperature [eV] - `t_i` -- ion temperature [eV] - - `omega` -- toroidal rotation frequency at the surface [rad/s] - `omega_e` -- electron diamagnetic frequency [rad/s] - `omega_i` -- ion diamagnetic frequency [rad/s] - `qval` -- safety factor q at the surface @@ -260,7 +259,7 @@ dispersion relation. """ function slayer_parameters(; n_e::Real, t_e::Real, t_i::Real, - omega::Real, omega_e::Real, omega_i::Real, + omega_e::Real, omega_i::Real, qval::Real, sval_r::Real, bt::Real, rs::Real, R0::Real, mu_i::Real, zeff::Real, chi_perp::Real, chi_tor::Real, diff --git a/src/Tearing/Dispersion/Coupled.jl b/src/Tearing/Dispersion/Coupled.jl index 4029943b9..43a4b7b87 100644 --- a/src/Tearing/Dispersion/Coupled.jl +++ b/src/Tearing/Dispersion/Coupled.jl @@ -15,10 +15,10 @@ # det = mc(Q::ComplexF64) # # At each evaluation, for k = 1 .. msing_max, the inner-layer Δ is computed -# at a Q rescaled by `tauk_ref / tauk_k`, then subtracted (with the dc -# offset) from the diagonal of an `msing_max × msing_max` upper-left -# submatrix of `dp_matrix`. The off-diagonal Δ' couplings are passed -# through unchanged. +# at a Q rescaled by `tauk_k / tauk_ref` and offset by that surface's real +# Doppler shift `q_shift`, then subtracted (with the dc offset) from the +# diagonal of an `msing_max × msing_max` upper-left submatrix of `dp_matrix`. +# The off-diagonal Δ' couplings are passed through unchanged. """ MultiSurfaceCoupling{V<:AbstractVector{<:SurfaceCoupling}} @@ -29,11 +29,16 @@ normalization), and the truncation `msing_max` (number of surfaces actually participating in the determinant). Calling `mc(Q)` returns `det(M(Q))` where ``` -M[k,k] = dp_matrix[k,k] - scale_k · Δ_inner_k(Q · tauk_ref / tauk_k) - dc_k +M[k,k] = dp_matrix[k,k] - scale_k · Δ_inner_k(Q · ratio_k + q_shift_k) - dc_k M[i,j] = dp_matrix[i,j] for i ≠ j (off-diagonal Δ' couplings) ``` A root of `mc` in the complex `Q` plane is a coupled tearing eigenvalue. + +`ratio_k = tauk_k/tauk_ref`: `Q` is defined as `tauk·ω`, so one shared physical +ω reaches surface `k` as `Q_k = tauk_k·ω`, on the same scaling as its `Q_e`/`Q_i`. + +`q_shift_k` is the real offset carried on each `SurfaceCoupling` (zero by default). """ struct MultiSurfaceCoupling{V<:AbstractVector{<:SurfaceCoupling}} surfaces::V @@ -91,7 +96,7 @@ function (mc::MultiSurfaceCoupling)(Q::Number) M = mc.dp_matrix[1:n, 1:n] @inbounds for k in 1:n sc = mc.surfaces[k] - Q_k = Qc * (ref_tauk / sc.tauk) + Q_k = Qc * (sc.tauk / ref_tauk) + sc.q_shift # m×m scalar coupling: use only the tearing channel. The # interchange (Glasser-stabilization) channel is carried in the # full 4m×4m dispersion in `CoupledFullMatch.jl`; this reduced diff --git a/src/Tearing/Dispersion/CoupledFullMatch.jl b/src/Tearing/Dispersion/CoupledFullMatch.jl index 9b0f136c9..8d309a9d7 100644 --- a/src/Tearing/Dispersion/CoupledFullMatch.jl +++ b/src/Tearing/Dispersion/CoupledFullMatch.jl @@ -71,16 +71,12 @@ studies use the reduced m × m `MultiSurfaceCoupling` instead. - `dp_raw::Matrix{ComplexF64}` — 2m × 2m outer-region matrix (side-major). - `ref_idx::Int` — reference surface for Q rescaling (1-based). - `msing_max::Int` — number of surfaces to include (truncates). - - `rotation::Vector{Float64}` — per-surface rotation frequencies (s⁻¹). - - `ntor::Int` — toroidal mode number `n` (default 1). """ struct MultiSurfaceCouplingFull{V<:AbstractVector{<:SurfaceCoupling},K<:NamedTuple} surfaces::V dp_raw::Matrix{ComplexF64} ref_idx::Int msing_max::Int - rotation::Vector{Float64} - ntor::Int inner_kwargs::K # kwargs forwarded to solve_inner; e.g. (pfac=0.1, nx=128, nq=5) end @@ -88,16 +84,14 @@ end multi_surface_coupling_full(surfaces, dp_raw; ref_idx=1, msing_max=length(surfaces), - rotation=zeros(length(surfaces)), - ntor=1) -> MultiSurfaceCouplingFull + inner_kwargs=NamedTuple()) -> MultiSurfaceCouplingFull Construct the 4m × 4m dispersion matrix driver. `dp_raw` must be the 2m × 2m matrix in side-major ordering (the `intr.delta_prime_raw` field populated by `ForceFreeStates.compute_delta_prime_matrix!` on the -Riccati path). `rotation[k]` is the per-surface rotation -frequency; it shifts the per-surface inner Q argument by -`i·ntor·rotation[k]`. Default zero rotation matches the static-equilibrium -case. +Riccati path). Surface k's inner layer is evaluated at `Q·ratio_k + q_shift_k`, with +`ratio_k` and the real E×B Doppler offset `q_shift_k` exactly as in the reduced +`multi_surface_coupling`. # Keyword arguments @@ -107,9 +101,6 @@ case. matching matrix becomes 4·msing_max × 4·msing_max, built from the corresponding 2·msing_max × 2·msing_max submatrix of `dp_raw`. Defaults to `length(surfaces)`. - - `rotation` — per-surface rotation frequencies in s⁻¹ (length m). - Defaults to all zero. - - `ntor` — toroidal mode number n. Defaults to 1. - `inner_kwargs` — NamedTuple of kwargs forwarded to `solve_inner` at every Q evaluation, e.g. `(pfac=0.1, xfac=10.0, nx=128, nq=5)` for Galerkin grid tuning. Defaults to `NamedTuple()`. @@ -118,8 +109,6 @@ function multi_surface_coupling_full(surfaces::AbstractVector{<:SurfaceCoupling} dp_raw::AbstractMatrix; ref_idx::Integer=1, msing_max::Integer=length(surfaces), - rotation::AbstractVector{<:Real}=zeros(length(surfaces)), - ntor::Integer=1, inner_kwargs::NamedTuple=NamedTuple()) m = length(surfaces) size(dp_raw) == (2m, 2m) || @@ -131,14 +120,9 @@ function multi_surface_coupling_full(surfaces::AbstractVector{<:SurfaceCoupling} 1 <= msing_max <= m || throw(ArgumentError("multi_surface_coupling_full: msing_max=$msing_max " * "out of range 1:$m")) - length(rotation) == m || - throw(ArgumentError("multi_surface_coupling_full: rotation length " * - "$(length(rotation)) ≠ $m")) return MultiSurfaceCouplingFull(surfaces, Matrix{ComplexF64}(dp_raw), Int(ref_idx), Int(msing_max), - Float64.(collect(rotation)), - Int(ntor), inner_kwargs) end @@ -165,10 +149,9 @@ function (mc::MultiSurfaceCouplingFull)(Q::Number) idx3 = idx1 + s2 # d^k_+ idx4 = idx2 + s2 # d^k_- - # Per-surface Q shift: guess_modify = Q + i·n·rotation[k]. - # Also apply ref_tauk / sc.tauk rescaling (we keep the SurfaceCoupling - # tauk normalization that SLAYER needs; GGJ has tauk=1 so it's a no-op). - Q_k = Qc * (ref_tauk / sc.tauk) + 1im * mc.ntor * mc.rotation[k] + # Map the shared scanned Q onto this surface's normalization, then Doppler it into + # the surface's E×B frame (GGJ carries tauk = 1, so the ratio is a no-op there). + Q_k = Qc * (sc.tauk / ref_tauk) + sc.q_shift resp = solve_inner(sc.model, sc.params, Q_k; mc.inner_kwargs...) # delta1 = interchange (parity −), delta2 = tearing (parity +); named diff --git a/src/Tearing/Dispersion/GrowthRateExtraction.jl b/src/Tearing/Dispersion/GrowthRateExtraction.jl index 42c2c9681..45f5bc74d 100644 --- a/src/Tearing/Dispersion/GrowthRateExtraction.jl +++ b/src/Tearing/Dispersion/GrowthRateExtraction.jl @@ -22,6 +22,7 @@ using Contour using DelaunayTriangulation +using Random: Xoshiro # --------------------------------------------------------------------- # Public result struct + main entry point. @@ -955,7 +956,9 @@ function _extract_growth_rates_amr(Q::Vector{ComplexF64}, throw(ArgumentError("_extract_growth_rates_amr: need ≥ 3 points to triangulate")) pts = [(real(q), imag(q)) for q in Q] - tri = triangulate(pts) + # AMR grids are full of cocircular points, where the Delaunay triangulation is not unique and + # the randomized insertion order picks the diagonals; a fixed seed makes the contours reproducible. + tri = triangulate(pts; rng=Xoshiro(0)) # Segment types (carry complementary-field value at each endpoint) re_segs = NamedTuple{(:p1, :p2, :a1, :a2), diff --git a/src/Tearing/Dispersion/SurfaceCoupling.jl b/src/Tearing/Dispersion/SurfaceCoupling.jl index 0d57207ae..f243febec 100644 --- a/src/Tearing/Dispersion/SurfaceCoupling.jl +++ b/src/Tearing/Dispersion/SurfaceCoupling.jl @@ -21,13 +21,15 @@ """ SurfaceCoupling{M<:InnerLayerModel, P} -Per-surface dispersion data: `(model, params, dp_diag, dc, scale, tauk)`. -Calling `sc(Q)` returns the complex residual +Per-surface dispersion data: `(model, params, dp_diag, dc, scale, tauk, q_shift)`. Calling `sc(Q)` returns the complex residual ``` -r(Q) = dp_diag - scale * solve_inner(model, params, Q).tearing - dc +r(Q) = dp_diag - scale * solve_inner(model, params, Q + q_shift).tearing - dc ``` +`q_shift` is a real offset on the layer's Q argument, zero by default. The Tearing runner +sets it to the surface's E×B Doppler shift. + A root of `sc` in the complex `Q` plane is a **tearing** eigenvalue at this surface in the *uncoupled* approximation (only the tearing channel of the inner-layer response appears — the interchange channel enters the @@ -42,16 +44,18 @@ struct SurfaceCoupling{M<:InnerLayerModel,P} dc::Float64 scale::Float64 tauk::Float64 + q_shift::Float64 end function (sc::SurfaceCoupling)(Q::Number) - Δ = solve_inner(sc.model, sc.params, ComplexF64(Q)).tearing + Δ = solve_inner(sc.model, sc.params, ComplexF64(Q) + sc.q_shift).tearing return sc.dp_diag - sc.scale * Δ - sc.dc end """ surface_coupling(model::SLAYERModel, params::SLAYERParameters, - dp_diag::Number; dc::Real=0.0) -> SurfaceCoupling + dp_diag::Number; dc::Real=0.0, q_shift::Real=0.0) + -> SurfaceCoupling 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 @@ -63,9 +67,9 @@ 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) + dp_diag::Number; dc::Real=0.0, q_shift::Real=0.0) return SurfaceCoupling(model, params, ComplexF64(dp_diag), - Float64(dc), params.lu^(1 / 3), params.tauk) + Float64(dc), params.lu^(1 / 3), params.tauk, Float64(q_shift)) end """ @@ -83,25 +87,24 @@ interchange channel, which provides Glasser (Mercier) stabilization natively. A Δ_crit proxy (χ_parallel-matching offset on the diagonal) is meaningful only for tearing-only slab-layer approximations like SLAYER; for GGJ it would double-count the interchange physics. The `SurfaceCoupling` -struct's `dc` field is hard-wired to 0 here. +struct's `dc` field is hard-wired to 0 here, and so is `q_shift`. """ -function surface_coupling(model::GGJModel, params::GGJParameters, - dp_diag::Number) +function surface_coupling(model::GGJModel, params::GGJParameters, dp_diag::Number) return SurfaceCoupling(model, params, ComplexF64(dp_diag), - 0.0, 1.0, 1.0) + 0.0, 1.0, 1.0, 0.0) end """ surface_coupling(model::InnerLayerModel, params, dp_diag::Number; - dc::Real=0.0, scale::Real=1.0, tauk::Real=1.0) - -> SurfaceCoupling + dc::Real=0.0, scale::Real=1.0, tauk::Real=1.0, + q_shift::Real=0.0) -> SurfaceCoupling Generic fallback constructor. Use this when wiring a new inner-layer model into the dispersion solver — pass the appropriate inner→outer-units `scale` and per-surface `tauk` explicitly. """ function surface_coupling(model::InnerLayerModel, params, dp_diag::Number; - dc::Real=0.0, scale::Real=1.0, tauk::Real=1.0) + dc::Real=0.0, scale::Real=1.0, tauk::Real=1.0, q_shift::Real=0.0) return SurfaceCoupling(model, params, ComplexF64(dp_diag), - Float64(dc), Float64(scale), Float64(tauk)) + Float64(dc), Float64(scale), Float64(tauk), Float64(q_shift)) end diff --git a/src/Tearing/Runner/HDF5Output.jl b/src/Tearing/Runner/HDF5Output.jl index b0fc8d2e4..b5395bfde 100644 --- a/src/Tearing/Runner/HDF5Output.jl +++ b/src/Tearing/Runner/HDF5Output.jl @@ -44,6 +44,9 @@ function write_slayer_hdf5!(parent::Union{HDF5.File,HDF5.Group}, end _write_per_surface!(g, result.params, result.dp_matrix) + # E×B rotation resolved per surface, and the Doppler offset it applies to the inner-layer Q. + g["PerSurface/omega_E"] = result.omega_E + g["PerSurface/q_shift"] = result.q_shift # Surface identity (absent when the analysis was built from bare parameters). isempty(result.rational_psi) || (g["PerSurface/rational_psi"] = result.rational_psi) isempty(result.rational_q) || (g["PerSurface/rational_q"] = result.rational_q) @@ -82,6 +85,10 @@ const TEARING_H5_ANNOTATIONS = [ "PerSurface/iota_e" => (; long_name="electron fraction ι_e = Q_e/(Q_e − Q_i) per surface", dims=("surface",)), "PerSurface/tau_k" => (; long_name="Q-normalization time S^(1/3)·τ_H per surface (Q = τ_k·ω; diamagnetic inputs Q_e, Q_i carry the opposite sign by convention)", units="s", dims=("surface",)), + "PerSurface/omega_E" => + (; long_name="E×B angular frequency Ω_E per unit n on each surface, from the kinetic file", units="rad/s", dims=("surface",)), + "PerSurface/q_shift" => + (; 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",)), @@ -106,11 +113,21 @@ 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",)), + (; + long_name="mode rotation angular frequency ω = Re(Q)/τ_k of each root: each surface's E×B frame when uncoupled, the lab frame when coupled", + units="rad/s", + dims=("surface",) + ), "Roots/gamma" => (; long_name="growth rate γ = Im(Q)/τ_k of each root, positive = unstable (an e-folding rate, so no 2π distinction applies)", units="1/s", dims=("surface",)), "Roots/no_root" => (; long_name="flag: no usable dispersion root found (Q_root is NaN, ω/γ are placeholders)", dims=("surface",)), diff --git a/src/Tearing/Runner/Result.jl b/src/Tearing/Runner/Result.jl index 60a728869..3a018a38d 100644 --- a/src/Tearing/Runner/Result.jl +++ b/src/Tearing/Runner/Result.jl @@ -25,9 +25,10 @@ downstream inspection and HDF5 output. `delta_prime_to_rs_reference`, written as `PerSurface/Delta_prime_matrix_rs`); GGJ path: the ψ_N matrix unchanged (written as `PerSurface/Delta_prime_matrix`) - `Q_root` -- tearing eigenvalue(s) in normalized Q - * length `nsurfaces` in `:uncoupled` mode - * length `1` in `:coupled` mode (global eigenvalue normalized by - `params[1].tauk`) + + + length `nsurfaces` in `:uncoupled` mode + + length `1` in `:coupled` mode (global eigenvalue normalized by + `params[1].tauk`) - `omega_Hz`, `gamma_Hz` -- physical rotation frequency / growth rate - `per_surface_extraction` -- `Vector{GrowthRateResult}` of length `nsurfaces` in uncoupled mode (each includes polelines, pole list, @@ -37,6 +38,10 @@ downstream inspection and HDF5 output. - `layer_widths` -- `Vector{LayerWidths}`, one per surface: the resistive layer thickness (in meters) from the `del_s` Riccati solve plus FKR / visco-resistive sanity scales. Empty when disabled. + - `omega_E` -- E×B angular frequency Ω_E per unit n [rad/s] on each surface, + from the kinetic file. Recorded in both coupling modes + - `q_shift` -- real Doppler offset `−τ_k·n·Ω_E` applied to each surface's + inner-layer Q; zero outside `:coupled` mode, where rotation does not enter - `scan_data` -- scan results (per-surface in uncoupled, single entry in coupled). Empty unless `control.store_scan == true`. """ @@ -53,17 +58,20 @@ struct SLAYERResult per_surface_extraction::Vector{GrowthRateResult} coupled_extraction::Union{Nothing,GrowthRateResult} layer_widths::Vector{LayerWidths} + omega_E::Vector{Float64} + q_shift::Vector{Float64} scan_data::Vector{Union{ScanResult,AMRResult}} end # Empty result (enabled=false path) function empty_slayer_result(control::SLAYERControl) return SLAYERResult(false, control, - SLAYERParameters[], - Float64[], Float64[], - zeros(ComplexF64, 0, 0), - ComplexF64[], Float64[], Float64[], - GrowthRateResult[], nothing, - LayerWidths[], - Union{ScanResult,AMRResult}[]) + SLAYERParameters[], + Float64[], Float64[], + zeros(ComplexF64, 0, 0), + ComplexF64[], Float64[], Float64[], + GrowthRateResult[], nothing, + LayerWidths[], + Float64[], Float64[], + Union{ScanResult,AMRResult}[]) end diff --git a/src/Tearing/Runner/run_slayer.jl b/src/Tearing/Runner/run_slayer.jl index 2459087cf..0556feeb2 100644 --- a/src/Tearing/Runner/run_slayer.jl +++ b/src/Tearing/Runner/run_slayer.jl @@ -34,14 +34,12 @@ function _load_profiles(control::SLAYERControl, dir_path::AbstractString) error("run_slayer: kinetic file '$path' is missing required " * "dataset '$name' for the SLAYER inner layer.") end - # ω_*e/ω_*i are recomputed per-surface from equilibrium gradients - # (compute_omega_star), so the diamagnetic-frequency inputs here are - # placeholders; `omega` carries the ExB rotation when present. + # `omega_E` carries the E×B rotation Ω_E (per unit n) when the file has it; ω_*e/ω_*i are + # derived per surface from the density and temperature splines by `build_slayer_inputs`. npsi = length(data.psi) - omega = data.omega_E === nothing ? zeros(npsi) : data.omega_E + omega_E = data.omega_E === nothing ? zeros(npsi) : data.omega_E profiles = KineticProfiles(; psi=data.psi, n_e=data.n_e, T_e=data.T_e, - T_i=data.T_i, omega=omega, - omega_e=zeros(npsi), omega_i=zeros(npsi)) + T_i=data.T_i, omega_E=omega_E) # χ⊥(ψ)/χ_φ(ψ) splines from the file. A χ array that is absent OR all-zero # is treated as "not provided" — χ must be positive (χ=0 ⇒ τ_⊥→∞), and the @@ -139,17 +137,37 @@ end # --------------------------------------------------------------------- # SLAYER: scale = lu^(1/3), tauk from the surface, dc from the χ‖ proxy. function _build_surface_coupling(model::SLAYERModel, params::SLAYERParameters, - dp_diag) - return surface_coupling(model, params, dp_diag; dc=params.dc_tmp) + dp_diag, q_shift::Real) + return surface_coupling(model, params, dp_diag; dc=params.dc_tmp, q_shift=q_shift) end # GGJ: scale = 1.0 (rescale_delta applied inside solve_inner), tauk = 1.0, # dc = 0 (the 4m×4m Pletzer-Dewar residual carries interchange stabilization -# natively). See the GGJ `surface_coupling` method. -function _build_surface_coupling(model::GGJModel, params::GGJParameters, - dp_diag) - return surface_coupling(model, params, dp_diag) +# natively), and no E×B shift. See the GGJ `surface_coupling` method. +_build_surface_coupling(model::GGJModel, params::GGJParameters, dp_diag, ::Real) = + surface_coupling(model, params, dp_diag) + +""" + _q_shift(p::SLAYERParameters, Ω_E) -> Float64 + +Real offset added to a surface's scanned inner-layer `Q` for E×B rotation `Ω_E` per unit n [rad/s]. +The layer's `Q` is the mode frequency in the local E×B frame (Burgess et al. 2026, Table 1), while +the coupled scan runs in the lab frame, so + +``` +Q_layer = τ_k·ω_lab − τ_k·n·Ω_E ⇒ q_shift = −τ_k·n·Ω_E +``` + +The conjugated layer Δ and the opposite-sign diamagnetic inputs (`Q_e = −τ_k·ω_*e`) flip together +and leave this sign unchanged. +""" +function _q_shift(p::SLAYERParameters, Ω_E::Real) + iszero(Ω_E) || p.n >= 1 || + throw(ArgumentError("run_slayer: E×B rotation Ω_E=$(Ω_E) rad/s needs the toroidal mode number, but SLAYERParameters has n=$(p.n)")) + return -p.tauk * p.n * Ω_E end +# GGJ carries no time normalization (tauk = 1), so a physical rotation has no Q-space image. +_q_shift(::InnerLayerParameters, ::Real) = 0.0 # --------------------------------------------------------------------- # Reference-length conversion of the outer Δ' for the slab layer @@ -199,19 +217,23 @@ end """ run_slayer_from_inputs(params::Vector{SLAYERParameters}, dp_matrix::AbstractMatrix, - control::SLAYERControl) -> SLAYERResult + control::SLAYERControl; + rational_psi, rational_q, omega_E) -> SLAYERResult Run the SLAYER tearing analysis given pre-built per-surface `SLAYERParameters` and the outer-region Δ' matrix. Bypasses the equilibrium-driven `build_slayer_inputs` step — use this when the parameters are already known (e.g. in unit tests or when rebuilding -from cached HDF5 output). +from cached HDF5 output). `omega_E` is the E×B angular frequency per unit n +at each surface [rad/s]; empty means no rotation. In coupled mode it Doppler-shifts +each surface's inner-layer Q; in both modes it is reported as `result.omega_E`. """ function run_slayer_from_inputs(params::AbstractVector{<:InnerLayerParameters}, dp_matrix::AbstractMatrix, control::SLAYERControl; rational_psi::Vector{Float64}=Float64[], - rational_q::Vector{Float64}=Float64[]) + rational_q::Vector{Float64}=Float64[], + omega_E::Vector{Float64}=Float64[]) validate(control) control.enabled || return empty_slayer_result(control) isempty(params) && return empty_slayer_result(control) @@ -272,8 +294,28 @@ function run_slayer_from_inputs(params::AbstractVector{<:InnerLayerParameters}, ) end - # Per-surface SurfaceCoupling objects - scs = [_build_surface_coupling(model, params[k], dp[k, k]) for k in 1:n] + # E×B rotation only Doppler-shifts the coupled determinant: each uncoupled layer is already + # solved in its own plasma frame, so a shift there would only relabel the reported frequency. + coupled = control.coupling_mode === :coupled + isempty(omega_E) || length(omega_E) == n || + throw(ArgumentError("run_slayer: omega_E has $(length(omega_E)) entries but $n rational surfaces were analysed.")) + Ω_E = isempty(omega_E) ? zeros(Float64, n) : omega_E + coupled && _is_ggj(model) && any(!iszero, Ω_E) && @warn( + "SLAYER: E×B rotation is not applied to GGJ surfaces, which carry no time normalization.") + q_shifts = coupled ? Float64[_q_shift(params[k], Ω_E[k]) for k in 1:n] : zeros(Float64, n) + scs = [_build_surface_coupling(model, params[k], dp[k, k], q_shifts[k]) for k in 1:n] + # Coupled Q is scanned in the reference surface's normalization. + ref_idx = 1 + # In the scanned Q, surface k's layer response moves by τ_ref·n·Ω_E,k. Heuristic check only: + # the root also carries the diamagnetic Re(Q), so a window can pass and still miss it. + if any(!iszero, q_shifts) + m_use = min(control.msing_max, n) + doppler = [-q_shifts[k] * scs[ref_idx].tauk / scs[k].tauk for k in 1:m_use] + lo, hi = control.Q_re_range + any(d -> !(lo <= d <= hi), doppler) && @warn( + "SLAYER: E×B Doppler offsets Re(Q)=$(round.(doppler; digits=3)) fall outside " * + "Q_re_range=$(control.Q_re_range); widen it or the coupled root may be missed.") + end # Per-surface resistive layer thickness [m] via the del_s Riccati solve. # Independent of the dispersion scan / coupling mode — a pure diagnostic. @@ -335,16 +377,15 @@ function run_slayer_from_inputs(params::AbstractVector{<:InnerLayerParameters}, elseif control.coupling_mode === :coupled m_use = min(control.msing_max, n) - mc = multi_surface_coupling(scs, dp; ref_idx=1, msing_max=m_use) + mc = multi_surface_coupling(scs, dp; ref_idx=ref_idx, msing_max=m_use) scan = _run_scan(mc, control) pthr = _pole_threshold_for(scan) - ref_tauk = scs[1].tauk + ref_tauk = scs[ref_idx].tauk # Coupled path: no root polishing / validity gate. The m×m coupled # determinant det(D'−D(Q)) is ill-conditioned (its magnitude floors well # above zero), so |det|-based polishing is a no-op and the residual-scale - # gate is unreliable. A σ_min-based coupled refinement is a dev follow-up; - # for now the coupled determinant uses the raw contour extraction (the - # BLAS pin in the scan still makes it thread-deterministic). + # gate is unreliable. The coupled root is therefore the raw contour + # intersection and carries the triangulation's interpolation error. gr = find_growth_rates(scan, ref_tauk; pole_threshold=pthr, filter_above_poles=control.filter_above_poles, @@ -364,7 +405,7 @@ function run_slayer_from_inputs(params::AbstractVector{<:InnerLayerParameters}, return SLAYERResult(true, control, params, rational_psi, rational_q, dp, Q_root, omega_Hz, gamma_Hz, per_surface_extraction, coupled_extraction, - layer_widths, scan_data_list) + layer_widths, Ω_E, q_shifts, scan_data_list) end # --------------------------------------------------------------------- @@ -503,5 +544,7 @@ function run_slayer(equil, surfaces::AbstractVector, delta_prime_matrix::Abstrac rational_psi = Float64[surfaces[p.ising].psifac for p in params] rational_q = Float64[surfaces[p.ising].q for p in params] - return run_slayer_from_inputs(params, dp, control; rational_psi=rational_psi, rational_q=rational_q) + # E×B rotation per unit n at each surface, for the coupled determinant's Doppler shifts. + omega_E = Float64[profiles.omega_E(ψ) for ψ in rational_psi] + return run_slayer_from_inputs(params, dp, control; rational_psi=rational_psi, rational_q=rational_q, omega_E=omega_E) end diff --git a/src/Utilities/KineticProfiles.jl b/src/Utilities/KineticProfiles.jl index 70f1a2974..5b900b859 100644 --- a/src/Utilities/KineticProfiles.jl +++ b/src/Utilities/KineticProfiles.jl @@ -11,48 +11,43 @@ using FastInterpolations """ KineticProfiles -Radial kinetic-profile container. All six profiles are 1D cubic splines of +Radial kinetic-profile container. All four profiles are 1D cubic splines of the normalized poloidal flux ψ ∈ [0, 1]. -| field | meaning | units | -|:--------- |:--------------------------------------- |:----- | -| `n_e` | electron density | m⁻³ | -| `T_e` | electron temperature | eV | -| `T_i` | ion temperature | eV | -| `omega` | toroidal rotation | rad/s | -| `omega_e` | electron diamagnetic frequency ω\\_\\*e | rad/s | -| `omega_i` | ion diamagnetic frequency ω\\_\\*i | rad/s | +| field | meaning | units | +|:--------- |:---------------------------- |:----- | +| `n_e` | electron density | m⁻³ | +| `T_e` | electron temperature | eV | +| `T_i` | ion temperature | eV | +| `omega_E` | E×B rotation Ω_E, per unit n | rad/s | -Construct via the keyword constructor `KineticProfiles(; psi, n_e, T_e, T_i, omega, omega_e, omega_i)` with matched-length vectors. The SLAYER -runner builds this object from a standardized kinetic-profile file via -`Equilibrium.read_kinetic_file`. +Construct via the keyword constructor `KineticProfiles(; psi, n_e, T_e, T_i, omega_E)` with +matched-length vectors. The SLAYER runner builds this object from a standardized +kinetic-profile file via `Equilibrium.read_kinetic_file`. Diamagnetic frequencies are not +carried here: `build_slayer_inputs` derives ω_*e and ω_*i per surface from the density and +temperature splines. Evaluate all profiles at a given ψ via the call operator: ```julia -vals = kp(0.5) # NamedTuple(n_e=..., T_e=..., ..., omega_i=...) +vals = kp(0.5) # NamedTuple(n_e=..., T_e=..., T_i=..., omega_E=...) ``` """ struct KineticProfiles{S} n_e::S T_e::S T_i::S - omega::S - omega_e::S - omega_i::S + omega_E::S end function KineticProfiles(; psi::AbstractVector{<:Real}, n_e::AbstractVector{<:Real}, T_e::AbstractVector{<:Real}, T_i::AbstractVector{<:Real}, - omega::AbstractVector{<:Real}, - omega_e::AbstractVector{<:Real}, - omega_i::AbstractVector{<:Real}) + omega_E::AbstractVector{<:Real}) xs = collect(Float64.(psi)) for (name, v) in (("n_e", n_e), ("T_e", T_e), ("T_i", T_i), - ("omega", omega), ("omega_e", omega_e), - ("omega_i", omega_i)) + ("omega_E", omega_E)) length(v) == length(xs) || throw(ArgumentError("KineticProfiles: length($name) = $(length(v)) " * "≠ length(psi) = $(length(xs))")) @@ -60,22 +55,18 @@ function KineticProfiles(; psi::AbstractVector{<:Real}, return KineticProfiles(cubic_interp(xs, Float64.(n_e)), cubic_interp(xs, Float64.(T_e)), cubic_interp(xs, Float64.(T_i)), - cubic_interp(xs, Float64.(omega)), - cubic_interp(xs, Float64.(omega_e)), - cubic_interp(xs, Float64.(omega_i))) + cubic_interp(xs, Float64.(omega_E))) end """ (kp::KineticProfiles)(psi::Real) -> NamedTuple Evaluate all profiles at `psi` and return them as a NamedTuple with fields -`(n_e, T_e, T_i, omega, omega_e, omega_i)`. +`(n_e, T_e, T_i, omega_E)`. """ (kp::KineticProfiles)(psi::Real) = ( n_e=kp.n_e(psi), T_e=kp.T_e(psi), T_i=kp.T_i(psi), - omega=kp.omega(psi), - omega_e=kp.omega_e(psi), - omega_i=kp.omega_i(psi) + omega_E=kp.omega_E(psi) ) diff --git a/test/runtests.jl b/test/runtests.jl index e1256e56d..0079af5e4 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -52,6 +52,7 @@ else include("./runtests_dispersion_residual.jl") include("./runtests_dispersion_coupled.jl") include("./runtests_dispersion_coupled_full.jl") + include("./runtests_dispersion_rotation.jl") include("./runtests_dispersion_scan.jl") include("./runtests_dispersion_amr.jl") include("./runtests_dispersion_polish.jl") diff --git a/test/runtests_dispersion_amr.jl b/test/runtests_dispersion_amr.jl index 014f3d019..b2acd1065 100644 --- a/test/runtests_dispersion_amr.jl +++ b/test/runtests_dispersion_amr.jl @@ -13,7 +13,7 @@ # Small 2×2 initial grid → 9 unique corners amr = amr_scan(counting_f, (-1.0, 1.0), (-1.0, 1.0); - nre0=2, nim0=2, passes=0) + nre0=2, nim0=2, passes=0) @test amr isa AMRResult @test length(amr.cells) == 4 # 2×2 cells # Dedup: 9 unique corners (3×3) @@ -27,7 +27,7 @@ amr0 = amr_scan(f, (-1.0, 1.0), (-1.0, 1.0); nre0=4, nim0=4, passes=0) amr3 = amr_scan(f, (-1.0, 1.0), (-1.0, 1.0); nre0=4, nim0=4, passes=3) @test length(amr3.cells) > length(amr0.cells) - @test length(amr3.Q) > length(amr0.Q) + @test length(amr3.Q) > length(amr0.Q) # A 4×4 coarse grid is 16 cells; adding 3 refinement passes must # leave the total bounded by exponential growth of only the cells # bracketing the root (roughly linear in the path length). @@ -36,26 +36,26 @@ @testset "amr_scan: argument validation" begin @test_throws ArgumentError amr_scan(identity, (0.0, 1.0), (0.0, 1.0); - nre0=0, nim0=2, passes=1) + nre0=0, nim0=2, passes=1) @test_throws ArgumentError amr_scan(identity, (0.0, 1.0), (0.0, 1.0); - nre0=2, nim0=0, passes=1) + nre0=2, nim0=0, passes=1) @test_throws ArgumentError amr_scan(identity, (0.0, 1.0), (0.0, 1.0); - nre0=2, nim0=2, passes=-1) + nre0=2, nim0=2, passes=-1) end @testset "amr_scan: max_cells safety cap fires" begin # A pathological f that forces every cell to subdivide every pass f(Q) = 0.0 + 0.0im # identically zero → every cell crosses @test_throws ErrorException amr_scan(f, (-1.0, 1.0), (-1.0, 1.0); - nre0=4, nim0=4, passes=10, - max_cells=100) + nre0=4, nim0=4, passes=10, + max_cells=100) end @testset "find_growth_rates(AMR): single isolated root" begin Q_root = 0.42 + 0.27im f(Q) = ComplexF64(Q) - Q_root amr = amr_scan(f, (-1.0, 1.5), (-0.5, 1.0); - nre0=8, nim0=6, passes=4) + nre0=8, nim0=6, passes=4) result = find_growth_rates(amr, 1.0) @test result isa GrowthRateResult @test abs(result.Q_root - Q_root) < 1e-3 # AMR-resolution limited @@ -63,12 +63,25 @@ @test length(result.valid_roots) == 1 end + @testset "find_growth_rates(AMR): extraction does not draw from the global RNG" begin + # The AMR grid's cocircular points make the Delaunay triangulation non-unique, so the + # insertion order picks the diagonals and hence the contours; it must not come from the + # global RNG, which is seeded from entropy at startup. + using Random + amr = amr_scan(Q -> ComplexF64(Q) - (0.42 + 0.27im), (-1.0, 1.5), (-0.5, 1.0); nre0=8, nim0=6, passes=2) + Random.seed!(7) + find_growth_rates(amr, 1.0) + after = rand() + Random.seed!(7) + @test after == rand() + end + @testset "find_growth_rates(AMR): higher-γ root selected" begin Q1 = 0.3 + 0.5im # higher γ Q2 = -0.4 + 0.1im f(Q) = (ComplexF64(Q) - Q1) * (ComplexF64(Q) - Q2) amr = amr_scan(f, (-1.0, 1.0), (-0.3, 0.8); - nre0=10, nim0=8, passes=4) + nre0=10, nim0=8, passes=4) result = find_growth_rates(amr, 1.0) @test length(result.valid_roots) == 2 @test abs(result.Q_root - Q1) < 1e-2 @@ -79,7 +92,7 @@ Q_p = -0.5 + 0.6im f(Q) = (ComplexF64(Q) - Q_r) / (ComplexF64(Q) - Q_p) amr = amr_scan(f, (-1.5, 1.5), (-0.5, 1.5); - nre0=10, nim0=8, passes=5) + nre0=10, nim0=8, passes=5) result = find_growth_rates(amr, 1.0; pole_threshold=10.0) @test length(result.poles) >= 1 @test any(p -> abs(p - Q_p) < 0.05, result.poles) @@ -90,7 +103,7 @@ Q_root = 1.0 + 2.0im f(Q) = ComplexF64(Q) - Q_root amr = amr_scan(f, (-2.0, 3.0), (-1.0, 4.0); - nre0=8, nim0=8, passes=4) + nre0=8, nim0=8, passes=4) tauk = 5e-5 result = find_growth_rates(amr, tauk) @test result.omega_Hz ≈ real(result.Q_root) / tauk @@ -119,13 +132,13 @@ Q_root = 0.5 + 0.3im f(Q) = ComplexF64(Q) - Q_root scan = brute_force_scan(f, (-1.0, 1.0), (-0.5, 1.0); - nre=80, nim=60, threaded=false) - amr = amr_scan(f, (-1.0, 1.0), (-0.5, 1.0); - nre0=8, nim0=6, passes=4) + nre=80, nim=60, threaded=false) + amr = amr_scan(f, (-1.0, 1.0), (-0.5, 1.0); + nre0=8, nim0=6, passes=4) r_grid = find_growth_rates(scan, 1.0) - r_amr = find_growth_rates(amr, 1.0) + r_amr = find_growth_rates(amr, 1.0) @test abs(r_grid.Q_root - Q_root) < 1e-3 - @test abs(r_amr.Q_root - Q_root) < 1e-3 + @test abs(r_amr.Q_root - Q_root) < 1e-3 @test abs(r_grid.Q_root - r_amr.Q_root) < 5e-3 end @@ -139,23 +152,23 @@ InnerLayerResponse(m.a + m.b * ComplexF64(Q), zero(ComplexF64)) Q_pin = 0.7 - 0.3im - sc = surface_coupling(LinModel(0.0im, 1.0+0im), nothing, - Q_pin; scale=1.0, tauk=1.0) + sc = surface_coupling(LinModel(0.0im, 1.0 + 0im), nothing, + Q_pin; scale=1.0, tauk=1.0) amr = amr_scan(sc, (-0.5, 1.5), (-1.0, 0.5); - nre0=8, nim0=6, passes=4) + nre0=8, nim0=6, passes=4) r = find_growth_rates(amr, sc.tauk) @test abs(r.Q_root - Q_pin) < 1e-2 # Multi-surface coupled scan through AMR Q_a, Q_b = 0.7 - 0.3im, -0.4 + 0.5im - sc1 = surface_coupling(LinModel(0.0im, 1.0+0im), nothing, - ComplexF64(0); scale=1.0, tauk=1.0) - sc2 = surface_coupling(LinModel(0.0im, 1.0+0im), nothing, - ComplexF64(0); scale=1.0, tauk=1.0) + sc1 = surface_coupling(LinModel(0.0im, 1.0 + 0im), nothing, + ComplexF64(0); scale=1.0, tauk=1.0) + sc2 = surface_coupling(LinModel(0.0im, 1.0 + 0im), nothing, + ComplexF64(0); scale=1.0, tauk=1.0) dp = ComplexF64[Q_a 0.0; 0.0 Q_b] mc = multi_surface_coupling([sc1, sc2], dp) amr_c = amr_scan(mc, (-1.0, 1.5), (-1.0, 1.0); - nre0=10, nim0=8, passes=4) + nre0=10, nim0=8, passes=4) r_c = find_growth_rates(amr_c, mc.surfaces[mc.ref_idx].tauk) @test abs(r_c.Q_root - Q_b) < 1e-2 # higher-γ root end @@ -175,15 +188,15 @@ # otherwise crosses the entire real axis). f(Q) = (ComplexF64(Q) - 0.0) / (ComplexF64(Q) - (-50.0)) + (1.0 + 1.0im) boxes = [((-75.0, -25.0), (-25.0, 25.0)), - ((-25.0, 25.0), (-25.0, 25.0)), - (( 25.0, 75.0), (-25.0, 25.0))] + ((-25.0, 25.0), (-25.0, 25.0)), + ((25.0, 75.0), (-25.0, 25.0))] result = multi_box_amr_scan(f, boxes; - pole_magnitude_threshold=10.0, - prescreen_nre=25, prescreen_nim=25, - nre0=25, nim0=25, passes=2, - max_cells=100_000, - max_cells_action=:warn_truncate, - parallel=false) + pole_magnitude_threshold=10.0, + prescreen_nre=25, prescreen_nim=25, + nre0=25, nim0=25, passes=2, + max_cells=100_000, + max_cells_action=:warn_truncate, + parallel=false) @test result isa MultiBoxAMRResult @test length(result.box_results) == 3 @test length(result.box_activity) == 3 @@ -212,15 +225,15 @@ # independent of sign-change tests. g(Q) = 1000.0 / (ComplexF64(Q) - (-50.0))^2 + (5.0 + 5.0im) boxes = [((-75.0, -25.0), (-25.0, 25.0)), - ((-25.0, 25.0), (-25.0, 25.0)), - (( 25.0, 75.0), (-25.0, 25.0))] + ((-25.0, 25.0), (-25.0, 25.0)), + ((25.0, 75.0), (-25.0, 25.0))] result = multi_box_amr_scan(g, boxes; - pole_magnitude_threshold=50.0, - prescreen_nre=25, prescreen_nim=25, - nre0=25, nim0=25, passes=1, - max_cells=100_000, - max_cells_action=:warn_truncate, - parallel=false) + pole_magnitude_threshold=50.0, + prescreen_nre=25, prescreen_nim=25, + nre0=25, nim0=25, passes=1, + max_cells=100_000, + max_cells_action=:warn_truncate, + parallel=false) @test result.box_activity[1] != NoActivity @test result.box_activity[2] == NoActivity @test result.box_activity[3] == NoActivity diff --git a/test/runtests_dispersion_coupled.jl b/test/runtests_dispersion_coupled.jl index 1b950016e..2c40e5c0c 100644 --- a/test/runtests_dispersion_coupled.jl +++ b/test/runtests_dispersion_coupled.jl @@ -19,19 +19,19 @@ InnerLayerResponse(m.a + m.b * ComplexF64(Q), zero(ComplexF64)) function _slayer_ref() - return slayer_parameters( + return slayer_parameters(; n_e=5.0e19, t_e=1000.0, t_i=1000.0, - omega=0.0, omega_e=-1.0e4, omega_i=5.0e3, + omega_e=-1.0e4, omega_i=5.0e3, qval=2.0, sval_r=1.0, bt=2.0, rs=0.5, R0=1.7, mu_i=2.0, zeff=1.0, chi_perp=1.0, chi_tor=1.0, m=2, n=1) end @testset "Constructor validation" begin - sc1 = surface_coupling(LinTestModel(0.0im, 1.0+0im), nothing, - 1.0+0im; scale=1.0, tauk=1.0) - sc2 = surface_coupling(LinTestModel(0.0im, 1.0+0im), nothing, - 2.0+0im; scale=1.0, tauk=1.0) + sc1 = surface_coupling(LinTestModel(0.0im, 1.0 + 0im), nothing, + 1.0 + 0im; scale=1.0, tauk=1.0) + sc2 = surface_coupling(LinTestModel(0.0im, 1.0 + 0im), nothing, + 2.0 + 0im; scale=1.0, tauk=1.0) good_dp = ComplexF64[1.0 0.1; 0.1 2.0] mc = multi_surface_coupling([sc1, sc2], good_dp) @@ -40,25 +40,25 @@ @test size(mc.dp_matrix) == (2, 2) # 3-surface default also caps at 3 (min(3, 3) = 3) - sc3 = surface_coupling(LinTestModel(0.0im, 1.0+0im), nothing, - 3.0+0im; scale=1.0, tauk=1.0) + sc3 = surface_coupling(LinTestModel(0.0im, 1.0 + 0im), nothing, + 3.0 + 0im; scale=1.0, tauk=1.0) good_dp3 = ComplexF64[1.0 0.1 0.0; 0.1 2.0 0.0; 0.0 0.0 3.0] mc3 = multi_surface_coupling([sc1, sc2, sc3], good_dp3) @test mc3.msing_max == 3 # 4-surface case caps at 3 (the design default — Δ' beyond 3 surfaces # tends to be erratic in practice) - sc4 = surface_coupling(LinTestModel(0.0im, 1.0+0im), nothing, - 4.0+0im; scale=1.0, tauk=1.0) + sc4 = surface_coupling(LinTestModel(0.0im, 1.0 + 0im), nothing, + 4.0 + 0im; scale=1.0, tauk=1.0) good_dp4 = ComplexF64[1.0 0.0 0.0 0.0; - 0.0 2.0 0.0 0.0; - 0.0 0.0 3.0 0.0; - 0.0 0.0 0.0 4.0] + 0.0 2.0 0.0 0.0; + 0.0 0.0 3.0 0.0; + 0.0 0.0 0.0 4.0] mc4 = multi_surface_coupling([sc1, sc2, sc3, sc4], good_dp4) @test mc4.msing_max == 3 # default capped at 3 # Caller can opt in to all 4 mc4_full = multi_surface_coupling([sc1, sc2, sc3, sc4], good_dp4; - msing_max=4) + msing_max=4) @test mc4_full.msing_max == 4 # Mismatched dp size @@ -69,32 +69,32 @@ # Out-of-range ref_idx @test_throws ArgumentError multi_surface_coupling([sc1, sc2], good_dp; - ref_idx=3) + ref_idx=3) @test_throws ArgumentError multi_surface_coupling([sc1, sc2], good_dp; - ref_idx=0) + ref_idx=0) # Out-of-range msing_max @test_throws ArgumentError multi_surface_coupling([sc1, sc2], good_dp; - msing_max=3) + msing_max=3) @test_throws ArgumentError multi_surface_coupling([sc1, sc2], good_dp; - msing_max=0) + msing_max=0) end @testset "Diagonal Δ' factorizes (det = ∏ per-surface residuals)" begin # When dp_matrix is diagonal, no off-diagonal coupling exists and # the coupled determinant should reduce exactly to the product of # per-surface residuals. - sc1 = surface_coupling(LinTestModel(1.0+0im, 1.0+0im), nothing, - 5.0+0im; scale=1.0, tauk=1.0) - sc2 = surface_coupling(LinTestModel(2.0+0im, 1.0+0im), nothing, - 7.0+0im; scale=1.0, tauk=1.0) - sc3 = surface_coupling(LinTestModel(0.5+0im, 0.5+0im), nothing, - 3.0+0im; scale=1.0, tauk=1.0) + sc1 = surface_coupling(LinTestModel(1.0 + 0im, 1.0 + 0im), nothing, + 5.0 + 0im; scale=1.0, tauk=1.0) + sc2 = surface_coupling(LinTestModel(2.0 + 0im, 1.0 + 0im), nothing, + 7.0 + 0im; scale=1.0, tauk=1.0) + sc3 = surface_coupling(LinTestModel(0.5 + 0im, 0.5 + 0im), nothing, + 3.0 + 0im; scale=1.0, tauk=1.0) dp = ComplexF64[5.0 0.0 0.0; - 0.0 7.0 0.0; - 0.0 0.0 3.0] + 0.0 7.0 0.0; + 0.0 0.0 3.0] mc = multi_surface_coupling([sc1, sc2, sc3], dp) - for Q in (0.5+0im, 2.0+0.3im, -1.0-0.5im, 4.5+1.0im) + for Q in (0.5 + 0im, 2.0 + 0.3im, -1.0 - 0.5im, 4.5 + 1.0im) @test mc(Q) ≈ sc1(Q) * sc2(Q) * sc3(Q) rtol = 1e-12 end end @@ -103,23 +103,23 @@ # With Δ_inner(Q) = b·Q and dp_diag = b·Q_root for each surface, # the coupled determinant has its roots exactly at the union of # single-surface roots. - Q1, Q2 = 0.5+0.0im, 2.0+0.0im - sc1 = surface_coupling(LinTestModel(0.0im, 1.0+0im), nothing, - Q1; scale=1.0, tauk=1.0) - sc2 = surface_coupling(LinTestModel(0.0im, 1.0+0im), nothing, - Q2; scale=1.0, tauk=1.0) + Q1, Q2 = 0.5 + 0.0im, 2.0 + 0.0im + sc1 = surface_coupling(LinTestModel(0.0im, 1.0 + 0im), nothing, + Q1; scale=1.0, tauk=1.0) + sc2 = surface_coupling(LinTestModel(0.0im, 1.0 + 0im), nothing, + Q2; scale=1.0, tauk=1.0) dp = ComplexF64[real(Q1) 0.0; 0.0 real(Q2)] mc = multi_surface_coupling([sc1, sc2], dp) @test abs(mc(Q1)) < 1e-12 @test abs(mc(Q2)) < 1e-12 - @test abs(mc(0.0+0.0im)) > 0 + @test abs(mc(0.0 + 0.0im)) > 0 end @testset "Off-diagonal coupling shifts the roots away from the diagonal" begin - sc1 = surface_coupling(LinTestModel(0.0im, 1.0+0im), nothing, - 0.5+0im; scale=1.0, tauk=1.0) - sc2 = surface_coupling(LinTestModel(0.0im, 1.0+0im), nothing, - 2.0+0im; scale=1.0, tauk=1.0) + sc1 = surface_coupling(LinTestModel(0.0im, 1.0 + 0im), nothing, + 0.5 + 0im; scale=1.0, tauk=1.0) + sc2 = surface_coupling(LinTestModel(0.0im, 1.0 + 0im), nothing, + 2.0 + 0im; scale=1.0, tauk=1.0) # Coupling-free baseline dp_diag = ComplexF64[0.5 0.0; 0.0 2.0] mc_diag = multi_surface_coupling([sc1, sc2], dp_diag) @@ -137,61 +137,58 @@ end @testset "msing_max truncation uses upper-left submatrix" begin - sc1 = surface_coupling(LinTestModel(0.0im, 1.0+0im), nothing, - 1.0+0im; scale=1.0, tauk=1.0) - sc2 = surface_coupling(LinTestModel(0.0im, 1.0+0im), nothing, - 2.0+0im; scale=1.0, tauk=1.0) - sc3 = surface_coupling(LinTestModel(0.0im, 1.0+0im), nothing, - 3.0+0im; scale=1.0, tauk=1.0) + sc1 = surface_coupling(LinTestModel(0.0im, 1.0 + 0im), nothing, + 1.0 + 0im; scale=1.0, tauk=1.0) + sc2 = surface_coupling(LinTestModel(0.0im, 1.0 + 0im), nothing, + 2.0 + 0im; scale=1.0, tauk=1.0) + sc3 = surface_coupling(LinTestModel(0.0im, 1.0 + 0im), nothing, + 3.0 + 0im; scale=1.0, tauk=1.0) dp = ComplexF64[1.0 0.0 0.0; - 0.0 2.0 0.0; - 0.0 0.0 3.0] + 0.0 2.0 0.0; + 0.0 0.0 3.0] # msing_max = 1 reduces to sc1(Q) alone mc1 = multi_surface_coupling([sc1, sc2, sc3], dp; msing_max=1) - for Q in (0.0+0im, 1.0+0im, 2.0+0im) + for Q in (0.0 + 0im, 1.0 + 0im, 2.0 + 0im) @test mc1(Q) ≈ sc1(Q) end # msing_max = 2 uses the upper-left 2×2 → sc1·sc2 mc2 = multi_surface_coupling([sc1, sc2, sc3], dp; msing_max=2) - for Q in (0.0+0im, 0.5+0.5im) + for Q in (0.0 + 0im, 0.5 + 0.5im) @test mc2(Q) ≈ sc1(Q) * sc2(Q) end # msing_max = 3 (default for ≥3 surfaces) uses the full 3×3 → sc1·sc2·sc3 mc3 = multi_surface_coupling([sc1, sc2, sc3], dp) @test mc3.msing_max == 3 # min(3, 3) = 3 - for Q in (0.5+0.5im, 1.5-0.5im) + for Q in (0.5 + 0.5im, 1.5 - 0.5im) @test mc3(Q) ≈ sc1(Q) * sc2(Q) * sc3(Q) end end - @testset "Per-surface Q rescaling via tauk_ref / tauk_k" begin - # Each surface evaluates its inner Δ at Q_k = Q · (tauk_ref/tauk_k). - # With Δ(Q) = Q (b=1, a=0), the diagonal modification is - # M[k,k] = dp_diag_k - scale·Q·(tauk_ref/tauk_k) - # Verify against an explicit closed form with mismatched tauks. - sc1 = surface_coupling(LinTestModel(0.0im, 1.0+0im), nothing, - 0.0+0im; scale=1.0, tauk=2.0) # ref tauk - sc2 = surface_coupling(LinTestModel(0.0im, 1.0+0im), nothing, - 0.0+0im; scale=1.0, tauk=4.0) # half rate + @testset "Per-surface Q rescaling direction" begin + # Q is defined as tauk·ω, so one shared physical eigenvalue reaches + # surface k as Q_k = Q · (tauk_k/tauk_ref). + sc1 = surface_coupling(LinTestModel(0.0im, 1.0 + 0im), nothing, + 0.0 + 0im; scale=1.0, tauk=2.0) # ref tauk + sc2 = surface_coupling(LinTestModel(0.0im, 1.0 + 0im), nothing, + 0.0 + 0im; scale=1.0, tauk=4.0) dp = ComplexF64[0.0 0.0; 0.0 0.0] + mc = multi_surface_coupling([sc1, sc2], dp; ref_idx=1) - for Q in (1.0+0im, 0.5+0.3im) - # M[1,1] = 0 - Q · (2/2) = -Q - # M[2,2] = 0 - Q · (2/4) = -Q/2 - # det = M[1,1] · M[2,2] = Q·Q/2 = Q²/2 - @test mc(Q) ≈ Q^2 / 2 rtol = 1e-12 + for Q in (1.0 + 0im, 0.5 + 0.3im) + # M[1,1] = -Q·(2/2) = -Q ; M[2,2] = -Q·(4/2) = -2Q ; det = 2Q² + @test mc(Q) ≈ 2 * Q^2 rtol = 1e-12 end - # Switch ref_idx to surface 2 + # Reference surface is a labelling choice, not physics: the ratio + # tauk_k/tauk_ref rescales the scanned Q, so switching ref_idx + # rescales the whole determinant by a known power of the tauk ratio. mc2 = multi_surface_coupling([sc1, sc2], dp; ref_idx=2) - for Q in (1.0+0im, 0.5+0.3im) - # M[1,1] = -Q · (4/2) = -2Q - # M[2,2] = -Q · (4/4) = -Q - # det = 2Q · Q = 2Q² - @test mc2(Q) ≈ 2 * Q^2 rtol = 1e-12 + for Q in (1.0 + 0im, 0.5 + 0.3im) + # M[1,1] = -Q·(2/4) = -Q/2 ; M[2,2] = -Q·(4/4) = -Q ; det = Q²/2 + @test mc2(Q) ≈ Q^2 / 2 rtol = 1e-12 end end @@ -199,28 +196,35 @@ # Build a 2-surface SLAYER MultiSurfaceCoupling, evaluate at # Q_pin, and back-fill dp_matrix so that det(M(Q_pin)) = 0 # exactly. + # A second surface with different n_e and T_e, so tauk differs and the rescale is exercised. p_a = _slayer_ref() - p_b = _slayer_ref() + p_b = slayer_parameters(; + n_e=2.0e19, t_e=500.0, t_i=500.0, + omega_e=-1.0e4, omega_i=5.0e3, + qval=3.0, sval_r=1.5, bt=2.0, + rs=0.6, R0=1.7, mu_i=2.0, zeff=1.0, + chi_perp=1.0, chi_tor=1.0, m=3, n=1) m = SLAYERModel() - sc1 = surface_coupling(m, p_a, 0.0+0im) - sc2 = surface_coupling(m, p_b, 0.0+0im) + sc1 = surface_coupling(m, p_a, 0.0 + 0im) + sc2 = surface_coupling(m, p_b, 0.0 + 0im) + @test !(sc2.tauk ≈ sc1.tauk) Q_pin = 0.3 + 0.4im ref_tauk = sc1.tauk # Compute the diagonal modifications at Q_pin - Δ1 = solve_inner(m, p_a, Q_pin * (ref_tauk/sc1.tauk)).tearing * sc1.scale - Δ2 = solve_inner(m, p_b, Q_pin * (ref_tauk/sc2.tauk)).tearing * sc2.scale + Δ1 = solve_inner(m, p_a, Q_pin * (sc1.tauk / ref_tauk)).tearing * sc1.scale + Δ2 = solve_inner(m, p_b, Q_pin * (sc2.tauk / ref_tauk)).tearing * sc2.scale # Build dp such that M(Q_pin) is exactly singular. # Choose off-diagonal couplings, then set diagonals so M[k,k]=Δ_k # makes the matrix singular by setting M[1,1]·M[2,2] = M[1,2]·M[2,1]. - c12, c21 = 0.05+0im, 0.05+0im + c12, c21 = 0.05 + 0im, 0.05 + 0im # Pick M[1,1] arbitrarily, solve for M[2,2]: M11 = 0.7 + 0.0im M22 = (c12 * c21) / M11 - dp = ComplexF64[M11+Δ1 c12; - c21 M22+Δ2] + dp = ComplexF64[M11+Δ1 c12; + c21 M22+Δ2] mc = multi_surface_coupling([sc1, sc2], dp) # The constructed M(Q_pin) is exactly singular by construction @@ -232,8 +236,8 @@ @testset "GGJ surfaces flow through the coupled API" begin p = glasser_wang_2020_eq55() - sc1 = surface_coupling(GGJModel(solver=:shooting), p, -1.0+0im) - sc2 = surface_coupling(GGJModel(solver=:shooting), p, -2.0+0im) + sc1 = surface_coupling(GGJModel(; solver=:shooting), p, -1.0 + 0im) + sc2 = surface_coupling(GGJModel(; solver=:shooting), p, -2.0 + 0im) dp = ComplexF64[-1.0 0.1; 0.1 -2.0] mc = multi_surface_coupling([sc1, sc2], dp) @test mc isa MultiSurfaceCoupling @@ -243,14 +247,14 @@ @testset "Broadcast over a 2D Q grid" begin # Coupled residual must be broadcast-compatible for PR 5/6 scans. - sc1 = surface_coupling(LinTestModel(0.0im, 1.0+0im), nothing, - 0.0+0im; scale=1.0, tauk=1.0) - sc2 = surface_coupling(LinTestModel(0.0im, 1.0+0im), nothing, - 0.0+0im; scale=1.0, tauk=1.0) + sc1 = surface_coupling(LinTestModel(0.0im, 1.0 + 0im), nothing, + 0.0 + 0im; scale=1.0, tauk=1.0) + sc2 = surface_coupling(LinTestModel(0.0im, 1.0 + 0im), nothing, + 0.0 + 0im; scale=1.0, tauk=1.0) dp = ComplexF64[0.0 0.0; 0.0 0.0] mc = multi_surface_coupling([sc1, sc2], dp) - Q_grid = [(qr + qi*im) for qr in -1.0:0.5:1.0, qi in -1.0:0.5:1.0] + Q_grid = [(qr + qi * im) for qr in -1.0:0.5:1.0, qi in -1.0:0.5:1.0] det_grid = mc.(Q_grid) @test size(det_grid) == size(Q_grid) @test all(d -> d isa ComplexF64, det_grid) diff --git a/test/runtests_dispersion_coupled_full.jl b/test/runtests_dispersion_coupled_full.jl index d471650f8..3b4143f5b 100644 --- a/test/runtests_dispersion_coupled_full.jl +++ b/test/runtests_dispersion_coupled_full.jl @@ -31,8 +31,6 @@ @test size(mc.dp_raw) == (4, 4) @test mc.msing_max == 2 @test mc.ref_idx == 1 - @test mc.rotation == [0.0, 0.0] - @test mc.ntor == 1 # Wrong outer dim @test_throws ArgumentError multi_surface_coupling_full([sc1, sc2], @@ -45,9 +43,6 @@ dp_raw; msing_max=0) @test_throws ArgumentError multi_surface_coupling_full([sc1, sc2], dp_raw; msing_max=3) - # Wrong rotation length - @test_throws ArgumentError multi_surface_coupling_full([sc1, sc2], - dp_raw; rotation=[0.0]) end @testset "1-surface 4×4 det matches hand computation" begin @@ -120,20 +115,9 @@ @test det_jl ≈ det(M) atol = 1e-12 * abs(det(M)) end - @testset "Rotation shift applies i·ntor·rotation to inner Q argument" begin - # Ensure the per-surface rotation enters the inner-layer argument. - # Use a linear Δ_t model so Q-dependence is tractable. - dp_raw = ComplexF64[1.0 0; 0 1.0] - # Δ_t(Q) = Q (pure linear), Δ_i(Q) = 0 - sc = surface_coupling(_LinearInnerF(0 + 0im, 1 + 0im, 0 + 0im, 0 + 0im), - nothing, 0 + 0im; scale=1.0, tauk=1.0, dc=0.0) - # Case A: rotation=0, Q=2+0im → inner sees 2+0im → Δ_t=2, Δ_i=0 - mc0 = multi_surface_coupling_full([sc], dp_raw; rotation=[0.0], ntor=1) - # Case B: rotation=3, Q=2+0im → inner sees 2 + 1j*1*3 = 2+3i → Δ_t=2+3i - mcR = multi_surface_coupling_full([sc], dp_raw; rotation=[3.0], ntor=1) - @test mc0(2.0 + 0.0im) ≠ mcR(2.0 + 0.0im) - - # Check by hand. Both with the same outer matrix: + @testset "Each surface sees Q·tauk_k/tauk_ref plus its real q_shift" begin + # Linear Δ_t(Q) = Q, Δ_i = 0. With block-diagonal dp_raw the 8×8 matching matrix is a + # symmetric permutation of two per-surface 4×4 blocks, so det factorizes. function detAt(Δ_t, Δ_i) M = ComplexF64[ 1 0 -1 1; @@ -142,8 +126,15 @@ 0 1 -Δ_i -Δ_t] return det(M) end - @test mc0(2.0 + 0.0im) ≈ detAt(2.0 + 0.0im, 0.0 + 0.0im) - @test mcR(2.0 + 0.0im) ≈ detAt(2.0 + 3.0im, 0.0 + 0.0im) + lin = _LinearInnerF(0 + 0im, 1 + 0im, 0 + 0im, 0 + 0im) + s1, s2 = 0.25, -0.75 + sc1 = surface_coupling(lin, nothing, 0 + 0im; scale=1.0, tauk=2.0, dc=0.0, q_shift=s1) + sc2 = surface_coupling(lin, nothing, 0 + 0im; scale=1.0, tauk=4.0, dc=0.0, q_shift=s2) + dp_raw = Matrix{ComplexF64}(I, 4, 4) + Q = 1.5 + 0.5im + + mc = multi_surface_coupling_full([sc1, sc2], dp_raw) + @test mc(Q) ≈ detAt(Q * (2.0 / 2.0) + s1, 0im) * detAt(Q * (4.0 / 2.0) + s2, 0im) end @testset "SurfaceCoupling scale multiplies both inner channels" begin diff --git a/test/runtests_dispersion_residual.jl b/test/runtests_dispersion_residual.jl index e8af06da4..9f283b071 100644 --- a/test/runtests_dispersion_residual.jl +++ b/test/runtests_dispersion_residual.jl @@ -21,7 +21,7 @@ function _slayer_ref() return slayer_parameters( n_e=5.0e19, t_e=1000.0, t_i=1000.0, - omega=0.0, omega_e=-1.0e4, omega_i=5.0e3, + omega_e=-1.0e4, omega_i=5.0e3, qval=2.0, sval_r=1.0, bt=2.0, rs=0.5, R0=1.7, mu_i=2.0, zeff=1.0, chi_perp=1.0, chi_tor=1.0, m=2, n=1) diff --git a/test/runtests_dispersion_rotation.jl b/test/runtests_dispersion_rotation.jl new file mode 100644 index 000000000..488fa95b8 --- /dev/null +++ b/test/runtests_dispersion_rotation.jl @@ -0,0 +1,87 @@ +@testset "Dispersion per-surface rotation shift" begin + using GeneralizedPerturbedEquilibrium.InnerLayer + using GeneralizedPerturbedEquilibrium.InnerLayer: InnerLayerModel, solve_inner + using GeneralizedPerturbedEquilibrium.Dispersion + using GeneralizedPerturbedEquilibrium.Tearing.Runner: _q_shift, _build_surface_coupling + using LinearAlgebra + + # Linear inner layer Δ(Q) = a + b·Q makes the applied Q offset readable + # straight off the residual. + struct RotTestModel <: InnerLayerModel + a::ComplexF64 + b::ComplexF64 + end + GeneralizedPerturbedEquilibrium.InnerLayer.solve_inner( + m::RotTestModel, params, Q::Number) = + InnerLayerResponse(m.a + m.b * ComplexF64(Q), zero(ComplexF64)) + + model = RotTestModel(0.0im, 1.0 + 0im) + + @testset "q_shift offsets the layer Q argument on the scalar residual" begin + shift = 0.75 + sc0 = surface_coupling(model, nothing, 1.0 + 0im; scale=1.0, tauk=1.0) + scs = surface_coupling(model, nothing, 1.0 + 0im; scale=1.0, tauk=1.0, + q_shift=shift) + @test scs.q_shift == shift + for Q in (0.0 + 0im, 2.0 - 1.0im, -3.5 + 0.25im) + @test scs(Q) ≈ sc0(Q + shift) + end + end + + @testset "Coupled determinant applies each surface's own shift" begin + # Diagonal Δ', so det = Π_k (dp_kk - Δ_k(Q·tauk_k/tauk_ref + shift_k)). + s1, s2 = 0.5, -1.25 + sc1 = surface_coupling(model, nothing, 1.0 + 0im; scale=1.0, tauk=1.0, + q_shift=s1) + sc2 = surface_coupling(model, nothing, 2.0 + 0im; scale=1.0, tauk=2.0, + q_shift=s2) + dp = ComplexF64[1.0 0.0; 0.0 2.0] + mc = multi_surface_coupling([sc1, sc2], dp) + + Q = 1.5 + 0.5im + expected = (1.0 - (Q * (1.0 / 1.0) + s1)) * (2.0 - (Q * (2.0 / 1.0) + s2)) + @test mc(Q) ≈ expected + end + + @testset "Real SLAYER surfaces: the lab-frame root moves with the E×B rotation" begin + # Q_layer = τ_k·ω_lab − τ_k·n·Ω_E, so the rotating determinant at Q + τ_ref·n·Ω_E equals the + # static one at Q: every root keeps its γ and moves by ω_lab(Ω_E) − ω_lab(0) = +n·Ω_E. + mk(; qval, rs, m) = slayer_parameters(; n_e=5.0e19, t_e=1000.0, t_i=1000.0, + omega_e=1.0e4, omega_i=5.0e3, + qval=qval, sval_r=1.0, bt=2.0, rs=rs, R0=1.7, + mu_i=2.0, zeff=1.0, chi_perp=1.0, chi_tor=1.0, m=m, n=2) + params = [mk(; qval=1.5, rs=0.5, m=3), mk(; qval=2.0, rs=0.6, m=4)] + slayer = SLAYERModel(; variant=:fitzpatrick) + dp = ComplexF64[-2.0 0.5; 0.5 -3.0] + lab_shift = 0.5 # τ_ref·n·Ω_E in the reference Q + Ω_rigid = lab_shift / (params[1].tauk * params[1].n) # rad/s per unit n, on both surfaces + build(Ω) = [_build_surface_coupling(slayer, params[k], dp[k, k], _q_shift(params[k], Ω[k])) for k in 1:2] + static = multi_surface_coupling(build(zeros(2)), dp; ref_idx=1, msing_max=2) + rotating = multi_surface_coupling(build(fill(Ω_rigid, 2)), dp; ref_idx=1, msing_max=2) + for Q in (0.3 + 0.2im, -0.4 + 0.6im, 0.1 + 1.1im) + # The layer arguments agree to rounding, which the adaptive layer ODE amplifies to ≲ 3e-6. + @test rotating(Q + lab_shift) ≈ static(Q) rtol = 1e-4 + # The opposite Doppler sign (mode counter-rotating with the plasma) must not match. + @test !isapprox(rotating(Q - lab_shift), static(Q); rtol=1e-3) + end + end + + @testset "Coupled roots do not depend on the reference surface" begin + # Diagonal Δ′ decouples the surfaces, so each coupled root must be that surface's own + # root in physical units whichever surface normalizes Q. + d = ComplexF64[1.0+1.0im, 3.0+2.0im] + tauk = [2.0e-4, 5.0e-4] + scs = [surface_coupling(model, nothing, d[k]; scale=1.0, tauk=tauk[k]) for k in 1:2] + # Δ(Q) = Q makes surface k's own root Q_k = d_k, i.e. ω + iγ = d_k/τ_k. + expected = sort(d ./ tauk; by=real) + for ref in 1:2 + mc = multi_surface_coupling(scs, diagm(d); ref_idx=ref, msing_max=2) + # det is quadratic in Q: recover it from three samples and solve. + xs = ComplexF64[0, 1, 2] + c = [x^j for x in xs, j in 0:2] \ [mc(x) for x in xs] + disc = sqrt(c[2]^2 - 4c[3] * c[1]) + roots = [(-c[2] + disc) / (2c[3]), (-c[2] - disc) / (2c[3])] + @test sort(roots ./ tauk[ref]; by=real) ≈ expected rtol = 1e-10 + end + end +end diff --git a/test/runtests_kinetic_profiles.jl b/test/runtests_kinetic_profiles.jl index d7308b2d1..89e364945 100644 --- a/test/runtests_kinetic_profiles.jl +++ b/test/runtests_kinetic_profiles.jl @@ -10,9 +10,7 @@ "n_e" => fill(5.0e19, length(psi)), "T_e" => 1000.0 .* (1.0 .- 0.7 .* psi), "T_i" => 1200.0 .* (1.0 .- 0.6 .* psi), - "omega" => 1.0e4 .* psi, - "omega_e" => fill(1.0e4, length(psi)), - "omega_i" => fill(5.0e3, length(psi)) + "omega_E" => 1.0e4 .* psi ) ) end @@ -20,23 +18,20 @@ @testset "kwarg constructor + evaluation" begin psi, d = _synthetic() kp = KineticProfiles(; psi=psi, n_e=d["n_e"], T_e=d["T_e"], - T_i=d["T_i"], omega=d["omega"], - omega_e=d["omega_e"], omega_i=d["omega_i"]) + T_i=d["T_i"], omega_E=d["omega_E"]) # Exact recovery at a node vals = kp(0.5) @test vals.n_e ≈ 5.0e19 @test vals.T_e ≈ 1000.0 * (1 - 0.7 * 0.5) @test vals.T_i ≈ 1200.0 * (1 - 0.6 * 0.5) - @test vals.omega ≈ 1.0e4 * 0.5 - @test vals.omega_e ≈ 1.0e4 - @test vals.omega_i ≈ 5.0e3 + @test vals.omega_E ≈ 1.0e4 * 0.5 # Smooth interpolation between nodes vals2 = kp(0.25) @test vals2.T_e ≈ 1000.0 * (1 - 0.7 * 0.25) rtol = 1e-6 # NamedTuple fields - @test keys(vals) == (:n_e, :T_e, :T_i, :omega, :omega_e, :omega_i) + @test keys(vals) == (:n_e, :T_e, :T_i, :omega_E) end @testset "length mismatch raises" begin @@ -46,8 +41,6 @@ n_e=fill(1.0, length(psi) - 1), # wrong length T_e=fill(1000.0, length(psi)), T_i=fill(1000.0, length(psi)), - omega=fill(0.0, length(psi)), - omega_e=fill(0.0, length(psi)), - omega_i=fill(0.0, length(psi))) + omega_E=fill(0.0, length(psi))) end end diff --git a/test/runtests_resist_eval.jl b/test/runtests_resist_eval.jl index 8e334a67a..8ce0c4898 100644 --- a/test/runtests_resist_eval.jl +++ b/test/runtests_resist_eval.jl @@ -105,9 +105,7 @@ n_e=fill(5.0e19, length(psi_pts)), T_e=1000.0 .* (1.0 .- 0.7 .* psi_pts), T_i=1000.0 .* (1.0 .- 0.6 .* psi_pts), - omega=fill(0.0, length(psi_pts)), - omega_e=fill(1.0e4, length(psi_pts)), - omega_i=fill(5.0e3, length(psi_pts))) + omega_E=fill(0.0, length(psi_pts))) dq = deriv_view(equil.profiles.q_spline, 1) s1 = SingType(; psifac=0.3, rho=sqrt(0.3), m=[2], n=[1], @@ -149,7 +147,7 @@ n = length(psi_pts) profiles = KineticProfiles(; psi=psi_pts, n_e=fill(5.0e19, n), T_e=fill(1000.0, n), T_i=fill(1000.0, n), - omega=fill(0.0, n), omega_e=fill(1.0e4, n), omega_i=fill(5.0e3, n)) + omega_E=fill(0.0, n)) s_unpop = SingType(; psifac=0.5, rho=sqrt(0.5), m=[2], n=[1], q=2.0, q1=1.0, @@ -168,9 +166,7 @@ n_e=fill(5.0e19, length(psi_pts)), T_e=1000.0 .* (1.0 .- 0.7 .* psi_pts), T_i=fill(1000.0, length(psi_pts)), - omega=fill(0.0, length(psi_pts)), - omega_e=fill(0.0, length(psi_pts)), - omega_i=fill(0.0, length(psi_pts))) + omega_E=fill(0.0, length(psi_pts))) dq = deriv_view(equil.profiles.q_spline, 1) s1 = SingType(; psifac=0.3, rho=sqrt(0.3), m=[2], n=[1], diff --git a/test/runtests_slayer_inputs.jl b/test/runtests_slayer_inputs.jl index fce001036..fd1c10827 100644 --- a/test/runtests_slayer_inputs.jl +++ b/test/runtests_slayer_inputs.jl @@ -20,12 +20,10 @@ n_e=fill(5.0e19, length(psi_pts)), T_e=1000.0 .* (1.0 .- 0.7 .* psi_pts), T_i=1000.0 .* (1.0 .- 0.6 .* psi_pts), - omega=fill(0.0, length(psi_pts)), - omega_e=fill(1.0e4, length(psi_pts)), - omega_i=fill(5.0e3, length(psi_pts))) + omega_E=fill(0.0, 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,16 +70,12 @@ 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. - # compute_omega_star=false makes Q_e/Q_i pass through directly from profiles.omega_e/i - # rather than being recomputed from n_e/T_e/T_i gradients — required for the Q_e == - # -tauk·omega_e(ψ) identity check below. - sl = build_slayer_inputs(equil, sings, profiles; bt=2.0, dr_val=0.0, - compute_omega_star=false) + sl = build_slayer_inputs(equil, sings, profiles; bt=2.0, dr_val=0.0) @test length(sl) == 2 @test sl[1] isa SLAYERParameters @@ -109,9 +103,14 @@ @test sl[1].lu != sl[2].lu @test sl[1].tauk != sl[2].tauk - # Q_e, Q_i follow the SLAYER layerinputs sign convention - @test sl[1].Q_e == -sl[1].tauk * profiles.omega_e(0.3) - @test sl[1].Q_i == -sl[1].tauk * profiles.omega_i(0.3) + # Q_e, Q_i follow the SLAYER layerinputs sign convention, built on the ω_* derived from + # the profile splines. These fixtures are linear in ψ with constant n_e, so + # ω_*e = (2π/χ₁)·dT_e/dψ = -700/psio and ω_*i = -(2π/χ₁)·dT_i/dψ = +600/psio exactly. + ω_star_e = -700.0 / equil.psio + ω_star_i = 600.0 / equil.psio + @test sl[1].Q_e ≈ -sl[1].tauk * ω_star_e rtol = 1e-10 + @test sl[1].Q_i ≈ -sl[1].tauk * ω_star_i rtol = 1e-10 + @test sl[1].Q_e > 0 > sl[1].Q_i end @testset "build_slayer_inputs: bt defaults to the physical B_T, not a normalization" begin @@ -126,12 +125,12 @@ # 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) - sl_at_phys = build_slayer_inputs(equil, sings, profiles; bt=bt_phys, dr_val=0.0, compute_omega_star=false) - sl_double = build_slayer_inputs(equil, sings, profiles; bt=2 * bt_phys, dr_val=0.0, compute_omega_star=false) + sl_default = build_slayer_inputs(equil, sings, profiles; dr_val=0.0) + sl_at_phys = build_slayer_inputs(equil, sings, profiles; bt=bt_phys, dr_val=0.0) + sl_double = build_slayer_inputs(equil, sings, profiles; bt=2 * bt_phys, dr_val=0.0) # Leaving bt unset must be identical to passing the physical field explicitly. @test sl_default[1].lu ≈ sl_at_phys[1].lu rtol = 1e-12 @@ -143,14 +142,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 +166,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 +185,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; @@ -216,8 +215,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..8822e954b 100644 --- a/test/runtests_slayer_params.jl +++ b/test/runtests_slayer_params.jl @@ -9,7 +9,7 @@ function _ref_kwargs(; dr_val=0.0, dc_type=:none) return ( n_e=5.0e19, t_e=1000.0, t_i=1000.0, - omega=0.0, omega_e=-1.0e4, omega_i=5.0e3, + omega_e=-1.0e4, omega_i=5.0e3, qval=2.0, sval_r=1.0, bt=2.0, rs=0.5, R0=1.7, mu_i=2.0, zeff=1.0, chi_perp=1.0, chi_tor=1.0, diff --git a/test/runtests_slayer_riccati.jl b/test/runtests_slayer_riccati.jl index 9d4a95462..48a6263b2 100644 --- a/test/runtests_slayer_riccati.jl +++ b/test/runtests_slayer_riccati.jl @@ -15,7 +15,7 @@ function _ref_params_large_D() return slayer_parameters(; n_e=5.0e19, t_e=3000.0, t_i=3000.0, - omega=0.0, omega_e=-1.0e4, omega_i=5.0e3, + omega_e=-1.0e4, omega_i=5.0e3, qval=2.0, sval_r=1.0, bt=2.0, rs=0.5, R0=1.7, mu_i=2.0, zeff=1.0, chi_perp=1.0, chi_tor=1.0, diff --git a/test/runtests_slayer_runner.jl b/test/runtests_slayer_runner.jl index ee114e37f..de18d5623 100644 --- a/test/runtests_slayer_runner.jl +++ b/test/runtests_slayer_runner.jl @@ -136,6 +136,42 @@ end end + @testset "run_slayer_from_inputs: E×B rotation enters only the coupled determinant" begin + params = [_mk_params(; rs=0.5, lu=1.0e7, tauk=1.0e-4, m=2, n=1, ising=1), + _mk_params(; rs=0.6, lu=2.0e7, tauk=1.2e-4, m=3, n=1, ising=2)] + dpm = ComplexF64[-2.0 0.0; 0.0 -3.0] + Ω_E = [2.0e3, -1.0e3] + grid = (; enabled=true, scan_mode=:brute_force, Q_re_range=(-1.0, 1.0), Q_im_range=(-0.5, 0.8), + nre=8, nim=8, pole_threshold=1e5) + + # Uncoupled layers are solved in their own plasma frame: no shift, even when supplied. + r_unc = run_slayer_from_inputs(params, dpm, SLAYERControl(; coupling_mode=:uncoupled, grid...); omega_E=Ω_E) + @test r_unc.q_shift == [0.0, 0.0] + # The resolved rotation is still reported in both modes. + @test r_unc.omega_E == Ω_E + + # Coupled: the kinetic-file Ω_E Dopplers each surface by −τ_k·n·Ω_E. + r_cpl = run_slayer_from_inputs(params, dpm, SLAYERControl(; coupling_mode=:coupled, grid...); omega_E=Ω_E) + @test r_cpl.q_shift ≈ [-1.0e-4 * 1 * Ω_E[1], -1.2e-4 * 1 * Ω_E[2]] + @test r_cpl.omega_E == Ω_E + # The mode sees n·Ω_E: the same rotation shifts an n = 2 surface twice as far. + r_n2 = run_slayer_from_inputs([_mk_params(; tauk=1.0e-4, m=3, n=2)], ComplexF64[-2.0;;], + SLAYERControl(; coupling_mode=:coupled, grid...); omega_E=[2.0e3]) + @test r_n2.q_shift ≈ [-1.0e-4 * 2 * 2.0e3] + @test_throws "rational surfaces were analysed" run_slayer_from_inputs(params, dpm, + SLAYERControl(; coupling_mode=:coupled, grid...); omega_E=[1.0]) + + # A Doppler offset beyond the scan box is flagged rather than silently losing the root: + # τ_ref·n·Ω_E = 1e-4 · 3e4 = 3 lies outside Re(Q) ∈ [-1, 1]. + @test_logs (:warn, r"outside Q_re_range") match_mode = :any run_slayer_from_inputs(params, dpm, + SLAYERControl(; coupling_mode=:coupled, grid...); omega_E=[3.0e4, 0.0]) + + # Rotation without a toroidal mode number would silently vanish, so it is an error. + params_n0 = [_mk_params(; rs=0.5, lu=1.0e7, tauk=1.0e-4, m=2, n=0, ising=1)] + @test_throws "needs the toroidal mode number" run_slayer_from_inputs(params_n0, ComplexF64[-2.0;;], + SLAYERControl(; coupling_mode=:coupled, grid...); omega_E=[2.0e3]) + end + @testset "run_slayer_from_inputs: disabled path is a no-op" begin c = SLAYERControl(; enabled=false) params = [_mk_params()] @@ -182,7 +218,7 @@ 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, alpha_mercier=0.5) - dp1 = ComplexF64[5.0+0im;;] + 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