From aded8ea839a37c28bbcd4f6c0137bd2a9b70960c Mon Sep 17 00:00:00 2001 From: StuartBenjamin Date: Sun, 18 Oct 2026 03:38:52 -0700 Subject: [PATCH] Equilibrium - BUGFIX - End the original inverse grid at psilow and psihigh MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit grid_type = "original" ran from the first off-axis file node to the last file surface, ignoring psilow and psihigh. On an edge-packed i-file the domain then reached ψ_N 0.9999 instead of psihigh, so the dmlim edge moved (psilim 0.996 -> 1.000 on a DIII-D i-file) and every Δ′ changed. The grid is now psilow, the file's nodes strictly between psilow and psihigh, and psihigh, as in Fortran GPEC. Co-Authored-By: Claude Opus 5.5 --- src/Equilibrium/EquilibriumTypes.jl | 2 +- src/Equilibrium/InverseEquilibrium.jl | 4 ++-- test/runtests_equil.jl | 5 +++-- 3 files changed, 6 insertions(+), 5 deletions(-) diff --git a/src/Equilibrium/EquilibriumTypes.jl b/src/Equilibrium/EquilibriumTypes.jl index bde720254..6846904ca 100644 --- a/src/Equilibrium/EquilibriumTypes.jl +++ b/src/Equilibrium/EquilibriumTypes.jl @@ -26,7 +26,7 @@ specified in the input. - `grid_type::String` - Grid type for flux surface discretization ("auto" — two-pass measured-curvature refinement when mpsi=0, three-region log layout when mpsi>0; "ldp", "pow1", "uniform"; inverse equilibria also "rho" (uniform in √ψ from psilow to psihigh) and "original" (the input - file's own ψ nodes off the axis; ignores psilow and mpsi); "log_asymptotic" is a legacy alias for "auto") + file's own ψ nodes between psilow and psihigh; ignores mpsi); "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; the value the equilibrium is actually formed on is `EquilibriumParameters.psihigh_resolved`. diff --git a/src/Equilibrium/InverseEquilibrium.jl b/src/Equilibrium/InverseEquilibrium.jl index 6453ae38d..cd60f44c3 100644 --- a/src/Equilibrium/InverseEquilibrium.jl +++ b/src/Equilibrium/InverseEquilibrium.jl @@ -164,8 +164,8 @@ function equilibrium_solver(input::InverseRunInput; override_psi_nodes::Union{No # Uniform in √ψ from psilow to psihigh (Fortran's rho ignores psilow) sq_xs = [(sqrt(psilow) + (sqrt(psihigh) - sqrt(psilow)) * i / mpsi)^2 for i in 0:mpsi] elseif grid_type in ("original", "orig") - # The input's own ψ nodes off the axis; psilow and mpsi are ignored, psihigh still sets psilim - sq_xs = filter(>(0.0), rz_in_xs) + # The input's own ψ nodes, ended by psilow and psihigh as the other grids are; mpsi is ignored + sq_xs = [psilow; filter(x -> psilow < x < psihigh, rz_in_xs); psihigh] mpsi = length(sq_xs) - 1 elseif grid_type == "uniform" if mpsi == 0 diff --git a/test/runtests_equil.jl b/test/runtests_equil.jl index 9b56986f6..a7872eff6 100644 --- a/test/runtests_equil.jl +++ b/test/runtests_equil.jl @@ -181,9 +181,10 @@ splines(eq) = (eq.rzphi_rsquared, eq.rzphi_offset, eq.rzphi_nu, eq.rzphi_jac, eq.eqfun_B, eq.eqfun_metric1, eq.eqfun_metric2) finite(eq) = all(isfinite, eq.profiles.dVdpsi_spline.y) && all(all(isfinite, s.nodal_derivs.partials) for s in splines(eq)) - # original: the file's own ψ nodes off the axis (this file starts at the axis), whatever psilow/mpsi + # original: the file's own ψ nodes between psilow and psihigh, whatever mpsi eq_o = Eq.setup_equilibrium(cfg("original")) - @test eq_o.rzphi_xs == Eq.read_chease_binary(cfg("original")).rz_in_xs[2:end] + xs_in = Eq.read_chease_binary(cfg("original")).rz_in_xs + @test eq_o.rzphi_xs == [0.01; filter(x -> 0.01 < x < 0.99, xs_in); 0.99] @test finite(eq_o) # rho: starts at psilow (Fortran's rho starts at psihigh/(mpsi+1)²)