GPEC - IMPROVEMENT - Make all input control structs immutable

.pre-commit-config.yaml

@@ -69,5 +69,5 @@ repos:
       - id: toml-no-deprecated-keys
         name: 'TOML conventions: no deprecated config keys'
         language: pygrep
-        entry: '^(mer_flag|force_wv_symmetry|ode_flag|cyl_flag|mat_flag|use_riccati|use_parallel|parallel_threads|populate_dense_xi|power_bp|power_b|power_r|power_rc)\s*='
+        entry: '^(mer_flag|force_wv_symmetry|ode_flag|cyl_flag|mat_flag|reform_eq_with_psilim|use_riccati|use_parallel|parallel_threads|populate_dense_xi|power_bp|power_b|power_r|power_rc)\s*='
         files: ^(examples/.*\.toml|test/test_data/.*\.toml)$

benchmarks/benchmark_q_vs_iota_edge.jl

@@ -20,11 +20,12 @@ 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
+    # 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]
 

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)
 

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)
 

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)
 

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)
 

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)
 

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)
 

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)
 

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)
 

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)
 

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)
 

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)
 

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]

src/Equilibrium/DirectEquilibrium.jl

@@ -478,7 +478,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)
@@ -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
@@ -651,6 +654,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"]