Skip to content
Open
Show file tree
Hide file tree
Changes from all commits
Commits
File filter

Filter by extension

Filter by extension

Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
3 changes: 3 additions & 0 deletions .gitignore
Original file line number Diff line number Diff line change
Expand Up @@ -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-*
173 changes: 173 additions & 0 deletions regression-harness/cases/diiid_slayer_n2.toml
Original file line number Diff line number Diff line change
@@ -0,0 +1,173 @@
# Regression case: DIII-D-like H-mode n=2 SLAYER tearing-mode growth rates.
# Reuses the n=1 SLAYER deck via [overrides] at toroidal mode number 2, where the deck's
# qmin = 1.204 puts the 3/2 surface in the plasma. Exists to guard the n-dependence of the
# layer build — in particular the diamagnetic frequencies ω_* ∝ n feeding Q_e/Q_i, which
# no n=1 case can constrain — with the inner rationals 3/2, 4/2, 5/2 as the tracked roots.
# Each [quantities.*] block names an HDF5 path in the run output, how to extract it,
# and the noise floor below which a difference is treated as zero.
[case]
name = "diiid_slayer_n2"

Copy link
Copy Markdown
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

So the entire purpose of this regression case is to guard the n-dependence of the layer module... that seems like overkill for a regression test that takes 5 minutes to run. I speak selfishly because my (and I'm assuming other's) regress running is automatically chosen by Claude. So if I change anything related to this module's code, Claude will queue up both of these cases to test something relatively small.

So personally, I am against regression case bloat and don't want them to take the place of unit tests. So my question here for your consideration - if the n1 and n2 case here actually needed and worth the additional regression harness runtime, which can already get quite long? Or is this something where we could keep a single regression case (perhaps n2 to try to catch n-dependent things like this in the future) and supplement with unit tests that check for specifics to slayer more cheaply?

description = "DIII-D-like H-mode equilibrium, n=2, SLAYER tearing-mode analysis guarding the 3/2 surface and the n-scaling of Q_e/Q_i (reuses the n=1 deck via overrides)"
example_dir = "examples/DIIID-like_SLAYER_example"

# Run the shared SLAYER deck at n = 2 (no separate n=2 example directory).
[overrides]
"ForceFreeStates.nn_low" = 2
"ForceFreeStates.nn_high" = 2

# Per-surface SLAYER layer parameters (geometry + dimensionless)
[quantities.slayer_ising]
h5path = "Tearing/PerSurface/rational_index"
type = "real_vector"
extract = "all_real"
label = "SLAYER surface indices"
noise_threshold = 0
order = 10

[quantities.slayer_m]
h5path = "Tearing/PerSurface/m"
type = "real_vector"
extract = "all_real"
label = "SLAYER poloidal m"
noise_threshold = 0
order = 11

[quantities.slayer_n]
h5path = "Tearing/PerSurface/n"
type = "real_vector"
extract = "all_real"
label = "SLAYER toroidal n"
noise_threshold = 0
order = 12

[quantities.slayer_rs]
h5path = "Tearing/PerSurface/rs"
type = "real_vector"
extract = "all_real"
label = "SLAYER minor radius rs"
noise_threshold = 1e-10
order = 13

[quantities.slayer_sval_r]
h5path = "Tearing/PerSurface/sval_r"
type = "real_vector"
extract = "all_real"
label = "SLAYER r-based shear"
noise_threshold = 1e-10
order = 14

[quantities.slayer_lu]
h5path = "Tearing/PerSurface/lu"
type = "real_vector"
extract = "all_real"
label = "SLAYER Lundquist S"
noise_threshold = 1e-8
order = 15

[quantities.slayer_D_norm]
h5path = "Tearing/PerSurface/D_norm"
type = "real_vector"
extract = "all_real"
label = "SLAYER D_norm"
noise_threshold = 1e-10
order = 16

[quantities.slayer_P_perp]
h5path = "Tearing/PerSurface/P_perp"
type = "real_vector"
extract = "all_real"
label = "SLAYER P_perp"
noise_threshold = 1e-8
order = 17

[quantities.slayer_tauk]
h5path = "Tearing/PerSurface/tau_k"
type = "real_vector"
extract = "all_real"
label = "SLAYER tauk"
noise_threshold = 1e-12
order = 18

[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

[quantities.slayer_iota_e]
h5path = "Tearing/PerSurface/iota_e"
type = "real_vector"
extract = "all_real"
label = "SLAYER iota_e"
noise_threshold = 1e-12
order = 19

# 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/ω/γ
# thresholds are absolute and modestly loose; re-pin intentionally if a
# solver/AMR change shifts the root inventory.
# Growth-rate / frequency / root outputs. Pinned for the inner three rational
# surfaces only (3/2, 4/2, 5/2) via `first_3`: the Δ'/γ contour search is
# numerically unreliable on the outermost surfaces near the edge, so those are
# deliberately not golden-tracked.
[quantities.slayer_Q]
h5path = "Tearing/Roots/Q_root"
type = "complex_vector"
extract = "first_3_complex"
label = "SLAYER Q_root [3/2,4/2,5/2]"
noise_threshold = 1e-4
order = 30

[quantities.slayer_omega_Hz]
h5path = "Tearing/Roots/omega"
type = "real_vector"
extract = "first_3"
label = "SLAYER ω_Hz [3/2,4/2,5/2]"
noise_threshold = 1.0
order = 32

[quantities.slayer_gamma_Hz]
h5path = "Tearing/Roots/gamma"
type = "real_vector"
extract = "first_3"
label = "SLAYER γ_Hz [3/2,4/2,5/2]"
noise_threshold = 1e-1
order = 33

# no_root flag (1 = extraction failed). Pinned for the inner three surfaces
# so the 3/2, 4/2, 5/2 root inventory cannot drift silently.
[quantities.slayer_no_root]
h5path = "Tearing/Roots/no_root"
type = "real_vector"
extract = "first_3"
label = "SLAYER no_root flags [3/2,4/2,5/2]"
noise_threshold = 0
order = 34

# 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

[quantities.runtime]
h5path = ""
type = "runtime"
extract = "value"
label = "Runtime (s)"
noise_threshold = 0.0
order = 999
20 changes: 12 additions & 8 deletions regression-harness/src/runner.jl
Original file line number Diff line number Diff line change
Expand Up @@ -27,15 +27,17 @@ 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 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")
rundir = joinpath(dirname(example_path), ".regress-override-" * basename(example_path))
rm(rundir; recursive=true, force=true) # stale leftover from a crashed run
cp(example_path, rundir)
rm(joinpath(rundir, "gpec.h5"); force=true) # drop any stale output copied along
cfg = TOML.parsefile(joinpath(rundir, "gpec.toml"))
Expand Down Expand Up @@ -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
Expand Down Expand Up @@ -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)
Expand Down
19 changes: 11 additions & 8 deletions src/InnerLayer/SLAYER/LayerInputs.jl
Original file line number Diff line number Diff line change
Expand Up @@ -217,20 +217,23 @@ 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)
# Override ω_*e, ω_*i with spline-derivative values when requested. In flux coordinates

Copy link
Copy Markdown
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

When I see override, I think of something that is already set; in this case its just another branch of how to set this value.

I think more true to form would be

ω_e_use, ω_i_use = (prof.omega_e, prof.omega_i)
# Override ω_*e, ω_*i with spline-derivative values when requested
compute_omega_star && (ω_e_use, ω_i_use = _omega_star_at(psi, n))

or something like that. The PR comment does suggest that prof.omega_e/i are never actually used, but I didn't factor that in here since I am unfamiliar with the rest of the relevant code

# ω_* = n·(dp/dψ)/(e·n_e); `_omega_star_at` returns the n = 1 value, so restore the
# factor n here. Values supplied through `profiles` are taken to be physical already.
ω_e_use, ω_i_use = if compute_omega_star

Copy link
Copy Markdown
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Personally, I think the something = an entire if-else block is kinda confusing code, but it seems like that's in a lot of this file so its probably ok.

But if I am understanding this PR correctly - the actual omega* has an included n, so its a little bit deceptive to have _omega_star_at return something that isn't actually omega* and needs to be multiplied by n. I think its more clear to build it into the helper even though Claude argues against it _omega_star_at(psi, n). Then you can remove this 3 line comment up here and just leave it as is

ωe1, ωi1 = _omega_star_at(psi)
(n_res * ωe1, n_res * ωi1)
else
(prof.omega_e, prof.omega_i)
end

# 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
Expand Down
16 changes: 16 additions & 0 deletions test/runtests_slayer_inputs.jl
Original file line number Diff line number Diff line change
Expand Up @@ -183,4 +183,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
Loading