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 <noreply@anthropic.com>

Author: jhalpern30

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

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,

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)
 

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"]

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]

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,

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"

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

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.

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.