InnerLayer.SLAYER - BUGFIX - Restore the toroidal mode number in the diamagnetic frequencies

.gitignore

@@ -38,3 +38,6 @@ scratch/
 
 # Local profiling scratch (one-off study scripts, not part of the package)
 profiling/
+
+# Override-case rundirs materialized beside their example by the regression harness
+.regress-override-*

regression-harness/cases/diiid_slayer_n1.toml

@@ -1,6 +1,8 @@
 # Regression case: DIII-D-like H-mode n=1 SLAYER tearing-mode growth rates.
 # Uncoupled per-surface analysis with the AMR Q-plane scan; tracks the per-surface
-# layer parameters and the extracted roots (Q, ω, γ) for the inner three rationals.
+# layer parameters, the diamagnetic drive Q_e/Q_i, and the extracted roots (Q, ω, γ)
+# for the inner three rationals. The n-scaling of Q_e/Q_i, which no single-n case can
+# constrain, is covered by a unit test in test/runtests_slayer_inputs.jl instead.
 # 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]
@@ -89,6 +91,24 @@ label = "SLAYER iota_e"
 noise_threshold = 1e-12
 order = 19
 
+# Normalized diamagnetic frequencies. iota_e above is their ratio, so it is blind
+# to a common factor on both; these pin the drive itself.
+[quantities.slayer_Q_e]
+h5path = "Tearing/PerSurface/Q_e"
+type = "real_vector"
+extract = "all_real"
+label = "SLAYER Q_e"
+noise_threshold = 1e-12
+order = 20
+
+[quantities.slayer_Q_i]
+h5path = "Tearing/PerSurface/Q_i"
+type = "real_vector"
+extract = "all_real"
+label = "SLAYER Q_i"
+noise_threshold = 1e-12
+order = 21
+
 # 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/ω/γ

regression-harness/src/runner.jl

@@ -27,16 +27,18 @@ end
 Materialize the directory GPEC will actually run in.
 
 With no `overrides`, the example deck runs in place (returns it untouched). With overrides,
-the deck is copied to a throwaway temp dir and the named gpec.toml keys are patched there,
-so one shared example can serve several cases (e.g. a collisionless variant via
-`"KineticForces.nutype" => "zero"`). Override keys are dotted paths into the TOML; missing
-intermediate tables are created. Returns `(rundir, is_temp)`; the caller removes the temp
-tree when `is_temp`.
+the deck is copied to a uniquely named throwaway sibling of the example dir and the named gpec.toml keys are
+patched there, so one shared example can serve several cases (e.g. a collisionless variant via
+`"KineticForces.nutype" => "zero"`). The copy sits beside the original so relative file
+references in the deck (e.g. `SLAYER.profile_file = "../<other_example>/..."`) still resolve.
+Override keys are dotted paths into the TOML; missing intermediate tables are created.
+Returns `(rundir, is_temp)`; the caller removes `rundir` when `is_temp`.
 """
 function _materialize_rundir(example_path::String, overrides::Dict{String,Any})
     isempty(overrides) && return (example_path, false)
-    rundir = joinpath(mktempdir(), "deck")
-    cp(example_path, rundir)
+    # Unique per invocation, so concurrent runs from one checkout cannot share (and race on) a rundir.
+    rundir = mktempdir(dirname(example_path); prefix=".regress-override-" * basename(example_path) * "-", cleanup=false)
+    cp(example_path, rundir; force=true)
     rm(joinpath(rundir, "gpec.h5"); force=true)   # drop any stale output copied along
     cfg = TOML.parsefile(joinpath(rundir, "gpec.toml"))
     for (dotted, val) in overrides
@@ -523,7 +525,8 @@ function run_local(db::SQLite.DB, case_spec::CaseSpec, repo_root::String;
             rm(runinfo_file; force=true)
         end
         if rundir_is_temp && rundir !== nothing
-            rm(dirname(rundir); recursive=true, force=true)
+            # rundir sits beside the example; removing its parent would delete examples/ itself
+            rm(rundir; recursive=true, force=true)
         end
     end
 end
@@ -647,7 +650,8 @@ function run_at_commit(db::SQLite.DB, commit_hash::String, ref_name::String,
             rm(runinfo_file; force=true)
         end
         if rundir_is_temp && rundir !== nothing
-            rm(dirname(rundir); recursive=true, force=true)
+            # rundir sits beside the example; removing its parent would delete examples/ itself
+            rm(rundir; recursive=true, force=true)
         end
         if own_worktree && worktree_path !== nothing
             remove_worktree(worktree_path, repo_root)

src/InnerLayer/SLAYER/LayerInputs.jl

@@ -257,21 +257,20 @@ function build_slayer_inputs(equil, sings, profiles::KineticProfiles;
 
     _rs_at, _da_dpsi_at = radial_label(equil; rs_method=rs_method, theta=theta)
 
-    # Per-surface ω_*e, ω_*i (diamagnetic frequencies) from spline
-    # derivatives. When `compute_omega_star=true` we override any ω_*e/ω_*i
-    # carried in `profiles`. Main-ion density is
-    # taken equal to the electron density (quasi-neutrality, matching the
-    # staging step).
+    # Per-surface ω_*e, ω_*i (diamagnetic frequencies) from spline derivatives.
+    # In flux coordinates ω_* = n·(dp/dψ)/(e·n_e), so the surface's toroidal mode
+    # number enters here. Main-ion density is taken equal to the electron density
+    # (quasi-neutrality, matching the staging step).
     chi1 = 2π * equil.psio
-    _omega_star_at(ψ) = begin
+    _omega_star_at(ψ, n_tor) = begin
         n_e = Float64(profiles.n_e(ψ))
         dn_e = Float64(profiles.n_e(ψ; deriv=DerivOp(1)))
         T_e = Float64(profiles.T_e(ψ))
         dT_e = Float64(profiles.T_e(ψ; deriv=DerivOp(1)))
         T_i = Float64(profiles.T_i(ψ))
         dT_i = Float64(profiles.T_i(ψ; deriv=DerivOp(1)))
-        ω_star_e = (2π / chi1) * (T_e * dn_e / n_e + dT_e)
-        ω_star_i = -(2π / (Float64(z_i) * chi1)) * (T_i * dn_e / n_e + dT_i)
+        ω_star_e = n_tor * (2π / chi1) * (T_e * dn_e / n_e + dT_e)
+        ω_star_i = -n_tor * (2π / (Float64(z_i) * chi1)) * (T_i * dn_e / n_e + dT_i)
         return (ω_star_e, ω_star_i)
     end
 
@@ -285,20 +284,17 @@ function build_slayer_inputs(equil, sings, profiles::KineticProfiles;
         da_dpsi = _da_dpsi_at(psi)
         sval_r = r_based_shear(rs, q, q1, da_dpsi)
 
-        prof = profiles(psi)
-        # Override ω_*e, ω_*i with spline-derivative values when requested.
-        ω_e_use, ω_i_use = if compute_omega_star
-            _omega_star_at(psi)
-        else
-            (prof.omega_e, prof.omega_i)
-        end
-
         # Resonant (m, n): take the first element of the mode-number vectors.
         # Parallel-FM `sing.m`/`sing.n` hold exactly one entry each; ideal
         # DCON may hold multiple — we pick the first and document the choice.
         m_res = sing.m[1]
         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)
+
         # Pull geometric trapped-fraction inputs from ResistGeometry when
         # available (populated by ForceFreeStates.resist_eval_all!); else
         # fall back to nothing and let slayer_parameters compute them from

test/runtests_slayer_inputs.jl

@@ -210,4 +210,20 @@
         @test sl isa Vector{SLAYERParameters}
         @test isempty(sl)
     end
+
+    @testset "build_slayer_inputs: omega_star carries the toroidal mode number" begin
+        # In flux coordinates ω_* = n·(dp/dψ)/(e·n_e): resolving the same surface at n = 1
+        # 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)]
+        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)
+
+        @test sl2[1].Q_e / sl2[1].tauk ≈ 2 * (sl1[1].Q_e / sl1[1].tauk) rtol = 1e-12
+        @test sl2[1].Q_i / sl2[1].tauk ≈ 2 * (sl1[1].Q_i / sl1[1].tauk) rtol = 1e-12
+        @test sl2[1].iota_e ≈ sl1[1].iota_e rtol = 1e-12
+    end
+
 end