From 73e85a93e9617e6d302cf7c5ef496122024baeb7 Mon Sep 17 00:00:00 2001 From: Jake Halpern Date: Thu, 13 Aug 2026 09:07:11 -0400 Subject: [PATCH 1/9] GPEC - REFACTOR - Freeze the PerturbedEquilibrium, KineticForces and ForcingTerms control structs Continues the issue #139 struct cleanup that already froze ForceFreeStatesControl: a control struct built from the input TOML should be built once and then left alone, so a reader at a call site can tell nothing behind it changes. - PerturbedEquilibriumControl and KineticForcesControl had no writers in src/ at all. - ForcingTermsControl was filled in two steps because [[ForcingTerms.coil_set]] parses to a Vector{Dict} that cannot ride along in the scalar kwarg splat; coil_sets_raw is now passed as an explicit keyword after the splat instead of assigned afterwards. - SLAYERControl was already immutable; no change needed. - The one test that mutated a control (check_psi_quadrature_convergence) builds a second control with the nonzero atol_psi rather than rewriting the first. No intended numerical change. Co-Authored-By: Claude Opus 5 --- src/ForcingTerms/ForcingTerms.jl | 2 +- src/GeneralizedPerturbedEquilibrium.jl | 5 ++--- src/KineticForces/KineticForcesStructs.jl | 2 +- src/PerturbedEquilibrium/PerturbedEquilibriumStructs.jl | 2 +- test/runtests_kinetic.jl | 4 ++-- 5 files changed, 7 insertions(+), 8 deletions(-) diff --git a/src/ForcingTerms/ForcingTerms.jl b/src/ForcingTerms/ForcingTerms.jl index 1869053f9..88f87fdb5 100644 --- a/src/ForcingTerms/ForcingTerms.jl +++ b/src/ForcingTerms/ForcingTerms.jl @@ -28,7 +28,7 @@ Coil settings (used when `forcing_data_format = "coil"`): - `nzeta_coil::Int` - Toroidal grid resolution; 0 = auto (32 × n) - `coil_sets_raw::Vector{Dict{String,Any}}` - Parsed `[[ForcingTerms.coil_set]]` TOML blocks """ -Base.@kwdef mutable struct ForcingTermsControl +@kwdef struct ForcingTermsControl # Forcing data file settings (ascii/hdf5 formats) forcing_data_file::String = "forcing.dat" forcing_data_format::String = "ascii" diff --git a/src/GeneralizedPerturbedEquilibrium.jl b/src/GeneralizedPerturbedEquilibrium.jl index 94a7e0e04..6d7266f3b 100755 --- a/src/GeneralizedPerturbedEquilibrium.jl +++ b/src/GeneralizedPerturbedEquilibrium.jl @@ -556,13 +556,12 @@ function main_from_inputs( if "ForcingTerms" in keys(inputs) forcing_raw = inputs["ForcingTerms"] # [[ForcingTerms.coil_set]] becomes a Vector{Dict} — must be excluded from - # kwarg splatting and handled separately via coil_sets_raw field + # kwarg splatting and passed as the explicit coil_sets_raw keyword coil_sets_raw = Vector{Dict{String,Any}}(get(forcing_raw, "coil_set", Dict{String,Any}[])) scalar_forcing = filter(p -> p.first != "coil_set", forcing_raw) ft_ctrl = ForcingTerms.ForcingTermsControl(; - (Symbol(k) => v for (k, v) in scalar_forcing)... + (Symbol(k) => v for (k, v) in scalar_forcing)..., coil_sets_raw=coil_sets_raw ) - ft_ctrl.coil_sets_raw = coil_sets_raw else ft_ctrl = ForcingTerms.ForcingTermsControl() # Use defaults end diff --git a/src/KineticForces/KineticForcesStructs.jl b/src/KineticForces/KineticForcesStructs.jl index 7e533d8e9..83519909b 100644 --- a/src/KineticForces/KineticForcesStructs.jl +++ b/src/KineticForces/KineticForcesStructs.jl @@ -60,7 +60,7 @@ Constructed via keyword arguments or from a TOML dict: ctrl = KineticForcesControl(; (Symbol(k) => v for (k, v) in inputs["KineticForces"])...) ``` """ -@kwdef mutable struct KineticForcesControl +@kwdef struct KineticForcesControl # Moment type moment::String = "pressure" # "heat" or "pressure" diff --git a/src/PerturbedEquilibrium/PerturbedEquilibriumStructs.jl b/src/PerturbedEquilibrium/PerturbedEquilibriumStructs.jl index 96f982386..1ce0895d6 100644 --- a/src/PerturbedEquilibrium/PerturbedEquilibriumStructs.jl +++ b/src/PerturbedEquilibrium/PerturbedEquilibriumStructs.jl @@ -29,7 +29,7 @@ Regularization: # High Priority (MWE) - `reg_spot::Float64` - Regularization width for singular surface smoothing (default: 0.05). Set to 0 to disable. Must be ≥ 0. """ -@kwdef mutable struct PerturbedEquilibriumControl +@kwdef struct PerturbedEquilibriumControl # High Priority (MWE) fixed_boundary::Bool = false output_eigenmodes::Bool = true diff --git a/test/runtests_kinetic.jl b/test/runtests_kinetic.jl index f23d10315..29955c2c7 100644 --- a/test/runtests_kinetic.jl +++ b/test/runtests_kinetic.jl @@ -411,8 +411,8 @@ # Hit maxevals: error above tolerance @test_logs (:warn, r"maxevals_psi") KF.check_psi_quadrature_convergence(total, 0.5, ctrl, "fgar") # Nonzero atol_psi dominating a small torque: the silent-garbage scenario - ctrl.atol_psi = 1e-2 - @test_logs (:warn, r"atol_psi") KF.check_psi_quadrature_convergence(1e-3 + 0.0im, 1e-3, ctrl, "fgar") + ctrl_atol = KF.KineticForcesControl(; atol_psi=1e-2) + @test_logs (:warn, r"atol_psi") KF.check_psi_quadrature_convergence(1e-3 + 0.0im, 1e-3, ctrl_atol, "fgar") end @testset "METHOD_REGISTRY" begin From a0520f47fdefc91feeb0bb0cf6173cdf34e5af91 Mon Sep 17 00:00:00 2001 From: Jake Halpern Date: Fri, 14 Aug 2026 14:28:10 -0400 Subject: [PATCH 2/9] EQUIL - REFACTOR - Freeze EquilibriumConfig by resolving psihigh onto the run inputs EquilibriumConfig is now an immutable struct: it holds the input request and is never written to after construction. The separatrix clamp previously wrote its result back into config.psihigh, which was load-bearing through aliasing (PlasmaEquilibrium.config is raw_profile.config). The resolved value now lives on the run-input structs as psihigh_resolved, defaulted to config.psihigh by an outer constructor so existing readers construct unchanged. resolve_psihigh! replaces the two duplicated clamp blocks in setup_equilibrium, the three solvers form their psi grid from psihigh_resolved, and the efit_by_inversion mid-solve InverseRunInput rebuild forwards it explicitly instead of relying on the config alias. EquilibriumParameters.psihigh_resolved carries it onto the finished equilibrium for sing_lim! and the driver. Two other config mutations had to go for the struct to freeze: the eq_filename path resolution now folds into the dict before construction, and the rerun path clears eq_filename on a dict copy so the re-serialized inputs keep the user's original. EquilibriumConfig(path::String) delegates to the dict constructor rather than duplicating it. setup_equilibrium's signature is unchanged, so no caller, harness or benchmark entry point moves. No intended numerical change. Co-Authored-By: Claude Opus 5 --- benchmarks/benchmark_q_vs_iota_edge.jl | 6 +- src/Equilibrium/DirectEquilibrium.jl | 3 +- .../DirectEquilibriumByInversion.jl | 24 ++++++-- src/Equilibrium/Equilibrium.jl | 23 ++------ src/Equilibrium/EquilibriumTypes.jl | 57 +++++++++++++------ src/Equilibrium/InverseEquilibrium.jl | 4 +- src/ForceFreeStates/Sing.jl | 2 +- src/GeneralizedPerturbedEquilibrium.jl | 7 +-- src/Rerun.jl | 7 ++- test/runtests_equil.jl | 25 ++++++++ test/runtests_sing.jl | 4 +- 11 files changed, 108 insertions(+), 54 deletions(-) diff --git a/benchmarks/benchmark_q_vs_iota_edge.jl b/benchmarks/benchmark_q_vs_iota_edge.jl index e258fb457..2dfb89eb6 100644 --- a/benchmarks/benchmark_q_vs_iota_edge.jl +++ b/benchmarks/benchmark_q_vs_iota_edge.jl @@ -20,9 +20,9 @@ const EXAMPLE_DIR = joinpath(@__DIR__, "..", "examples", "DIIID-like_ideal_examp # Dense ldp reference equilibrium: treat its q(ψ) as ground truth function reference_q() - _, eq_config, additional_input = GPE.build_inputs_from_toml(EXAMPLE_DIR) - eq_config.grid_type = "ldp" - eq_config.mpsi = 1024 + inputs, _, additional_input = GPE.build_inputs_from_toml(EXAMPLE_DIR) + equil_dict = merge(inputs["Equilibrium"], Dict{String,Any}("grid_type" => "ldp", "mpsi" => 1024)) + eq_config = GPE.Equilibrium.EquilibriumConfig(equil_dict, EXAMPLE_DIR) equil = GPE.Equilibrium.setup_equilibrium(eq_config, additional_input) return equil, eq_config end diff --git a/src/Equilibrium/DirectEquilibrium.jl b/src/Equilibrium/DirectEquilibrium.jl index b98511ab9..adf298651 100644 --- a/src/Equilibrium/DirectEquilibrium.jl +++ b/src/Equilibrium/DirectEquilibrium.jl @@ -474,7 +474,7 @@ robustness. psio = raw_profile.psio mtheta = equil_params.mtheta psilow = equil_params.psilow - psihigh = equil_params.psihigh + psihigh = raw_profile.psihigh_resolved # Locate the magnetic axis and separatrix for the field-line integrations ro, zo, _, rs2 = direct_position!(raw_profile) @@ -647,6 +647,7 @@ robustness. params = EquilibriumParameters() params.bt_sign = raw_profile.bt_sign + params.psihigh_resolved = psihigh return PlasmaEquilibrium(raw_profile.config, params, profiles, geometry, rzphi_xs, rzphi_ys, diff --git a/src/Equilibrium/DirectEquilibriumByInversion.jl b/src/Equilibrium/DirectEquilibriumByInversion.jl index 252ae515b..136d8e6ac 100644 --- a/src/Equilibrium/DirectEquilibriumByInversion.jl +++ b/src/Equilibrium/DirectEquilibriumByInversion.jl @@ -73,12 +73,12 @@ end """ clamp_psihigh_to_separatrix(raw_profile) -> (clamped_psihigh, was_adjusted) -Binary-searches for the highest psihigh ≤ raw_profile.config.psihigh at which the +Binary-searches for the highest psihigh ≤ `raw_profile.psihigh_resolved` at which the ψ level set is still a closed curve in the EFIT grid. Returns the safe value and a Bool indicating whether any clamping occurred. """ function clamp_psihigh_to_separatrix(raw_profile::DirectRunInput) - psihigh = raw_profile.config.psihigh + psihigh = raw_profile.psihigh_resolved ψ_coarse = raw_profile.psi_in.nodal_derivs.partials[1, :, :] has_closed_contour(ψ_high) = any( @@ -99,6 +99,22 @@ function clamp_psihigh_to_separatrix(raw_profile::DirectRunInput) return (lo, true) end +""" + resolve_psihigh!(raw_profile) -> raw_profile + +Clamp `raw_profile.psihigh_resolved` to the outermost closed flux surface, warning if the +requested `psihigh` had none. Mutates and returns the input; `config.psihigh` keeps the request. +""" +function resolve_psihigh!(raw_profile::DirectRunInput) + psihigh_safe, adjusted = clamp_psihigh_to_separatrix(raw_profile) + if adjusted + @warn "psihigh=$(raw_profile.psihigh_resolved) has no closed flux surface in EFIT grid; " * + "clamped to $(round(psihigh_safe; sigdigits=7))" + raw_profile.psihigh_resolved = psihigh_safe + end + return raw_profile +end + """ resample_contour_to_theta_grid!(R_out, Z_out, curve, ro, zo, theta_grid) @@ -408,7 +424,7 @@ function equilibrium_solver_by_inversion( psio = raw_profile.psio mtheta = equil_params.mtheta psilow = equil_params.psilow - psihigh = equil_params.psihigh + psihigh = raw_profile.psihigh_resolved # Locate the magnetic axis and separatrix for the contour tracing ro, zo, _, rs2 = direct_position!(raw_profile) @@ -672,7 +688,7 @@ function equilibrium_solver_by_inversion( # Intermediate inverse input; the captured DirectIngest rides on the original eq_input, # which setup_equilibrium forwards onto the equilibrium, so this one carries ingest=nothing. inv_input = InverseRunInput(raw_profile.config, raw_profile.sq_in, - rz_in_xs, rz_in_ys, rz_in_R, rz_in_Z, ro, zo, psio, nothing) + rz_in_xs, rz_in_ys, rz_in_R, rz_in_Z, ro, zo, psio, nothing, raw_profile.psihigh_resolved) pe = equilibrium_solver(inv_input; override_psi_nodes) diff --git a/src/Equilibrium/Equilibrium.jl b/src/Equilibrium/Equilibrium.jl index 03935640e..2fa68d053 100644 --- a/src/Equilibrium/Equilibrium.jl +++ b/src/Equilibrium/Equilibrium.jl @@ -92,27 +92,16 @@ function setup_equilibrium(eq_config::EquilibriumConfig, additional_input=nothin if additional_input isa DirectRunInput eq_input = additional_input eq_input.config = eq_config - # Re-run the separatrix clamp for efit-family replays so an overridden - # psihigh from the rerun TOML is re-validated against the closed flux region. - if eq_type in EFIT_KINDS - psihigh_safe, adjusted = clamp_psihigh_to_separatrix(eq_input) - if adjusted - @warn "psihigh=$(eq_input.config.psihigh) has no closed flux surface in EFIT grid; " * - "clamped to $(round(psihigh_safe; sigdigits=7))" - eq_input.config.psihigh = psihigh_safe - end - end + # Reset before re-resolving: a pass-1 clamped value must not shadow a psihigh + # overridden in the rerun TOML. + eq_input.psihigh_resolved = eq_config.psihigh + eq_type in EFIT_KINDS && resolve_psihigh!(eq_input) elseif additional_input isa InverseRunInput eq_input = additional_input eq_input.config = eq_config + eq_input.psihigh_resolved = eq_config.psihigh elseif eq_type in EFIT_KINDS - eq_input = read_efit(eq_config) - psihigh_safe, adjusted = clamp_psihigh_to_separatrix(eq_input) - if adjusted - @warn "psihigh=$(eq_input.config.psihigh) has no closed flux surface in EFIT grid; " * - "clamped to $(round(psihigh_safe; sigdigits=7))" - eq_input.config.psihigh = psihigh_safe - end + eq_input = resolve_psihigh!(read_efit(eq_config)) elseif eq_type in ["chease2", "chease_ascii"] eq_input = read_chease_ascii(eq_config) elseif eq_type in ["chease", "chease_binary"] diff --git a/src/Equilibrium/EquilibriumTypes.jl b/src/Equilibrium/EquilibriumTypes.jl index 4bb29f8f2..5ec2a82ce 100644 --- a/src/Equilibrium/EquilibriumTypes.jl +++ b/src/Equilibrium/EquilibriumTypes.jl @@ -5,8 +5,8 @@ end """ EquilibriumConfig(...) -A mutable struct containing configuration parameters for equilibrium reconstruction. -Bundles all necessary settings originally specified in the equil fortran namelists. +An immutable struct containing configuration parameters for equilibrium reconstruction +specified in the input. ## Fields @@ -27,7 +27,11 @@ Bundles all necessary settings originally specified in the equil fortran namelis refinement when mpsi=0, three-region log layout when mpsi>0; "ldp", "pow1", "uniform"; "log_asymptotic" is a legacy alias for "auto") - `psilow::Float64` - Lower limit of normalized flux coordinate - - `psihigh::Float64` - Upper limit of normalized flux coordinate + - `psihigh::Float64` - Requested upper limit of normalized flux coordinate. For efit-family + equilibria this is the user's request, which may lie outside the closed-flux region; the + value the equilibrium is actually formed on is `DirectRunInput.psihigh_resolved` (carried + onto the equilibrium as `EquilibriumParameters.psihigh_resolved`). Read that, not this, + for anything downstream of `setup_equilibrium`. - `mpsi::Int` - Number of radial grid intervals; 0 with grid_type="auto" selects the two-pass auto grid: the main driver forms a coarse pass-1 equilibrium, measures its curvature, pins knots on rational surfaces, and re-forms on the refined grid. Standalone `setup_equilibrium` @@ -40,7 +44,7 @@ Bundles all necessary settings originally specified in the equil fortran namelis - `force_termination::Bool` - Terminate after equilibrium setup (skip stability calculations) - `use_galgrid::Bool` - Use the same grid as galerkin method """ -@kwdef mutable struct EquilibriumConfig +@kwdef struct EquilibriumConfig eq_type::String = "efit" eq_filename::String = "mypath" r0exp::Float64 = 1.0 @@ -176,15 +180,14 @@ function EquilibriumConfig(equil_dict::Dict{String,Any}, base_path::String="./") end end - # Construct validated struct - config = EquilibriumConfig(; symbolize_keys(config_data)...) - # Only resolve `eq_filename` against `base_path` if the user actually - # supplied one (otherwise leave the kwdef sentinel for the embedded path). - if haskey(config_data, "eq_filename") && !isabspath(config.eq_filename) - config.eq_filename = normpath(joinpath(base_path, config.eq_filename)) + # Only resolve `eq_filename` against `base_path` if the user actually supplied one + # (otherwise leave the kwdef sentinel for the embedded path). The empty string is the + # rerun path's "no input file" marker and must stay empty, not become `base_path`. + if haskey(config_data, "eq_filename") && !isempty(config_data["eq_filename"]) && !isabspath(config_data["eq_filename"]) + config_data["eq_filename"] = normpath(joinpath(base_path, config_data["eq_filename"])) end - return config + return EquilibriumConfig(; symbolize_keys(config_data)...) end """ @@ -208,12 +211,7 @@ function EquilibriumConfig(path::String) end # Construct validated struct - config = EquilibriumConfig(; symbolize_keys(config_data)...) - if !isabspath(config.eq_filename) - config.eq_filename = normpath(joinpath(dirname(path), config.eq_filename)) - end - - return config + return EquilibriumConfig(Dict{String,Any}(config_data), dirname(path)) end """ @@ -444,6 +442,10 @@ raw equilibrium data and preparing the initial splines. - `bt_sign::Int` — Sign of the toroidal field (+1 or -1); read from fpol sign in EFIT g-files - `ingest::EquilibriumIngest` — captured raw arrays for the `gpec.h5` rerun snapshot (a [`DirectIngest`](@ref) for file-based reads, or `nothing` for analytic equilibria) + - `psihigh_resolved::Float64` — outer flux limit the equilibrium is formed on: `config.psihigh` + clamped to the outermost closed flux surface by [`resolve_psihigh!`](@ref). Defaults to + `config.psihigh` and only differs for efit-family equilibria whose requested limit falls + outside the closed-flux region. The solvers build their ψ grid from this field. """ mutable struct DirectRunInput{S<:FastInterpolations.CubicSeriesInterpolant,I2D<:FastInterpolations.CubicInterpolantND} config::EquilibriumConfig @@ -458,8 +460,16 @@ mutable struct DirectRunInput{S<:FastInterpolations.CubicSeriesInterpolant,I2D<: psio::Float64 # The total flux difference |ψ_axis - ψ_boundary| [Weber / radian]. bt_sign::Int # Sign of the toroidal field: +1 or -1 (from fpol sign in g-file) ingest::EquilibriumIngest + psihigh_resolved::Float64 end +# Readers construct without a resolved psihigh; it starts at the request and `resolve_psihigh!` +# clamps it for efit-family equilibria. +DirectRunInput(config::EquilibriumConfig, sq_in, psi_in, psi_in_xs, psi_in_ys, + rmin, rmax, zmin, zmax, psio, bt_sign, ingest) = + DirectRunInput(config, sq_in, psi_in, psi_in_xs, psi_in_ys, + rmin, rmax, zmin, zmax, psio, bt_sign, ingest, config.psihigh) + """ InverseRunInput(...) @@ -478,6 +488,9 @@ A container struct for inputs to the `inverse_run` function. - `psio::Float64` - Total flux difference |ψ_axis - ψ_boundary| [Wb/rad] - `ingest::EquilibriumIngest` - captured raw arrays for the `gpec.h5` rerun snapshot (an [`InverseIngest`](@ref) for file-based reads, or `nothing` for analytic equilibria) + - `psihigh_resolved::Float64` - outer flux limit the equilibrium is formed on; see + [`DirectRunInput`](@ref). Equals `config.psihigh` for every inverse reader (CHEASE, + analytic); `efit_by_inversion` forwards the clamped value from its `DirectRunInput`. """ mutable struct InverseRunInput{S<:FastInterpolations.CubicSeriesInterpolant,I2D<:FastInterpolations.CubicInterpolantND} config::EquilibriumConfig @@ -490,8 +503,14 @@ mutable struct InverseRunInput{S<:FastInterpolations.CubicSeriesInterpolant,I2D< zo::Float64 # Z axis location psio::Float64 # Total flux difference |psi_axis - psi_boundary| ingest::EquilibriumIngest + psihigh_resolved::Float64 end +InverseRunInput(config::EquilibriumConfig, sq_in, rz_in_xs, rz_in_ys, rz_in_R, rz_in_Z, + ro, zo, psio, ingest) = + InverseRunInput(config, sq_in, rz_in_xs, rz_in_ys, rz_in_R, rz_in_Z, + ro, zo, psio, ingest, config.psihigh) + """ EquilibriumParameters @@ -502,6 +521,9 @@ A mutable struct containing computed equilibrium parameters and diagnostic flags - `ro::Union{Nothing,Float64}` - R-coordinate of the magnetic axis [m] - `zo::Union{Nothing,Float64}` - Z-coordinate of the magnetic axis [m] - `psio::Union{Nothing,Float64}` - Total flux difference |ψ_axis - ψ_boundary| [Wb/rad] + - `psihigh_resolved::Union{Nothing,Float64}` - Outer flux limit the equilibrium was formed on, + equal to the outermost ψ node. This is `config.psihigh` clamped to the outermost closed flux + surface; downstream code wanting the plasma edge must read this, not `config.psihigh`. - `rsep::Union{Nothing,Vector{Float64}}` - R-coordinates of the plasma boundary [m] - `zsep::Union{Nothing,Vector{Float64}}` - Z-coordinates of the plasma boundary [m] - `rext::Union{Nothing,Vector{Float64}}` - R-coordinates of the plasma edge [m] @@ -554,6 +576,7 @@ A mutable struct containing computed equilibrium parameters and diagnostic flags ro::Union{Nothing,Float64} = nothing # R-coordinate of the magnetic axis [m] zo::Union{Nothing,Float64} = nothing # Z-coordinate of the magnetic axis [m] psio::Union{Nothing,Float64} = nothing # Total flux difference |ψ_axis - ψ_boundary| [Wb/rad] + psihigh_resolved::Union{Nothing,Float64} = nothing # Outer flux limit actually formed on (clamped config.psihigh) rsep::Union{Nothing,Vector{Float64}} = nothing # R-coordinates of the plasma boundary [m] zsep::Union{Nothing,Vector{Float64}} = nothing # Z-coordinates of the plasma boundary [m] rext::Union{Nothing,Vector{Float64}} = nothing # R-coordinates of the plasma edge [m] diff --git a/src/Equilibrium/InverseEquilibrium.jl b/src/Equilibrium/InverseEquilibrium.jl index 01256e8ce..ff78c9a3d 100644 --- a/src/Equilibrium/InverseEquilibrium.jl +++ b/src/Equilibrium/InverseEquilibrium.jl @@ -59,7 +59,7 @@ function equilibrium_solver(input::InverseRunInput; override_psi_nodes::Union{No mpsi = config.mpsi mtheta = config.mtheta psilow = config.psilow - psihigh = config.psihigh + psihigh = input.psihigh_resolved newq0 = config.newq0 # c----------------------------------------------------------------------- @@ -386,7 +386,7 @@ function equilibrium_solver(input::InverseRunInput; override_psi_nodes::Union{No return PlasmaEquilibrium( input.config, - EquilibriumParameters(), + EquilibriumParameters(; psihigh_resolved=psihigh), profiles, geometry, rzphi_xs, rzphi_ys, diff --git a/src/ForceFreeStates/Sing.jl b/src/ForceFreeStates/Sing.jl index 60bf920b7..50c93f7de 100644 --- a/src/ForceFreeStates/Sing.jl +++ b/src/ForceFreeStates/Sing.jl @@ -115,7 +115,7 @@ function sing_lim!(intr::ForceFreeStatesInternal, ctrl::ForceFreeStatesControl, # Initial guesses based on equilibrium intr.qlim = min(equil.params.qmax, ctrl.qhigh) # equilibrium solve only goes up to qmax, so we're capped there intr.q1lim = profiles.q_deriv(profiles.xs[end]; hint=Ref(profiles.npts_minus_1)) - intr.psilim = equil.config.psihigh + intr.psilim = equil.params.psihigh_resolved # Optionally override qlim based on dmlim (Fortran sas_flag=t equivalent). # Multi-n runs (nn_low != nn_high) are not supported — the "outermost rational + dmlim/n" diff --git a/src/GeneralizedPerturbedEquilibrium.jl b/src/GeneralizedPerturbedEquilibrium.jl index 6d7266f3b..3dcebaa2d 100755 --- a/src/GeneralizedPerturbedEquilibrium.jl +++ b/src/GeneralizedPerturbedEquilibrium.jl @@ -318,12 +318,11 @@ function main_from_inputs( sing_lim!(intr, ctrl, equil) # If truncating before psihigh, reform equilibrium if desired - if intr.psilim != equil.config.psihigh && ctrl.reform_eq_with_psilim - @warn "Reforming equilibrium splines from psihigh to psilim not implemented yet. Proceeding with psihigh = $(equil.config.psihigh)." + if intr.psilim != equil.params.psihigh_resolved && ctrl.reform_eq_with_psilim + @warn "Reforming equilibrium splines from psihigh to psilim not implemented yet. Proceeding with psihigh = $(equil.params.psihigh_resolved)." # JMH - Nik please put the logic we discussed here # something like ? - # equil.config.psihigh = intr.psilim - # equil = set_up_equilibrium(equil.config) + # re-form the equilibrium with its outer limit set to intr.psilim end # Compute local stability (if desired). `locstab` holds `D_I` from the ballooning diff --git a/src/Rerun.jl b/src/Rerun.jl index 2a5a81237..b66d4e3ff 100644 --- a/src/Rerun.jl +++ b/src/Rerun.jl @@ -267,9 +267,10 @@ function build_inputs_from_h5(args::Vector{String}) " output: $(abspath(joinpath(output_dir, output_name)))\n$_BANNER" _drop_deprecated_keys!(inputs["Equilibrium"], _DEPRECATED_EQUIL_KEYS, "Equilibrium") - eq_config = Equilibrium.EquilibriumConfig(inputs["Equilibrium"], output_dir) - # Clear eq_filename: unused on replay, and a stale absolute path could mislead downstream code. - eq_config.eq_filename = "" + # Clear eq_filename on a copy: unused on replay, a stale absolute path could mislead + # downstream code, and `inputs` itself is re-serialized into the rerun's gpec_toml_raw. + equil_dict = merge(inputs["Equilibrium"], Dict{String,Any}("eq_filename" => "")) + eq_config = Equilibrium.EquilibriumConfig(equil_dict, output_dir) # Analytic kinds regenerate from their TOML section; file-based kinds rebuild splines from # the stored ingest. A file-based run with no ingest can only come from a pre-ingest gpec.h5. diff --git a/test/runtests_equil.jl b/test/runtests_equil.jl index dd1e49ce5..c69e1bd82 100644 --- a/test/runtests_equil.jl +++ b/test/runtests_equil.jl @@ -68,6 +68,31 @@ @test all(>(0), B_nodes) end + @testset "Resolved psihigh" begin + # The config holds the user's request and is never written to; the value the + # equilibrium is actually formed on rides on params.psihigh_resolved. + for eq in (plasma_eq_efit, plasma_eq_arclength, plasma_eq_inversion) + @test eq.params.psihigh_resolved == eq.rzphi_xs[end] + @test eq.params.psihigh_resolved <= eq.config.psihigh + end + # 0.994 sits inside the closed-flux region, so nothing is clamped + @test plasma_eq_efit.params.psihigh_resolved == 0.994 + + # Requesting the separatrix itself must leave the request intact on the config + edge_config = GeneralizedPerturbedEquilibrium.Equilibrium.EquilibriumConfig(; + eq_filename=joinpath(data_dir, "EQDSK_COCOS_02"), + eq_type="efit", + jac_type="boozer", + grid_type="ldp", + psilow=0.01, + psihigh=1.0 + ) + plasma_eq_edge = GeneralizedPerturbedEquilibrium.Equilibrium.setup_equilibrium(edge_config) + @test edge_config.psihigh == 1.0 + @test plasma_eq_edge.params.psihigh_resolved <= 1.0 + @test plasma_eq_edge.params.psihigh_resolved == plasma_eq_edge.rzphi_xs[end] + end + @testset "EFIT Method Consistency" begin # All three methods solve the same equilibrium — q-profiles should broadly agree. # Tolerance is 10% to allow for method-specific discretisation differences. diff --git a/test/runtests_sing.jl b/test/runtests_sing.jl index ff66bc035..90e521b9f 100644 --- a/test/runtests_sing.jl +++ b/test/runtests_sing.jl @@ -175,12 +175,12 @@ using FastInterpolations: cubic_interp, CubicFit, LinearBinarySearch, Series, Ex intr = GeneralizedPerturbedEquilibrium.ForceFreeStates.ForceFreeStatesInternal() GeneralizedPerturbedEquilibrium.ForceFreeStates.sing_lim!(intr, ctrl, equil) @test isapprox(intr.qlim, equil.params.qmax; atol=1e-12) - @test isapprox(intr.psilim, equil.config.psihigh; atol=1e-12) + @test isapprox(intr.psilim, equil.params.psihigh_resolved; atol=1e-12) ctrl = GeneralizedPerturbedEquilibrium.ForceFreeStates.ForceFreeStatesControl(; qhigh=max(equil.params.qmin + 0.1, equil.params.qmax - 0.5)) GeneralizedPerturbedEquilibrium.ForceFreeStates.sing_lim!(intr, ctrl, equil) @test intr.qlim < equil.params.qmax + 1e-12 - @test intr.psilim <= equil.config.psihigh + @test intr.psilim <= equil.params.psihigh_resolved q_at_psilim = equil.profiles.q_spline(intr.psilim) @test isapprox(q_at_psilim, intr.qlim; atol=1e-6) end From c4aba4eab4b4ea414a6d44b9dadd055e99f8cda1 Mon Sep 17 00:00:00 2001 From: Jake Halpern Date: Fri, 14 Aug 2026 14:45:21 -0400 Subject: [PATCH 3/9] FFS - CLEANUP - Remove the never-implemented reform_eq_with_psilim control The flag only ever gated a warning saying the reform was not implemented, so it could not change any result. Drops the field, its docstring entry, and the dead branch in the driver. Co-Authored-By: Claude Opus 5 --- src/ForceFreeStates/ForceFreeStatesStructs.jl | 2 -- src/GeneralizedPerturbedEquilibrium.jl | 8 -------- 2 files changed, 10 deletions(-) diff --git a/src/ForceFreeStates/ForceFreeStatesStructs.jl b/src/ForceFreeStates/ForceFreeStatesStructs.jl index 7c2a70d87..6b139d2df 100644 --- a/src/ForceFreeStates/ForceFreeStatesStructs.jl +++ b/src/ForceFreeStates/ForceFreeStatesStructs.jl @@ -244,7 +244,6 @@ gpec.toml. - `kinetic_source::String` - Kinetic matrix source: "fixed" (X-shaped test matrices scaled by kinetic_factor relative to ideal matrix Frobenius norms; Ak, Dk, Hk Hermitian, Bk, Ck, Ek non-Hermitian), "calculated" (PENTRC — not yet implemented) - `kinetic_factor::Float64` - Dimensionless scaling factor for kinetic matrices. Zero (the default) disables the kinetic path; any positive value enables it and scales the kinetic matrices: when kinetic_source="fixed", scales X-shaped test matrices relative to ideal matrix norms; when kinetic_source="calculated", applied as uniform post-hoc multiplier to W and T components. - `qlow::Float64` - Integration terminated at q limit determined by minimum of qlow and q0 from equil - - `reform_eq_with_psilim::Bool` - Reform equilibrium with computed psilim (not yet implemented) - `psiedge::Float64` - If less than psilim, records a dW(ψ) diagnostic scan over [psiedge, psilim] on odet.edge_scan. The integration domain (psilim) is always controlled by qhigh / psihigh and is not modified by this scan (unless `truncate_at_dW_peak=true`, see caveats below). - `truncate_at_dW_peak::Bool` - When `true` and `psiedge < psilim`, the edge-dW scan's peak location is adopted as the new physical plasma edge — `intr.psilim`/`intr.qlim`/`odet.u` are pulled back to the peak, AND the FM Δ' chunks/propagators are made self-consistent with the new boundary (the chunk that straddles the peak is rebuilt + re-integrated; any chunks past the peak are dropped). This reproduces the spirit of the original ode_record_edge heuristic from Fortran STRIDE while keeping Δ' and δW well-defined at the new boundary. The Δ' metric is still physically dependent on where the peak falls in the edge band, so use this flag deliberately when you mean to scan against the peak-defined edge (e.g. for studying edge-mode regimes); leave at `false` (default) for the full-domain Δ' at `qhigh` / `psihigh` / `dmlim`. - `diagnose::Bool` - Enable diagnostic output (not yet implemented) @@ -284,7 +283,6 @@ gpec.toml. kinetic_source::String = "fixed" kinetic_factor::Float64 = 0.0 qlow::Float64 = 0.0 - reform_eq_with_psilim::Bool = false psiedge::Float64 = 0.99 truncate_at_dW_peak::Bool = false # Edge-dW peak becomes new physical edge; Δ' BVP made self-consistent. See docstring. parallel_threads::Int = 2 diff --git a/src/GeneralizedPerturbedEquilibrium.jl b/src/GeneralizedPerturbedEquilibrium.jl index a2b514b11..e2053eff5 100755 --- a/src/GeneralizedPerturbedEquilibrium.jl +++ b/src/GeneralizedPerturbedEquilibrium.jl @@ -321,14 +321,6 @@ function main_from_inputs( # Determine psilim and qlim (where we will integrate to) sing_lim!(intr, ctrl, equil) - # If truncating before psihigh, reform equilibrium if desired - if intr.psilim != equil.params.psihigh_resolved && ctrl.reform_eq_with_psilim - @warn "Reforming equilibrium splines from psihigh to psilim not implemented yet. Proceeding with psihigh = $(equil.params.psihigh_resolved)." - # JMH - Nik please put the logic we discussed here - # something like ? - # re-form the equilibrium with its outer limit set to intr.psilim - end - # Compute local stability (if desired). `locstab` holds `D_I` from the ballooning # coefficient system and the local ballooning result; `nothing` when not computed. locstab = nothing From dd7e96eec349b8b3bfbd3fd12ad2a39e5d8db718 Mon Sep 17 00:00:00 2001 From: logan-nc Date: Sat, 15 Aug 2026 14:18:59 -0400 Subject: [PATCH 4/9] EQUIL - BUGFIX - Make newq0 a Float64 target resolved into a local newq0 is a target on-axis safety factor, not an index: the solver forms f0fac = f0^2*((newq0/q0)^2 - 1) from it. Declaring it ::Int admitted only integers plus the -1 sentinel and threw InexactError on any realistic request such as 1.05. The direct solver also wrote the resolved sentinel back into the config (equil_params.newq0 = -q0, where equil_params is raw_profile.config), which the immutability change turns into a setfield! error. Resolve into a local instead, matching equilibrium_solver(::InverseRunInput). No result can move: a non-integer newq0 threw at construction, so no working run used one, and TOML's `newq0 = 0` converts to 0.0 unchanged. Also note on DirectRunInput that IMAS equilibria are not in EFIT_KINDS and so are never clamped, since the psihigh_resolved name implies otherwise. Co-Authored-By: Claude Opus 5 (1M context) Claude-Session: https://claude.ai/code/session_01Dh1NCejnd3fYMmcRKoQRcG --- examples/DIIID-like_SLAYER_example/gpec.toml | 2 +- .../DIIID-like_gal_resistive_example/gpec.toml | 2 +- .../DIIID-like_gal_resistive_pe_example/gpec.toml | 2 +- examples/DIIID-like_ideal_example/gpec.toml | 2 +- examples/Solovev_ideal_example/gpec.toml | 2 +- examples/Solovev_ideal_example_3D/gpec.toml | 2 +- examples/Solovev_ideal_example_multi_n/gpec.toml | 2 +- .../single_n_1/gpec.toml | 2 +- .../single_n_2/gpec.toml | 2 +- examples/Solovev_kinetic_NTV_example/gpec.toml | 2 +- .../Solovev_kinetic_calculated_example/gpec.toml | 2 +- examples/a10_kinetic_example/gpec.toml | 2 +- src/Equilibrium/DirectEquilibrium.jl | 15 +++++++++------ src/Equilibrium/EquilibriumTypes.jl | 9 ++++++--- .../regression_solovev_ideal_example/gpec.toml | 2 +- .../gpec.toml | 2 +- .../gpec.toml | 2 +- .../regression_solovev_kinetic_example/gpec.toml | 2 +- .../regression_solovev_kinetic_multi_n/gpec.toml | 2 +- .../regression_solovev_kinetic_nuzero/gpec.toml | 2 +- 20 files changed, 33 insertions(+), 27 deletions(-) diff --git a/examples/DIIID-like_SLAYER_example/gpec.toml b/examples/DIIID-like_SLAYER_example/gpec.toml index a89b382f6..fd4060acd 100644 --- a/examples/DIIID-like_SLAYER_example/gpec.toml +++ b/examples/DIIID-like_SLAYER_example/gpec.toml @@ -14,7 +14,7 @@ psihigh = 0.9995 # Upper limit of normalized flux coordinate mpsi = 0 # Number of radial grid points (0 = auto-compute from psi_accuracy) psi_accuracy = 0.001 # Target absolute error in q for auto-mpsi mtheta = 256 # Number of poloidal grid points -newq0 = 0 # Override for on-axis safety factor (0 = use input value) +newq0 = 0 # Target on-axis safety factor q(0) (0 = use input value, -1 = flip axis extrapolation) etol = 1e-10 # Error tolerance for equilibrium solver force_termination = false # Terminate after equilibrium setup (skip stability calculations) diff --git a/examples/DIIID-like_gal_resistive_example/gpec.toml b/examples/DIIID-like_gal_resistive_example/gpec.toml index f84d8c52a..eb70d8792 100644 --- a/examples/DIIID-like_gal_resistive_example/gpec.toml +++ b/examples/DIIID-like_gal_resistive_example/gpec.toml @@ -14,7 +14,7 @@ psihigh = 0.993 # Upper limit of normalized flux coordinate (0.99 mpsi = 128 # Number of radial grid intervals (0 = two-pass auto grid from psi_accuracy) psi_accuracy = 0.001 # Target relative accuracy of splined profile derivatives for the auto grid mtheta = 256 # Number of poloidal grid points -newq0 = 0 # Override for on-axis safety factor (0 = use input value) +newq0 = 0 # Target on-axis safety factor q(0) (0 = use input value, -1 = flip axis extrapolation) etol = 1e-10 # Error tolerance for equilibrium solver force_termination = false # Terminate after equilibrium setup (skip stability calculations) diff --git a/examples/DIIID-like_gal_resistive_pe_example/gpec.toml b/examples/DIIID-like_gal_resistive_pe_example/gpec.toml index e5536c75e..12de94489 100644 --- a/examples/DIIID-like_gal_resistive_pe_example/gpec.toml +++ b/examples/DIIID-like_gal_resistive_pe_example/gpec.toml @@ -13,7 +13,7 @@ psihigh = 0.993 # Upper limit of normalized flux coordinate (0.99 mpsi = 128 # Number of radial grid intervals (0 = two-pass auto grid from psi_accuracy) psi_accuracy = 0.001 # Target relative accuracy of splined profile derivatives for the auto grid mtheta = 256 # Number of poloidal grid points -newq0 = 0 # Override for on-axis safety factor (0 = use input value) +newq0 = 0 # Target on-axis safety factor q(0) (0 = use input value, -1 = flip axis extrapolation) etol = 1e-10 # Error tolerance for equilibrium solver force_termination = false # Terminate after equilibrium setup (skip stability calculations) diff --git a/examples/DIIID-like_ideal_example/gpec.toml b/examples/DIIID-like_ideal_example/gpec.toml index 9a3bf9bc5..1aa48d680 100644 --- a/examples/DIIID-like_ideal_example/gpec.toml +++ b/examples/DIIID-like_ideal_example/gpec.toml @@ -13,7 +13,7 @@ psihigh = 0.995 # Upper limit of normalized poloidal flux (captur mpsi = 0 # Number of radial grid intervals (0 = two-pass auto grid from psi_accuracy) psi_accuracy = 0.001 # Target relative accuracy of splined profile derivatives for the auto grid mtheta = 256 # Number of poloidal grid points -newq0 = 0 # Override for on-axis safety factor (0 = use input value) +newq0 = 0 # Target on-axis safety factor q(0) (0 = use input value, -1 = flip axis extrapolation) etol = 1e-10 # Error tolerance for equilibrium solver force_termination = false # Terminate after equilibrium setup (skip stability calculations) diff --git a/examples/Solovev_ideal_example/gpec.toml b/examples/Solovev_ideal_example/gpec.toml index 66ba0d48b..3a101ce3d 100644 --- a/examples/Solovev_ideal_example/gpec.toml +++ b/examples/Solovev_ideal_example/gpec.toml @@ -11,7 +11,7 @@ psilow = 1e-4 # Lower limit of normalized poloidal flu psihigh = 0.9995 # Upper limit of normalized poloidal flux mpsi = 128 # Number of radial grid intervals (0 = two-pass auto grid from psi_accuracy) mtheta = 256 # Number of poloidal grid points -newq0 = 0 # Override for on-axis safety factor (0 = use input value) +newq0 = 0 # Target on-axis safety factor q(0) (0 = use input value, -1 = flip axis extrapolation) etol = 1e-7 # Error tolerance for equilibrium solver force_termination = false # Terminate after equilibrium setup (skip stability calculations) diff --git a/examples/Solovev_ideal_example_3D/gpec.toml b/examples/Solovev_ideal_example_3D/gpec.toml index 70123040b..7e530d042 100644 --- a/examples/Solovev_ideal_example_3D/gpec.toml +++ b/examples/Solovev_ideal_example_3D/gpec.toml @@ -10,7 +10,7 @@ psilow = 1e-4 # Lower limit of normalized poloidal flu psihigh = 0.9995 # Upper limit of normalized poloidal flux mpsi = 128 # Number of radial grid intervals (0 = two-pass auto grid from psi_accuracy) mtheta = 256 # Number of poloidal grid points -newq0 = 0 # Override for on-axis safety factor (0 = use input value) +newq0 = 0 # Target on-axis safety factor q(0) (0 = use input value, -1 = flip axis extrapolation) etol = 1e-7 # Error tolerance for equilibrium solver force_termination = false # Terminate after equilibrium setup (skip stability calculations) diff --git a/examples/Solovev_ideal_example_multi_n/gpec.toml b/examples/Solovev_ideal_example_multi_n/gpec.toml index 3fb3f7b06..592e74676 100644 --- a/examples/Solovev_ideal_example_multi_n/gpec.toml +++ b/examples/Solovev_ideal_example_multi_n/gpec.toml @@ -11,7 +11,7 @@ psilow = 1e-4 # Lower limit of normalized poloidal flu psihigh = 0.9995 # Upper limit of normalized poloidal flux mpsi = 128 # Number of radial grid intervals (0 = two-pass auto grid from psi_accuracy) mtheta = 256 # Number of poloidal grid points -newq0 = 0 # Override for on-axis safety factor (0 = use input value) +newq0 = 0 # Target on-axis safety factor q(0) (0 = use input value, -1 = flip axis extrapolation) etol = 1e-7 # Error tolerance for equilibrium solver force_termination = false # Terminate after equilibrium setup (skip stability calculations) diff --git a/examples/Solovev_ideal_example_multi_n/single_n_1/gpec.toml b/examples/Solovev_ideal_example_multi_n/single_n_1/gpec.toml index 8b2b75007..59760b231 100644 --- a/examples/Solovev_ideal_example_multi_n/single_n_1/gpec.toml +++ b/examples/Solovev_ideal_example_multi_n/single_n_1/gpec.toml @@ -10,7 +10,7 @@ psilow = 1e-4 # Lower limit of normalized poloidal flu psihigh = 0.9995 # Upper limit of normalized poloidal flux mpsi = 128 # Number of radial grid intervals (0 = two-pass auto grid from psi_accuracy) mtheta = 256 # Number of poloidal grid points -newq0 = 0 # Override for on-axis safety factor (0 = use input value) +newq0 = 0 # Target on-axis safety factor q(0) (0 = use input value, -1 = flip axis extrapolation) etol = 1e-7 # Error tolerance for equilibrium solver force_termination = false # Terminate after equilibrium setup (skip stability calculations) diff --git a/examples/Solovev_ideal_example_multi_n/single_n_2/gpec.toml b/examples/Solovev_ideal_example_multi_n/single_n_2/gpec.toml index c7a51c723..a80ac7741 100644 --- a/examples/Solovev_ideal_example_multi_n/single_n_2/gpec.toml +++ b/examples/Solovev_ideal_example_multi_n/single_n_2/gpec.toml @@ -10,7 +10,7 @@ psilow = 1e-4 # Lower limit of normalized poloidal flu psihigh = 0.9995 # Upper limit of normalized poloidal flux mpsi = 128 # Number of radial grid intervals (0 = two-pass auto grid from psi_accuracy) mtheta = 256 # Number of poloidal grid points -newq0 = 0 # Override for on-axis safety factor (0 = use input value) +newq0 = 0 # Target on-axis safety factor q(0) (0 = use input value, -1 = flip axis extrapolation) etol = 1e-7 # Error tolerance for equilibrium solver force_termination = false # Terminate after equilibrium setup (skip stability calculations) diff --git a/examples/Solovev_kinetic_NTV_example/gpec.toml b/examples/Solovev_kinetic_NTV_example/gpec.toml index 733a92c31..8458b627f 100644 --- a/examples/Solovev_kinetic_NTV_example/gpec.toml +++ b/examples/Solovev_kinetic_NTV_example/gpec.toml @@ -13,7 +13,7 @@ psilow = 1e-4 # Lower limit of normalized poloidal flu psihigh = 0.9995 # Upper limit of normalized poloidal flux mpsi = 128 # Number of radial grid points (0 = auto-compute from psi_accuracy) mtheta = 256 # Number of poloidal grid points -newq0 = 0 # Override for on-axis safety factor (0 = use input value) +newq0 = 0 # Target on-axis safety factor q(0) (0 = use input value, -1 = flip axis extrapolation) etol = 1e-7 # Error tolerance for equilibrium solver force_termination = false # Terminate after equilibrium setup (skip stability calculations) diff --git a/examples/Solovev_kinetic_calculated_example/gpec.toml b/examples/Solovev_kinetic_calculated_example/gpec.toml index 149f40af7..fa713eac7 100644 --- a/examples/Solovev_kinetic_calculated_example/gpec.toml +++ b/examples/Solovev_kinetic_calculated_example/gpec.toml @@ -11,7 +11,7 @@ psilow = 1e-4 # Lower limit of normalized flux coordin psihigh = 0.9995 # Upper limit of normalized flux coordinate mpsi = 16 # Number of radial grid points mtheta = 256 # Number of poloidal grid points -newq0 = 0 # Override for on-axis safety factor (0 = use input value) +newq0 = 0 # Target on-axis safety factor q(0) (0 = use input value, -1 = flip axis extrapolation) etol = 1e-7 # Error tolerance for equilibrium solver force_termination = false # Terminate after equilibrium setup (skip stability calculations) diff --git a/examples/a10_kinetic_example/gpec.toml b/examples/a10_kinetic_example/gpec.toml index d3f61b391..61c9d3b2e 100644 --- a/examples/a10_kinetic_example/gpec.toml +++ b/examples/a10_kinetic_example/gpec.toml @@ -15,7 +15,7 @@ psilow = 1e-3 # Lower limit of normalized poloidal flux psihigh = 0.99 # Upper limit of normalized poloidal flux mpsi = 16 # Number of radial grid points (low value for fast iteration) mtheta = 256 # Number of poloidal grid points -newq0 = 0 # Override for on-axis safety factor (0 = use input value) +newq0 = 0 # Target on-axis safety factor q(0) (0 = use input value, -1 = flip axis extrapolation) etol = 1e-7 # Error tolerance for equilibrium solver [Wall] diff --git a/src/Equilibrium/DirectEquilibrium.jl b/src/Equilibrium/DirectEquilibrium.jl index e0117c02e..961a0f7f5 100644 --- a/src/Equilibrium/DirectEquilibrium.jl +++ b/src/Equilibrium/DirectEquilibrium.jl @@ -545,15 +545,18 @@ robustness. if q0 <= 0.0 @warn "q0 extrapolation to axis gives q0 = $(@sprintf("%.3f", q0)) ≤ 0 — likely a spline artifact from psilow being too large; check psilow or use newq0 to override." end - if equil_params.newq0 == -1 - equil_params.newq0 = -q0 + # The -1 sentinel means "flip the extrapolated q0"; resolve it into a local so the + # config stays the user's request (matches equilibrium_solver(::InverseRunInput)). + newq0 = equil_params.newq0 + if newq0 == -1 + newq0 = -q0 end - if equil_params.newq0 != 0.0 - @info "Revising q-profile for newq0 = $(@sprintf("%.3f", equil_params.newq0))" + if newq0 != 0.0 + @info "Revising q-profile for newq0 = $(@sprintf("%.3f", newq0))" f0 = profiles.F_spline.y[1] - profiles.F_deriv(psi_nodes[1]; hint=Ref(1)) * psi_nodes[1] - f0fac = f0^2 * ((equil_params.newq0 / q0)^2 - 1.0) + f0fac = f0^2 * ((newq0 / q0)^2 - 1.0) for i in 1:(mpsi+1) - ffac = sqrt(1.0 + f0fac / profiles.F_spline.y[i]^2) * sign(equil_params.newq0) + ffac = sqrt(1.0 + f0fac / profiles.F_spline.y[i]^2) * sign(newq0) sq_nodes[i, 1] *= ffac sq_nodes[i, 4] *= ffac rzphi_nodes[i, :, 3] .*= ffac diff --git a/src/Equilibrium/EquilibriumTypes.jl b/src/Equilibrium/EquilibriumTypes.jl index 5ec2a82ce..a5ba8ed3c 100644 --- a/src/Equilibrium/EquilibriumTypes.jl +++ b/src/Equilibrium/EquilibriumTypes.jl @@ -39,7 +39,8 @@ specified in the input. - `psi_accuracy::Float64` - Target relative accuracy τ of splined profile derivatives for the two-pass auto grid (knot count scales as τ^(-1/3)) - `mtheta::Int` - Number of poloidal grid points - - `newq0::Int` - Override for on-axis safety factor (0 = use input value) + - `newq0::Float64` - Target on-axis safety factor q(0); the q and F profiles are rescaled to + meet it (0 = use input value, -1 = use the axis extrapolation with its sign flipped) - `etol::Float64` - Error tolerance for equilibrium solver - `force_termination::Bool` - Terminate after equilibrium setup (skip stability calculations) - `use_galgrid::Bool` - Use the same grid as galerkin method @@ -68,7 +69,7 @@ specified in the input. psi_accuracy::Float64 = 0.001 mtheta::Int = 512 - newq0::Int = 0 + newq0::Float64 = 0.0 etol::Float64 = 1e-10 force_termination::Bool = false @@ -445,7 +446,9 @@ raw equilibrium data and preparing the initial splines. - `psihigh_resolved::Float64` — outer flux limit the equilibrium is formed on: `config.psihigh` clamped to the outermost closed flux surface by [`resolve_psihigh!`](@ref). Defaults to `config.psihigh` and only differs for efit-family equilibria whose requested limit falls - outside the closed-flux region. The solvers build their ψ grid from this field. + outside the closed-flux region. The solvers build their ψ grid from this field. IMAS + equilibria are read into this struct but are not in `EFIT_KINDS`, so they are never + clamped and always keep the request. """ mutable struct DirectRunInput{S<:FastInterpolations.CubicSeriesInterpolant,I2D<:FastInterpolations.CubicInterpolantND} config::EquilibriumConfig diff --git a/test/test_data/regression_solovev_ideal_example/gpec.toml b/test/test_data/regression_solovev_ideal_example/gpec.toml index f2fc6d1f0..ecf933058 100644 --- a/test/test_data/regression_solovev_ideal_example/gpec.toml +++ b/test/test_data/regression_solovev_ideal_example/gpec.toml @@ -10,7 +10,7 @@ psilow = 1e-4 # Lower limit of normalized poloidal flu psihigh = 0.9995 # Upper limit of normalized poloidal flux mpsi = 16 # Number of radial grid intervals (0 = two-pass auto grid from psi_accuracy) mtheta = 256 # Number of poloidal grid points -newq0 = 0 # Override for on-axis safety factor (0 = use input value) +newq0 = 0 # Target on-axis safety factor q(0) (0 = use input value, -1 = flip axis extrapolation) etol = 1e-7 # Error tolerance for equilibrium solver force_termination = false # Terminate after equilibrium setup (skip stability calculations) diff --git a/test/test_data/regression_solovev_ideal_example_multi_n/gpec.toml b/test/test_data/regression_solovev_ideal_example_multi_n/gpec.toml index 0f0bc5c47..56131c9a3 100644 --- a/test/test_data/regression_solovev_ideal_example_multi_n/gpec.toml +++ b/test/test_data/regression_solovev_ideal_example_multi_n/gpec.toml @@ -10,7 +10,7 @@ psilow = 1e-4 # Lower limit of normalized poloidal flu psihigh = 0.9995 # Upper limit of normalized poloidal flux mpsi = 16 # Number of radial grid intervals (0 = two-pass auto grid from psi_accuracy) mtheta = 256 # Number of poloidal grid points -newq0 = 0 # Override for on-axis safety factor (0 = use input value) +newq0 = 0 # Target on-axis safety factor q(0) (0 = use input value, -1 = flip axis extrapolation) etol = 1e-7 # Error tolerance for equilibrium solver force_termination = false # Terminate after equilibrium setup (skip stability calculations) diff --git a/test/test_data/regression_solovev_kinetic_calculated/gpec.toml b/test/test_data/regression_solovev_kinetic_calculated/gpec.toml index 5b87267a9..47e8985e6 100644 --- a/test/test_data/regression_solovev_kinetic_calculated/gpec.toml +++ b/test/test_data/regression_solovev_kinetic_calculated/gpec.toml @@ -10,7 +10,7 @@ psilow = 1e-4 # Lower limit of normalized poloidal flu psihigh = 0.9995 # Upper limit of normalized poloidal flux mpsi = 16 # Number of radial grid intervals (0 = two-pass auto grid from psi_accuracy) mtheta = 256 # Number of poloidal grid points -newq0 = 0 # Override for on-axis safety factor (0 = use input value) +newq0 = 0 # Target on-axis safety factor q(0) (0 = use input value, -1 = flip axis extrapolation) etol = 1e-7 # Error tolerance for equilibrium solver force_termination = false # Terminate after equilibrium setup (skip stability calculations) diff --git a/test/test_data/regression_solovev_kinetic_example/gpec.toml b/test/test_data/regression_solovev_kinetic_example/gpec.toml index cef6ad546..4872c79e1 100644 --- a/test/test_data/regression_solovev_kinetic_example/gpec.toml +++ b/test/test_data/regression_solovev_kinetic_example/gpec.toml @@ -10,7 +10,7 @@ psilow = 1e-4 # Lower limit of normalized poloidal flu psihigh = 0.9995 # Upper limit of normalized poloidal flux mpsi = 16 # Number of radial grid intervals (0 = two-pass auto grid from psi_accuracy) mtheta = 256 # Number of poloidal grid points -newq0 = 0 # Override for on-axis safety factor (0 = use input value) +newq0 = 0 # Target on-axis safety factor q(0) (0 = use input value, -1 = flip axis extrapolation) etol = 1e-7 # Error tolerance for equilibrium solver force_termination = false # Terminate after equilibrium setup (skip stability calculations) diff --git a/test/test_data/regression_solovev_kinetic_multi_n/gpec.toml b/test/test_data/regression_solovev_kinetic_multi_n/gpec.toml index b11caab8b..20522dd36 100644 --- a/test/test_data/regression_solovev_kinetic_multi_n/gpec.toml +++ b/test/test_data/regression_solovev_kinetic_multi_n/gpec.toml @@ -10,7 +10,7 @@ psilow = 1e-4 # Lower limit of normalized poloidal flu psihigh = 0.9995 # Upper limit of normalized poloidal flux mpsi = 16 # Number of radial grid intervals (0 = two-pass auto grid from psi_accuracy) mtheta = 256 # Number of poloidal grid points -newq0 = 0 # Override for on-axis safety factor (0 = use input value) +newq0 = 0 # Target on-axis safety factor q(0) (0 = use input value, -1 = flip axis extrapolation) etol = 1e-7 # Error tolerance for equilibrium solver force_termination = false # Terminate after equilibrium setup (skip stability calculations) diff --git a/test/test_data/regression_solovev_kinetic_nuzero/gpec.toml b/test/test_data/regression_solovev_kinetic_nuzero/gpec.toml index ca6eeb003..8053523f1 100644 --- a/test/test_data/regression_solovev_kinetic_nuzero/gpec.toml +++ b/test/test_data/regression_solovev_kinetic_nuzero/gpec.toml @@ -11,7 +11,7 @@ psilow = 1e-4 # Lower limit of normalized flux coordin psihigh = 0.9995 # Upper limit of normalized flux coordinate mpsi = 16 # Number of radial grid points mtheta = 256 # Number of poloidal grid points -newq0 = 0 # Override for on-axis safety factor (0 = use input value) +newq0 = 0 # Target on-axis safety factor q(0) (0 = use input value, -1 = flip axis extrapolation) etol = 1e-7 # Error tolerance for equilibrium solver force_termination = false # Terminate after equilibrium setup (skip stability calculations) From a97e6d4292e487ba9ea960a4d0672c6fc0280884 Mon Sep 17 00:00:00 2001 From: logan-nc Date: Sat, 15 Aug 2026 14:19:16 -0400 Subject: [PATCH 5/9] GPEC - BUGFIX - Keep reform_eq_with_psilim parsing as a deprecated key The control was removed from ForceFreeStatesControl, but the struct is built by splatting the whole [ForceFreeStates] table as kwargs, so any existing gpec.toml still carrying the key now dies with an unknown-keyword MethodError -- as does replaying an older gpec.h5, whose stored TOML blob goes through the same path. _drop_deprecated_keys! runs inside main_from_inputs, so one tuple entry covers both. Co-Authored-By: Claude Opus 5 (1M context) Claude-Session: https://claude.ai/code/session_01Dh1NCejnd3fYMmcRKoQRcG --- src/GeneralizedPerturbedEquilibrium.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/GeneralizedPerturbedEquilibrium.jl b/src/GeneralizedPerturbedEquilibrium.jl index 5b85835bd..5da38bd5f 100755 --- a/src/GeneralizedPerturbedEquilibrium.jl +++ b/src/GeneralizedPerturbedEquilibrium.jl @@ -71,7 +71,7 @@ using .ForceFreeStates: find_kinetic_singular_surfaces! using .ForceFreeStates: eulerlagrange_integration, free_run, normalize_eigenfunctions! using .ForceFreeStates: galerkin_solve, write_galerkin!, GalerkinResult, gal_matched_odestate -const _DEPRECATED_FFS_KEYS = ("mer_flag", "force_wv_symmetry", "ode_flag", "cyl_flag", "mat_flag") +const _DEPRECATED_FFS_KEYS = ("mer_flag", "force_wv_symmetry", "ode_flag", "cyl_flag", "mat_flag", "reform_eq_with_psilim") const _DEPRECATED_EQUIL_KEYS = ("power_bp", "power_b", "power_r", "power_rc") # Drop deprecated keys from a parsed gpec.toml section so legacy files keep parsing From f78afce240040a9075bdf6340461c7f7b8cbf51a Mon Sep 17 00:00:00 2001 From: logan-nc Date: Sat, 15 Aug 2026 14:19:16 -0400 Subject: [PATCH 6/9] H5 - MINOR - Annotate the new Equilibrium/psihigh_resolved dataset EquilibriumParameters is dumped generically by field name, so adding psihigh_resolved adds a gpec.h5 dataset automatically. Without a metadata entry the schema walk in h5_metadata_check.jl fails it for missing long_name and units. No EQUIL_H5_NAMES entry is needed -- the writer falls back to the field name, which reads consistently next to Info/psilim -- and annotate! defaults units to "1". Co-Authored-By: Claude Opus 5 (1M context) Claude-Session: https://claude.ai/code/session_01Dh1NCejnd3fYMmcRKoQRcG --- src/HDF5Schema.jl | 1 + 1 file changed, 1 insertion(+) diff --git a/src/HDF5Schema.jl b/src/HDF5Schema.jl index a8beacffa..a1a47ef44 100644 --- a/src/HDF5Schema.jl +++ b/src/HDF5Schema.jl @@ -65,6 +65,7 @@ const MAIN_H5_ANNOTATIONS = [ "Equilibrium/R_axis" => (; long_name="R-coordinate of the magnetic axis", units="m"), "Equilibrium/Z_axis" => (; long_name="Z-coordinate of the magnetic axis", units="m"), "Equilibrium/psi_total" => (; long_name="total poloidal flux difference |ψ_axis − ψ_boundary|", units="Wb/rad"), + "Equilibrium/psihigh_resolved" => (; long_name="normalized poloidal flux at the outermost formed flux surface (requested psihigh clamped to the outermost closed surface)"), "Equilibrium/R_midplane" => (; long_name="R of the boundary at the inboard and outboard midplane crossings", units="m"), "Equilibrium/R_extremum" => (; long_name="R of the boundary at its upper and lower Z extrema", units="m"), "Equilibrium/Z_extremum" => (; long_name="Z of the boundary at its upper and lower Z extrema", units="m"), From 31df659f86203de4f505cf0da46285041379f8ff Mon Sep 17 00:00:00 2001 From: logan-nc Date: Sat, 15 Aug 2026 14:19:16 -0400 Subject: [PATCH 7/9] BENCH - BUGFIX - Rebuild the config per mpsi instead of mutating a frozen copy The mpsi scan still did deepcopy(eq_config) + cfg.mpsi = N, which errors now that EquilibriumConfig is immutable; only reference_q() had been converted. Carry the equilibrium dict out of reference_q() and build a fresh config per N, as reference_q() already does. Not covered by CI. Co-Authored-By: Claude Opus 5 (1M context) Claude-Session: https://claude.ai/code/session_01Dh1NCejnd3fYMmcRKoQRcG --- benchmarks/benchmark_q_vs_iota_edge.jl | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/benchmarks/benchmark_q_vs_iota_edge.jl b/benchmarks/benchmark_q_vs_iota_edge.jl index 2dfb89eb6..828f1bbb0 100644 --- a/benchmarks/benchmark_q_vs_iota_edge.jl +++ b/benchmarks/benchmark_q_vs_iota_edge.jl @@ -24,7 +24,8 @@ function reference_q() equil_dict = merge(inputs["Equilibrium"], Dict{String,Any}("grid_type" => "ldp", "mpsi" => 1024)) eq_config = GPE.Equilibrium.EquilibriumConfig(equil_dict, EXAMPLE_DIR) equil = GPE.Equilibrium.setup_equilibrium(eq_config, additional_input) - return equil, eq_config + # The dict rides along so the mpsi scan can rebuild configs; EquilibriumConfig is immutable. + return equil, equil_dict end # Rational surface ψ_s(q = m/n) and q'(ψ_s) recovered from a fitted spline by bisection @@ -44,7 +45,7 @@ function rational_locations(q_itp, dq_itp, psis, q_targets) end function main() - equil, eq_config = reference_q() + equil, equil_dict = reference_q() xs_ref = equil.profiles.xs q_ref_itp = equil.profiles.q_spline dq_ref_itp = equil.profiles.q_deriv @@ -64,8 +65,7 @@ function main() for N in (16, 32, 64, 128, 256) # Same log_asymptotic knot layout the auto grid uses for fixed mpsi - cfg = deepcopy(eq_config) - cfg.mpsi = N + cfg = GPE.Equilibrium.EquilibriumConfig(merge(equil_dict, Dict{String,Any}("mpsi" => N)), EXAMPLE_DIR) knots = GPE.Equilibrium._build_psi_grid(cfg, psilow, psihigh) q_nodes = [q_ref_itp(p) for p in knots] From 2873a81b28c72c6a08857654fb482f0ebf337baf Mon Sep 17 00:00:00 2001 From: logan-nc Date: Sat, 15 Aug 2026 14:19:17 -0400 Subject: [PATCH 8/9] EQUIL - TESTS - Cover newq0, deprecated keys and psihigh_resolved on inverse readers newq0 != 0 had no coverage at all, which is why the broken override stayed green. The -1 sentinel makes f0fac vanish, so ffac is exactly -1 and the revised q-profile is the negated baseline -- asserted exactly. A non-integer target covers the type fix. Also assert psihigh_resolved reaches params on the CHEASE inverse path (where the separatrix clamp never runs), and that a deprecated ForceFreeStates key warns and is dropped rather than killing the run. Co-Authored-By: Claude Opus 5 (1M context) Claude-Session: https://claude.ai/code/session_01Dh1NCejnd3fYMmcRKoQRcG --- test/runtests_equil.jl | 56 ++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 56 insertions(+) diff --git a/test/runtests_equil.jl b/test/runtests_equil.jl index c69e1bd82..19f7e9e38 100644 --- a/test/runtests_equil.jl +++ b/test/runtests_equil.jl @@ -93,6 +93,53 @@ @test plasma_eq_edge.params.psihigh_resolved == plasma_eq_edge.rzphi_xs[end] end + @testset "newq0 q-profile revision" begin + # newq0 is a target q(0), not an index: non-integer requests must survive the config, + # and neither the -1 sentinel nor an explicit target may write back to the frozen config. + newq0_config(newq0) = GeneralizedPerturbedEquilibrium.Equilibrium.EquilibriumConfig(; + eq_filename=joinpath(data_dir, "EQDSK_COCOS_02"), + eq_type="efit", + jac_type="boozer", + grid_type="ldp", + mpsi=32, + psilow=0.01, + psihigh=0.994, + newq0=newq0 + ) + + base_config = newq0_config(0) + eq_base = GeneralizedPerturbedEquilibrium.Equilibrium.setup_equilibrium(base_config) + + # -1 means "flip the sign of the axis extrapolation": f0fac vanishes, so ffac is + # exactly -1 and the revised q-profile is the negated baseline. + flip_config = newq0_config(-1) + eq_flip = GeneralizedPerturbedEquilibrium.Equilibrium.setup_equilibrium(flip_config) + @test flip_config.newq0 == -1.0 + @test eq_flip.profiles.q_spline.y == -eq_base.profiles.q_spline.y + + # A non-integer target used to throw InexactError at construction (newq0 was ::Int) + target_config = newq0_config(1.05) + @test target_config.newq0 == 1.05 + eq_target = GeneralizedPerturbedEquilibrium.Equilibrium.setup_equilibrium(target_config) + @test all(isfinite, eq_target.profiles.q_spline.y) + # Extrapolated q(0) meets the target; not exact because ffac uses the first node's F + # rather than the extrapolated axis value. + xs = eq_target.profiles.xs + q_axis = eq_target.profiles.q_spline.y[1] - eq_target.profiles.q_deriv(xs[1]; hint=Ref(1)) * xs[1] + @test isapprox(q_axis, 1.05; rtol=0.02) + end + + @testset "Deprecated TOML keys are dropped, not fatal" begin + # Removed control knobs must keep old gpec.toml decks (and older gpec.h5 replays, + # whose stored TOML blob goes through the same path) parsing with a warning. + ffs_table = Dict{String,Any}("reform_eq_with_psilim" => false, "nn_low" => 1) + @test_logs (:warn, r"reform_eq_with_psilim") GeneralizedPerturbedEquilibrium._drop_deprecated_keys!( + ffs_table, GeneralizedPerturbedEquilibrium._DEPRECATED_FFS_KEYS, "ForceFreeStates") + @test !haskey(ffs_table, "reform_eq_with_psilim") + @test GeneralizedPerturbedEquilibrium.ForceFreeStates.ForceFreeStatesControl(; + (Symbol(k) => v for (k, v) in ffs_table)...) isa GeneralizedPerturbedEquilibrium.ForceFreeStates.ForceFreeStatesControl + end + @testset "EFIT Method Consistency" begin # All three methods solve the same equilibrium — q-profiles should broadly agree. # Tolerance is 10% to allow for method-specific discretisation differences. @@ -142,6 +189,15 @@ @test plasma_eq_ascii isa GeneralizedPerturbedEquilibrium.Equilibrium.PlasmaEquilibrium end + @testset "Resolved psihigh (inverse readers)" begin + # CHEASE data already conforms to the plasma boundary, so the separatrix clamp never + # runs and the resolved value is the request — but it must still reach params. + for eq in (plasma_eq_binary, plasma_eq_ascii) + @test eq.params.psihigh_resolved == 0.994 + @test isapprox(eq.params.psihigh_resolved, eq.rzphi_xs[end]; atol=1e-12) + end + end + @testset "CHEASE Consistency (ASCII vs Binary)" begin # Both formats encode the same physical data; differences arise only from # floating-point text serialization in the ASCII format vs exact binary storage. From a3bdfd21853738ce32f1132272a12c2a563a4c31 Mon Sep 17 00:00:00 2001 From: Matthew Pharr Date: Sun, 16 Aug 2026 00:13:16 -0400 Subject: [PATCH 9/9] EQUIL - MINOR - address Nik's docs comments --- src/Equilibrium/EquilibriumTypes.jl | 12 ++++-------- test/runtests_equil.jl | 9 --------- 2 files changed, 4 insertions(+), 17 deletions(-) diff --git a/src/Equilibrium/EquilibriumTypes.jl b/src/Equilibrium/EquilibriumTypes.jl index a5ba8ed3c..344f5a60f 100644 --- a/src/Equilibrium/EquilibriumTypes.jl +++ b/src/Equilibrium/EquilibriumTypes.jl @@ -27,11 +27,8 @@ specified in the input. refinement when mpsi=0, three-region log layout when mpsi>0; "ldp", "pow1", "uniform"; "log_asymptotic" is a legacy alias for "auto") - `psilow::Float64` - Lower limit of normalized flux coordinate - - `psihigh::Float64` - Requested upper limit of normalized flux coordinate. For efit-family - equilibria this is the user's request, which may lie outside the closed-flux region; the - value the equilibrium is actually formed on is `DirectRunInput.psihigh_resolved` (carried - onto the equilibrium as `EquilibriumParameters.psihigh_resolved`). Read that, not this, - for anything downstream of `setup_equilibrium`. + - `psihigh::Float64` - Requested upper limit of normalized flux coordinate; the value the + equilibrium is actually formed on is `EquilibriumParameters.psihigh_resolved`. - `mpsi::Int` - Number of radial grid intervals; 0 with grid_type="auto" selects the two-pass auto grid: the main driver forms a coarse pass-1 equilibrium, measures its curvature, pins knots on rational surfaces, and re-forms on the refined grid. Standalone `setup_equilibrium` @@ -524,9 +521,8 @@ A mutable struct containing computed equilibrium parameters and diagnostic flags - `ro::Union{Nothing,Float64}` - R-coordinate of the magnetic axis [m] - `zo::Union{Nothing,Float64}` - Z-coordinate of the magnetic axis [m] - `psio::Union{Nothing,Float64}` - Total flux difference |ψ_axis - ψ_boundary| [Wb/rad] - - `psihigh_resolved::Union{Nothing,Float64}` - Outer flux limit the equilibrium was formed on, - equal to the outermost ψ node. This is `config.psihigh` clamped to the outermost closed flux - surface; downstream code wanting the plasma edge must read this, not `config.psihigh`. + - `psihigh_resolved::Union{Nothing,Float64}` - Outer flux limit the equilibrium was formed on + (the outermost ψ node); the plasma edge downstream of `setup_equilibrium`. - `rsep::Union{Nothing,Vector{Float64}}` - R-coordinates of the plasma boundary [m] - `zsep::Union{Nothing,Vector{Float64}}` - Z-coordinates of the plasma boundary [m] - `rext::Union{Nothing,Vector{Float64}}` - R-coordinates of the plasma edge [m] diff --git a/test/runtests_equil.jl b/test/runtests_equil.jl index 19f7e9e38..af3cd5743 100644 --- a/test/runtests_equil.jl +++ b/test/runtests_equil.jl @@ -189,15 +189,6 @@ @test plasma_eq_ascii isa GeneralizedPerturbedEquilibrium.Equilibrium.PlasmaEquilibrium end - @testset "Resolved psihigh (inverse readers)" begin - # CHEASE data already conforms to the plasma boundary, so the separatrix clamp never - # runs and the resolved value is the request — but it must still reach params. - for eq in (plasma_eq_binary, plasma_eq_ascii) - @test eq.params.psihigh_resolved == 0.994 - @test isapprox(eq.params.psihigh_resolved, eq.rzphi_xs[end]; atol=1e-12) - end - end - @testset "CHEASE Consistency (ASCII vs Binary)" begin # Both formats encode the same physical data; differences arise only from # floating-point text serialization in the ASCII format vs exact binary storage.