diff --git a/src/Tearing/Runner/Control.jl b/src/Tearing/Runner/Control.jl index 17128a3d6..1a12ff7e6 100644 --- a/src/Tearing/Runner/Control.jl +++ b/src/Tearing/Runner/Control.jl @@ -29,7 +29,9 @@ constructor. # Physics knobs - - `bt` -- toroidal field [T]. `nothing` → use `equil.config.b0exp` + - `bt` -- toroidal field `[T]`. `nothing` (default) resolves the physical + `B_T = F(ψ)/(2π·R₀)` per surface from the equilibrium's F-spline; a scalar or a + callable of `psi` overrides it - `mu_i` -- ion mass in proton-mass units (default 2.0 for D) - `zeff` -- effective charge - `chi_perp`, `chi_tor` -- fallback perpendicular / toroidal heat diff --git a/src/Tearing/Runner/run_slayer.jl b/src/Tearing/Runner/run_slayer.jl index b42980465..2a29c63a8 100644 --- a/src/Tearing/Runner/run_slayer.jl +++ b/src/Tearing/Runner/run_slayer.jl @@ -360,6 +360,9 @@ read from `control.profile_file` (relative to `dir_path`) through the shared profiles they set χ⊥(ψ)/χ_φ(ψ), otherwise the scalar `control.chi_perp`/ `chi_tor` fallbacks are used. +The toroidal field comes from `control.bt`; leaving it unset (the default) makes +`build_slayer_inputs` evaluate the physical `B_T = F(ψ)/(2π·R₀)` per surface. + Returns an `enabled=false` `SLAYERResult` when `control.enabled` is false. """ @@ -394,7 +397,12 @@ function run_slayer(equil, surfaces::AbstractVector, delta_prime_matrix::Abstrac resistivity_model=_build_resistivity_model(control.resistivity_model), lnLambda_form=control.lnLambda_form) else - bt = control.bt === nothing ? equil.config.b0exp : control.bt + # `equil.config.b0exp` is a NORMALIZATION (commonly exactly 1.0), not the toroidal + # field, so substituting it here silently ran the layer physics at B_T = 1 T. Pass the + # control value through instead: `nothing` makes build_slayer_inputs compute the + # physical B_T = F(psi)/(2*pi*R_0) per surface from the equilibrium's F-spline, which is + # what its docstring already prescribes. + bt = control.bt # χ⊥/χ_φ from the kinetic file when present, else the scalar fallbacks. chi_perp = loaded.chi_perp === nothing ? control.chi_perp : loaded.chi_perp chi_tor = loaded.chi_tor === nothing ? control.chi_tor : loaded.chi_tor diff --git a/test/runtests_slayer_inputs.jl b/test/runtests_slayer_inputs.jl index c0b7d2d7d..55d6124e9 100644 --- a/test/runtests_slayer_inputs.jl +++ b/test/runtests_slayer_inputs.jl @@ -106,6 +106,34 @@ @test sl[1].Q_i == -sl[1].tauk * profiles.omega_i(0.3) end + @testset "build_slayer_inputs: bt defaults to the physical B_T, not a normalization" begin + # `b0exp` is a normalization, not a field; passing it as `bt` ran the layer physics at the + # wrong toroidal field. Asserted behaviourally rather than by grepping the source, so a + # refactor cannot silently break the check and a reintroduction anywhere in the chain still + # fails: tau_h = R0*sqrt(mu0*rho)/(n*sval_r*bt), so lu is linear in the field actually used + # and the resolved bt is recoverable from the returned parameters. + # + # Note this fixture cannot exhibit the original bug: Solovev is normalized so that + # b0exp == 1.0 and F(psi)/(2*pi*R0) == 1.0 to roundoff, which is exactly why the defect + # 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)] + 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) + + # 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 + + # An explicit bt overrides, and lu tracks it linearly -- so the field that was actually + # used is what the parameters encode, not merely some field. + @test sl_double[1].lu / sl_default[1].lu ≈ 2.0 rtol = 1e-10 + @test sl_double[1].lu != sl_default[1].lu + 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)] diff --git a/test/runtests_slayer_runner.jl b/test/runtests_slayer_runner.jl index fb4b923dd..0eacc18ed 100644 --- a/test/runtests_slayer_runner.jl +++ b/test/runtests_slayer_runner.jl @@ -306,4 +306,5 @@ end end end + end