diff --git a/.gitignore b/.gitignore index fcb3cc96c..1d1bf86a2 100644 --- a/.gitignore +++ b/.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-* diff --git a/regression-harness/cases/diiid_slayer_n1.toml b/regression-harness/cases/diiid_slayer_n1.toml index 540ca7e50..90edf2e09 100644 --- a/regression-harness/cases/diiid_slayer_n1.toml +++ b/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/ω/γ diff --git a/regression-harness/src/runner.jl b/regression-harness/src/runner.jl index ed7e5f4d5..1a192a140 100644 --- a/regression-harness/src/runner.jl +++ b/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 = "..//..."`) 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) diff --git a/src/InnerLayer/SLAYER/LayerInputs.jl b/src/InnerLayer/SLAYER/LayerInputs.jl index cb1d14b85..252e81ae9 100644 --- a/src/InnerLayer/SLAYER/LayerInputs.jl +++ b/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 diff --git a/test/runtests_slayer_inputs.jl b/test/runtests_slayer_inputs.jl index a5a6da414..fce001036 100644 --- a/test/runtests_slayer_inputs.jl +++ b/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