From 5f32a39ce2f69864378152dbee12a826f24b1d1b Mon Sep 17 00:00:00 2001 From: logan-nc <6198372+logan-nc@users.noreply.github.com> Date: Fri, 31 Jul 2026 17:55:07 -0700 Subject: [PATCH 1/4] =?UTF-8?q?EQUIL=20-=20BUGFIX=20-=20Fix=20ideal=20?= =?UTF-8?q?=CE=94=E2=80=B2=20non-convergence=20on=20the=20auto=20psi-grid?= MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit The two-pass auto grid pinned a knot on each rational surface; the Δ′ asymptotic matching stencil [ψ_s ± singfac_min/|n·q′|] then straddled it, so the cubic 3rd derivative the singular extraction samples was inconsistent across the crossing and the Δ′ diagonal oscillated with psi_accuracy (7.2–10.1 on the DIII-D-like case) instead of converging. - Bracket each rational (bracket_mandatory_nodes): center ψ_s inside a clean, locally-uniform interval between a knot pair, never on a knot, so the matching stencil lies in one spline segment and q‴ is single-valued across it. - Fixed-width (τ-independent) near-rational resolution patch so Δ′ is accurate at any psi_accuracy — Δ′ is a property of the rational surface, so its resolution must not scale with the global accuracy target (that only sizes the pedestal/edge grid). - Global minimum knot-spacing floor (enforce_min_spacing) on the auto refined grid and the fixed ldp/pow1 grids, so no grid packs its edge below the field-line-integration noise scale at high mpsi. dpm[1,1] now converges τ-independently to ~9.2 (matching the uniform mpsi ladder) at the default psi_accuracy=1e-3, vs the prior erratic 7.2–10.1. Equilibrium scalars are unchanged (grid-placement only). The accuracy-driven density remains grid-dependent near the separatrix (a real q/geometry singularity); that is tracked separately for a self-consistent / one-pass grid redesign and does not affect this fix. Co-Authored-By: Claude Opus 4.8 Claude-Session: https://claude.ai/code/session_0125dY2qbGtm3gky8TPv7vh6 --- src/Equilibrium/DirectEquilibrium.jl | 14 +- src/Equilibrium/Equilibrium.jl | 2 +- src/Equilibrium/GridRefinement.jl | 181 +++++++++++++++++++------ src/Equilibrium/InverseEquilibrium.jl | 15 +- src/GeneralizedPerturbedEquilibrium.jl | 8 +- test/runtests_grid_refinement.jl | 58 +++++++- 6 files changed, 221 insertions(+), 57 deletions(-) diff --git a/src/Equilibrium/DirectEquilibrium.jl b/src/Equilibrium/DirectEquilibrium.jl index 02f18410e..24b78f218 100644 --- a/src/Equilibrium/DirectEquilibrium.jl +++ b/src/Equilibrium/DirectEquilibrium.jl @@ -415,10 +415,11 @@ solvers. function _build_psi_grid(equil_params, psilow, psihigh) mpsi = equil_params.mpsi if equil_params.grid_type in ("auto", "log_asymptotic") && mpsi == 0 && equil_params.psi_accuracy > 0 - # Two-pass auto grid: this is the coarse pass-1 layout; the driver measures the - # formed equilibrium's curvature and re-forms on a refined grid (GridRefinement.jl). - mpsi = 128 - @info "Auto psi grid: forming pass-1 equilibrium on coarse $(mpsi)-interval log_asymptotic grid pending curvature-based refinement" + # Two-pass auto grid: this is the pass-1 layout; the driver measures the formed + # equilibrium's curvature and re-forms on a refined grid (GridRefinement.jl). Sized to + # resolve the gradient structure accurately so pass 2 is well-provisioned (PASS1_INTERVALS). + mpsi = PASS1_INTERVALS + @info "Auto psi grid: forming pass-1 equilibrium on $(mpsi)-interval log_asymptotic grid pending curvature-based refinement" elseif mpsi == 0 mpsi = 128 end @@ -444,7 +445,10 @@ function _build_psi_grid(equil_params, psilow, psihigh) else error("Unsupported grid_type: $(equil_params.grid_type)") end - return psi_nodes + # Floor node spacing on the fixed grids: ldp/pow1 pack the edge as ~(π/2·mpsi)⁻² and at high + # mpsi drive it below the integration-noise scale (garbage curvature). The auto grid floors its + # refined pass-2 grid in refined_psi_grid, so it is left untouched here. + return equil_params.grid_type in ("auto", "log_asymptotic") ? psi_nodes : enforce_min_spacing(psi_nodes, MIN_KNOT_SPACING) end """ diff --git a/src/Equilibrium/Equilibrium.jl b/src/Equilibrium/Equilibrium.jl index 03935640e..cf08d8a4f 100644 --- a/src/Equilibrium/Equilibrium.jl +++ b/src/Equilibrium/Equilibrium.jl @@ -30,7 +30,7 @@ export setup_equilibrium, EquilibriumConfig, PlasmaEquilibrium, EquilibriumParam KineticProfileSplines, load_kinetic_profiles, KineticProfileData, read_kinetic_file, write_kinetic_h5 export flux_surface_metric, flux_surface_area -export wants_two_pass, refined_psi_grid, merge_mandatory_nodes, implied_knot_count +export wants_two_pass, refined_psi_grid, merge_mandatory_nodes, bracket_mandatory_nodes, enforce_min_spacing, implied_knot_count export compute_sqrt_jac_delpsi, compute_sqrtamat, rootarea_to_area_weight, area_to_rootarea_weight # --- Constants --- diff --git a/src/Equilibrium/GridRefinement.jl b/src/Equilibrium/GridRefinement.jl index 9478d96bc..e58eb81dc 100644 --- a/src/Equilibrium/GridRefinement.jl +++ b/src/Equilibrium/GridRefinement.jl @@ -57,12 +57,42 @@ const CORE_MODEL_PSI_MAX = 0.03 const EDGE_MODEL_PSI_MIN = 0.9 # θ-lines subsample stride for the 2D geometry channels const THETA_STRIDE = 8 -# Local density elevation around mandatory (rational) surfaces: knot spacing h_s = coef·τ^(1/3) -# at the surface, geometric growth away from it, applied within the given radius. This keeps the -# equidistributed base grid from starving the approach to each rational, where the ideal-MHD -# Δ′ extraction is most sensitive to the local equilibrium-spline resolution. -const SING_PACK_COEF = 0.06 -const SING_PACK_RADIUS = 0.05 +# Rational-surface bracketing (Δ′ robustness). The ideal-MHD Δ′ asymptotic matching samples the +# cubic equilibrium splines' 2nd/3rd derivatives across each rational ψ_s over the matching stencil +# [ψ_s − dpsi, ψ_s + dpsi], dpsi = singfac_min/|n·q′|. A cubic 3rd derivative is piecewise constant +# and jumps at every knot, so a knot inside that stencil makes q‴ inconsistent across the crossing +# and Δ′ non-convergent. We keep the stencil inside ONE clean spline interval by bracketing each +# surface with knots at ψ_s ± w (w = BRACKET_COEF·dpsi) and clearing the zone between them — knots +# stay near the surface (kinetic-layer resolution) but never on it. BRACKET_COEF > 1 leaves margin +# for the pass-1→pass-2 shift of ψ_s. +const BRACKET_COEF = 4.0 +# Global minimum knot spacing. No grid — auto refinement, the near-surface bracket, or a fixed +# ldp/pow1 grid packing its edge at high mpsi — may place nodes closer than this. Below it, cubic +# high derivatives are dominated by field-line integration noise rather than curvature (e.g. an +# ldp grid at mpsi=2048 packs the edge to ~6e-7, deep in noise). Chosen well under the spacing of a +# converged uniform reference (~5e-4 at mpsi=2048) so it clips only pathological clustering. +const MIN_KNOT_SPACING = 1.0e-4 +# Near-rational resolution patch. The Δ′ extraction samples the equilibrium splines' 3rd +# derivatives at each rational surface, which the smooth-q measured-curvature density starves at +# mid-radius (q looks smooth there, so few knots are placed, yet Δ′ still needs fine local +# resolution). Within RATIONAL_RES_RADIUS of each rational we floor the density at 1/RATIONAL_RES_SPACING +# so equidistribution lays a locally-uniform fine patch there — the resolution Δ′ needs — which +# `bracket_mandatory_nodes` then centers the surface within. The spacing is **fixed** (τ-independent): +# Δ′ is a property of the rational surface, so its accuracy must not depend on the global accuracy +# target τ (that only sizes the pedestal/edge grid) — a τ-dependent patch would leave Δ′ far from +# converged at the default τ. RATIONAL_RES_SPACING is set where the q‴ estimate converges (the ldp +# mpsi ladder converges the q=2 Δ′ by h ≈ 5e-4). +const RATIONAL_RES_SPACING = 5.0e-4 +const RATIONAL_RES_RADIUS = 5.0e-3 +# Cap on the refined pass-2 interval count (density integral is clamped to this). +const REFINED_N_CAP = 1024 +# Pass-1 interval count for the two-pass auto grid: the layout on which the formed equilibrium's +# curvature is measured to provision the refined pass-2 grid. Pass 1 is equilibrium-only (cheap; +# the ODE integration that dominates runtime is on pass 2). Sized to measure the gradient structure +# (pedestal, edge) accurately. NOTE: because the curvature estimate is currently grid-dependent +# (it grows with the measurement grid — see the dynamic-grid issue), a larger pass 1 measures more +# curvature and enlarges pass 2; revisit once the density is made grid-invariant. +const PASS1_INTERVALS = 1024 """ wants_two_pass(config::EquilibriumConfig) -> Bool @@ -142,7 +172,7 @@ function _density_from_curvature!(rho::Vector{Float64}, d4::Vector{Float64}, f_s end """ - _knot_density(equil::PlasmaEquilibrium; tau, kin=nothing, mandatory=Float64[], sing_pack_coef=SING_PACK_COEF) -> Vector{Float64} + _knot_density(equil::PlasmaEquilibrium; tau, kin=nothing) -> Vector{Float64} Knot density ρ(ψ) (knots per unit ψ_N) at the pass-1 nodes `equil.profiles.xs`, combining by max: @@ -159,14 +189,12 @@ by max: from uniform relative q′ error on q ≈ −A·ln(1−ψ), independent of A and normally inactive; and the core for ψ ≤ 0.03, geometric-in-log(ψ) from the power-law axis form, which replaces the (noise-dominated) measurement there. A global minimum - density applies everywhere; - - local packing around each `mandatory` (rational) surface: spacing `sing_pack_coef`·τ^(1/3) - at the surface with geometric growth away from it, within `SING_PACK_RADIUS`. + density applies everywhere. -A running max over ±1 node smooths single-stencil dropouts. +Rational surfaces are handled by `bracket_mandatory_nodes` after equidistribution, not by +elevating the density here. A running max over ±1 node smooths single-stencil dropouts. """ -function _knot_density(equil::PlasmaEquilibrium; tau::Float64, kin::Union{Nothing,KineticProfileSplines}=nothing, - mandatory::Vector{Float64}=Float64[], sing_pack_coef::Float64=SING_PACK_COEF) +function _knot_density(equil::PlasmaEquilibrium; tau::Float64, kin::Union{Nothing,KineticProfileSplines}=nothing) xs = equil.profiles.xs n = length(xs) # Curvature is measured against ρ = √ψ (regular at the axis); the ρ-space density @@ -229,15 +257,6 @@ function _knot_density(equil::PlasmaEquilibrium; tau::Float64, kin::Union{Nothin rho_s[i] = max(rho_s[i], 1.0 / H_TARGET_MAX) end - # Local packing around mandatory (rational) surfaces - h_s = max(sing_pack_coef * tau^(1 / 3), H_TARGET_MIN) - for psi_m in mandatory - @inbounds for i in 1:n - d = abs(xs[i] - psi_m) - d > SING_PACK_RADIUS && continue - rho_s[i] = max(rho_s[i], 1.0 / max(dlog * d, h_s)) - end - end return rho_s end @@ -321,45 +340,129 @@ function merge_mandatory_nodes(grid::Vector{Float64}, mandatory::Vector{Float64} return merged end +""" + enforce_min_spacing(grid, hmin) -> Vector{Float64} + +Drop interior nodes so no two kept nodes are closer than `hmin`, preserving both endpoints. +Guards against noise-level packing from any generator — the auto-grid near-surface bracket, or a +fixed `ldp`/`pow1` grid whose edge spacing collapses (`~(π/2·mpsi)⁻²` for `ldp`) at high mpsi, +below which cubic high derivatives are integration noise, not curvature. A no-op when the grid is +already coarser than `hmin` everywhere (e.g. a converged uniform reference). +""" +function enforce_min_spacing(grid::Vector{Float64}, hmin::Float64) + n = length(grid) + (n < 3 || hmin <= 0) && return copy(grid) + kept = Float64[grid[1]] + for i in 2:(n-1) + grid[i] - kept[end] >= hmin && push!(kept, grid[i]) + end + # Keep the far endpoint; if the last retained interior node crowds it, drop that node instead. + length(kept) > 1 && grid[n] - kept[end] < hmin && pop!(kept) + push!(kept, grid[n]) + return kept +end + +""" + bracket_mandatory_nodes(grid, centers, min_half_widths, min_spacing; collapse_atol=1e-7) -> Vector{Float64} + +Center each rational `centers[i]` inside a single clean spline interval by replacing the knots +straddling it with a symmetric pair at `centers[i] ± w`. The half-width `w` is the larger of +`min_half_widths[i]` (the Δ′-stencil floor `BRACKET_COEF·dpsi`) and **half the local grid +spacing**, so the bracket blends into the surrounding grid rather than punching a narrow interval +into it: the region stays locally uniform with the center at its midpoint. That local uniformity +is what keeps the cubic 3rd derivative the Δ′ extraction samples both *consistent across* the +surface (no knot-on-surface jump) and *stable across* grid refinement — a narrow bracket amid +wide neighbors would be consistent but noisy. Non-bracket knots within `w` of the center, or +within `min_spacing` of either new bracket knot, are dropped; endpoints always win; bracket knots +outside the open domain are dropped, and knots closer than `collapse_atol` are collapsed to keep +the grid strictly increasing. +""" +function bracket_mandatory_nodes(grid::Vector{Float64}, centers::Vector{Float64}, min_half_widths::Vector{Float64}, min_spacing::Float64; collapse_atol::Float64=1e-7) + isempty(centers) && return copy(grid) + length(centers) == length(min_half_widths) || error("bracket_mandatory_nodes: centers and min_half_widths length mismatch") + lo, hi = grid[1], grid[end] + order = sortperm(centers) + kept = Tuple{Float64,Bool}[(g, false) for g in grid] # (value, is_bracket_knot) + last_c = -Inf + for idx in order + c, hw = centers[idx], min_half_widths[idx] + (lo < c < hi) || continue + c - last_c < collapse_atol && continue # duplicate rational (same q via several (m,n)) + k = clamp(searchsortedlast(grid, c), 1, length(grid) - 1) + w = max(hw, 0.5 * (grid[k+1] - grid[k])) # blend to half the local spacing + left, right = c - w, c + w + # Drop non-bracket, non-endpoint knots inside the zone or crowding a new bracket knot. + filter!(kept) do t + t[2] || t[1] == lo || t[1] == hi || + !(left < t[1] < right || abs(t[1] - left) < min_spacing || abs(t[1] - right) < min_spacing) + end + left > lo && push!(kept, (left, true)) + right < hi && push!(kept, (right, true)) + last_c = c + end + sort!(kept; by=first) + merged = Float64[] + for (v, _) in kept + (isempty(merged) || v - merged[end] > collapse_atol) && push!(merged, v) + end + all(diff(merged) .> 0) || error("bracket_mandatory_nodes produced a non-increasing grid") + return merged +end + """ refined_psi_grid(equil::PlasmaEquilibrium; tau, kin=nothing, mandatory=Float64[], - delta_frac=0.25, N_cap=1024, sing_pack_coef=SING_PACK_COEF) -> Vector{Float64} + singfac_min=1e-4, n_min=1, bracket_coef=BRACKET_COEF, + min_spacing=MIN_KNOT_SPACING, N_cap=1024) -> Vector{Float64} Build the refined pass-2 ψ grid from a formed pass-1 equilibrium: measured-curvature knot -density (`_knot_density`), equidistribution, and mandatory-knot insertion -(`merge_mandatory_nodes`). `tau` is the target interpolation accuracy (`psi_accuracy`); -`kin` optionally supplies kinetic profiles whose pedestal gradients attract knots; -`mandatory` lists ψ values that must appear as knots (e.g. rational surfaces); -`sing_pack_coef` scales the knot spacing at those surfaces (convergence studies only — -the default is scan-calibrated). +density (`_knot_density`), equidistribution, a global minimum-spacing floor +(`enforce_min_spacing`), and rational-surface bracketing (`bracket_mandatory_nodes`). `tau` is +the target interpolation accuracy (`psi_accuracy`); `kin` optionally supplies kinetic profiles +whose pedestal gradients attract knots; `mandatory` lists rational-surface ψ values to bracket; +`singfac_min` and `n_min` (smallest |n| in the run) set each surface's matching half-stencil +`dpsi = singfac_min/(n_min·|q′|)`, and the bracket half-width is `bracket_coef·dpsi` (floored at +`min_spacing`). Rational surfaces are bracketed, not pinned: a knot on the surface would make the +Δ′ extraction's cubic 3rd derivative jump mid-stencil (see `BRACKET_COEF`). """ function refined_psi_grid(equil::PlasmaEquilibrium; tau::Float64, kin::Union{Nothing,KineticProfileSplines}=nothing, mandatory::Vector{Float64}=Float64[], - delta_frac::Float64=0.25, - N_cap::Int=1024, - sing_pack_coef::Float64=SING_PACK_COEF) + singfac_min::Float64=1e-4, + n_min::Int=1, + bracket_coef::Float64=BRACKET_COEF, + min_spacing::Float64=MIN_KNOT_SPACING, + N_cap::Int=REFINED_N_CAP) xs = equil.profiles.xs - rho = _knot_density(equil; tau, kin, mandatory, sing_pack_coef) + rho = _knot_density(equil; tau, kin) + # Floor the density to a fixed (τ-independent) locally-uniform fine patch around each rational + # so Δ′ has the resolution to sample 3rd derivatives there at any accuracy target (see + # RATIONAL_RES_SPACING). + if !isempty(mandatory) + h_rat = max(RATIONAL_RES_SPACING, min_spacing) + @inbounds for m in mandatory, i in eachindex(xs) + abs(xs[i] - m) <= RATIONAL_RES_RADIUS && (rho[i] = max(rho[i], 1.0 / h_rat)) + end + end M_total = sum(0.5 * (rho[i] + rho[i-1]) * (xs[i] - xs[i-1]) for i in 2:length(xs)) N = clamp(ceil(Int, M_total), 32, N_cap) N == N_cap && M_total > N_cap && @warn "refined_psi_grid: knot count capped at $N_cap (density integral wants $(ceil(Int, M_total))); psi_accuracy=$tau may not be attainable" - grid = _equidistribute(xs, rho, N) - return merge_mandatory_nodes(grid, mandatory; delta_frac) + grid = enforce_min_spacing(_equidistribute(xs, rho, N), min_spacing) + isempty(mandatory) && return grid + min_half_widths = [max(bracket_coef * singfac_min / (n_min * abs(equil.profiles.q_deriv(m))), min_spacing) for m in mandatory] + return bracket_mandatory_nodes(grid, mandatory, min_half_widths, min_spacing) end """ - implied_knot_count(equil::PlasmaEquilibrium; tau, kin=nothing, mandatory=Float64[]) -> Int + implied_knot_count(equil::PlasmaEquilibrium; tau, kin=nothing) -> Int Knot count the measured-curvature density model implies for a formed equilibrium. Used as a post-refinement consistency diagnostic: if this differs substantially from the actual grid size, the pass-1 grid under- or over-sampled a feature. """ -function implied_knot_count(equil::PlasmaEquilibrium; tau::Float64, kin::Union{Nothing,KineticProfileSplines}=nothing, - mandatory::Vector{Float64}=Float64[]) +function implied_knot_count(equil::PlasmaEquilibrium; tau::Float64, kin::Union{Nothing,KineticProfileSplines}=nothing) xs = equil.profiles.xs - rho = _knot_density(equil; tau, kin, mandatory) + rho = _knot_density(equil; tau, kin) return ceil(Int, sum(0.5 * (rho[i] + rho[i-1]) * (xs[i] - xs[i-1]) for i in 2:length(xs))) end diff --git a/src/Equilibrium/InverseEquilibrium.jl b/src/Equilibrium/InverseEquilibrium.jl index 01256e8ce..725391b15 100644 --- a/src/Equilibrium/InverseEquilibrium.jl +++ b/src/Equilibrium/InverseEquilibrium.jl @@ -140,10 +140,10 @@ function equilibrium_solver(input::InverseRunInput; override_psi_nodes::Union{No mpsi = length(sq_xs) - 1 elseif grid_type in ("auto", "log_asymptotic") if mpsi == 0 && config.psi_accuracy > 0 - # Two-pass auto grid: coarse pass-1 layout; the driver refines and re-forms - # on the measured-curvature grid (GridRefinement.jl). - mpsi = 128 - @info "Auto psi grid: forming pass-1 equilibrium on coarse $(mpsi)-interval log_asymptotic grid pending curvature-based refinement" + # Two-pass auto grid: pass-1 layout sized to resolve the gradient structure so pass 2 + # is well-provisioned; the driver refines and re-forms (GridRefinement.jl, PASS1_INTERVALS). + mpsi = PASS1_INTERVALS + @info "Auto psi grid: forming pass-1 equilibrium on $(mpsi)-interval log_asymptotic grid pending curvature-based refinement" elseif mpsi == 0 mpsi = 128 end @@ -173,6 +173,13 @@ function equilibrium_solver(input::InverseRunInput; override_psi_nodes::Union{No else error("Unsupported grid_type: $grid_type") end + # Floor node spacing on the fixed grids (ldp/pow1 pack the edge below the integration-noise + # scale at high mpsi); the auto grid floors its refined grid in refined_psi_grid, and an + # override grid arrives already floored. + if override_psi_nodes === nothing && !(grid_type in ("auto", "log_asymptotic")) + sq_xs = enforce_min_spacing(sq_xs, MIN_KNOT_SPACING) + mpsi = length(sq_xs) - 1 + end sq_fs = zeros(Float64, mpsi+1, 4) sq = cubic_interp(sq_xs, Series(sq_fs); extrap=ExtendExtrap()) diff --git a/src/GeneralizedPerturbedEquilibrium.jl b/src/GeneralizedPerturbedEquilibrium.jl index 2e10a1bd5..54d474264 100755 --- a/src/GeneralizedPerturbedEquilibrium.jl +++ b/src/GeneralizedPerturbedEquilibrium.jl @@ -236,8 +236,12 @@ function main_from_inputs( # from the in-memory input — no file re-read. if Equilibrium.wants_two_pass(eq_config) mandatory = ForceFreeStates.rational_psi_nodes(equil; nlow=ctrl.nn_low, nhigh=ctrl.nn_high) + # Smallest |n| in the run sets the widest matching half-stencil dpsi = singfac_min/(n_min·|q′|), + # so the rational-surface brackets clear a zone large enough for every mode. + n_min = minimum(abs(n) for n in ctrl.nn_low:ctrl.nn_high if n != 0) psi_nodes = Equilibrium.refined_psi_grid(equil; - tau=eq_config.psi_accuracy, kin=kinetic_profiles, mandatory=mandatory) + tau=eq_config.psi_accuracy, kin=kinetic_profiles, mandatory=mandatory, + singfac_min=ctrl.singfac_min, n_min=n_min) rerun_input = if additional_input !== nothing # Analytic *Config, IMAS dd, or prebuilt RunInput — all re-formable. The IMAS # path re-runs read_imas, which must resolve the same psihigh both passes; @@ -251,7 +255,7 @@ function main_from_inputs( nothing # fall back to re-reading the input file end equil = Equilibrium.setup_equilibrium(eq_config, rerun_input; override_psi_nodes=psi_nodes) - implied = Equilibrium.implied_knot_count(equil; tau=eq_config.psi_accuracy, kin=kinetic_profiles, mandatory) + implied = Equilibrium.implied_knot_count(equil; tau=eq_config.psi_accuracy, kin=kinetic_profiles) if implied > 1.5 * (length(psi_nodes) - 1) @warn "Two-pass psi grid: refined equilibrium implies $implied knots vs $(length(psi_nodes) - 1) used — " * "pass 1 may have under-sampled a feature; consider tightening psi_accuracy" diff --git a/test/runtests_grid_refinement.jl b/test/runtests_grid_refinement.jl index 31c803119..35c350c90 100644 --- a/test/runtests_grid_refinement.jl +++ b/test/runtests_grid_refinement.jl @@ -68,6 +68,48 @@ const GridRef = GeneralizedPerturbedEquilibrium.Equilibrium @test all(diff(merged) .> 0) end + @testset "enforce_min_spacing" begin + # No-op when already coarser than the floor + g = collect(range(0.0, 1.0; length=6)) # spacing 0.2 + @test GridRef.enforce_min_spacing(g, 0.1) == g + + # Thins to respect the floor, keeping both endpoints + g = collect(range(0.0, 1.0; length=11)) # spacing 0.1 + f = GridRef.enforce_min_spacing(g, 0.25) + @test f[1] == 0.0 && f[end] == 1.0 + @test minimum(diff(f)) >= 0.25 - 1e-12 + + # Clips an edge-packed (ldp-like) grid so no interval is below the floor + ldp = [sin((i / 100) * (π / 2))^2 for i in 0:100] # packs both ends ~ (π/2·100)⁻² ≈ 2.5e-4 + @test minimum(diff(ldp)) < 1e-3 + f = GridRef.enforce_min_spacing(ldp, 1e-3) + @test f[1] == ldp[1] && f[end] == ldp[end] + @test minimum(diff(f)) >= 1e-3 - 1e-12 + end + + @testset "bracket_mandatory_nodes" begin + grid = collect(range(0.0, 1.0; length=21)) # spacing 0.05 + + # Empty centers is the identity + @test GridRef.bracket_mandatory_nodes(grid, Float64[], Float64[], 1e-4) == grid + + # A rational is centered between a bracket pair, never on a knot, in a clean interval. + b = GridRef.bracket_mandatory_nodes(grid, [0.5], [0.01], 1e-4) + @test all(diff(b) .> 0) + @test !any(g -> abs(g - 0.5) < 1e-12, b) # 0.5 is not a knot + k = searchsortedlast(b, 0.5) + @test b[k] < 0.5 < b[k+1] # interior to one interval + # bracket knots straddle the center symmetrically at ~half the local spacing (blended) + @test isapprox(0.5 - b[k], b[k+1] - 0.5; rtol=0.2) + + # Two well-separated centers each get a clean interval; grid stays strictly increasing + b = GridRef.bracket_mandatory_nodes(grid, [0.3, 0.7], [0.01, 0.01], 1e-4) + @test all(diff(b) .> 0) + for c in (0.3, 0.7) + @test !any(g -> abs(g - c) < 1e-12, b) + end + end + @testset "_validate_psi_nodes" begin grid = collect(range(0.01, 0.99; length=10)) @test GridRef._validate_psi_nodes(grid, 0.01, 0.99) === grid @@ -121,15 +163,19 @@ const GridRef = GeneralizedPerturbedEquilibrium.Equilibrium @test isapprox(q, round(2q) / 2; atol=1e-6) # m/n with n in 1:2 is a half-integer end - grid = GridRef.refined_psi_grid(equil1; tau=eq_config.psi_accuracy, mandatory=rat) + grid = GridRef.refined_psi_grid(equil1; tau=eq_config.psi_accuracy, mandatory=rat, singfac_min=1e-4, n_min=1) @test all(diff(grid) .> 0) @test grid[1] == eq_config.psilow && grid[end] == eq_config.psihigh - @test all(p -> p in grid, rat) - # Every rational surface appears as a mandatory knot in the refined grid, and the - # min-spacing snap in merge_mandatory_nodes keeps the grid strictly increasing even - # when multi-n (n=1,2) rationals cluster. - @test all(p -> any(g -> abs(g - p) < 1e-12, grid), rat) + # Rationals are bracketed, not pinned: each sits strictly interior to a clean interval — + # never on a knot — with its Δ′ matching stencil [p ± singfac_min/|n·q′|] free of knots, so + # the cubic 3rd derivative the singular extraction samples is single-valued across it. + for p in rat + dpsi = 1e-4 / abs(equil1.profiles.q_deriv(p)) + @test minimum(abs.(grid .- p)) > dpsi # no knot inside the matching stencil + k = searchsortedlast(grid, p) + @test 1 <= k < length(grid) # p interior to (grid[k], grid[k+1]) + end @test minimum(diff(grid)) > 0 # Pass 2 honors the refined grid exactly From 8f5254edf142e415d0b9a0762d8827eb7dc0167c Mon Sep 17 00:00:00 2001 From: d-burg Date: Thu, 13 Aug 2026 11:32:45 -0400 Subject: [PATCH 2/4] TEST - IMPROVEMENT - Re-pin q=2 Delta' on a converged ldp ladder (5x tighter) The q=2 pin was stale after the grid change, as its own docstring anticipated ("expect to re-pin whenever the grid generator changes"). It is now a CONVERGED value rather than a mid-spread snapshot, which closes the "replace this with a converged-Delta' pin on a fixed ldp grid" item the previous pin carried. Three independent grid families on the fixture agree to ~1%: ldp mpsi 512/1024/2048 : 8.478 / 9.142 / 9.244 uniform mpsi 1024/2048 : 9.146 / 9.143 auto tau 2e-3..2.5e-4 : 9.198 / 9.239 / 9.206 / 9.157 (flat to 0.9%) -> dpm[1,1] -> ~9.2. The old pin 6.9237 sits below even ldp-128. The platform spread also collapsed with the fix (matching stencil no longer straddles a knot): 9.239 macOS/aarch64, 9.333 linux/julia1.11, 9.356 linux/julia1.x = 1.3% end to end vs ~25% before. So rtol goes 1.5e-1 -> 3e-2, a 5x TIGHTENING. q=3/q=4 unchanged: ladder puts them at -5.716 / -16.033 (ldp 2048), inside existing tolerances. Local: 106/106 pass (runtests_parallel_integration). Co-Authored-By: Claude Opus 5 --- test/runtests_parallel_integration.jl | 36 ++++++++++++++++----------- 1 file changed, 21 insertions(+), 15 deletions(-) diff --git a/test/runtests_parallel_integration.jl b/test/runtests_parallel_integration.jl index 307615f2d..a1dd7f8f3 100644 --- a/test/runtests_parallel_integration.jl +++ b/test/runtests_parallel_integration.jl @@ -546,12 +546,18 @@ using TOML # Pinned diagonal `delta_prime_matrix` REAL parts, PEST3-convention self-response Δ' from # the STRIDE BVP with vacuum coupling, on the two-pass measured-curvature grid (rational - # surfaces pinned as mandatory knots). These pin one point at the default psi_accuracy so - # the test catches unintended changes; they are NOT converged Δ'. The extraction is - # intrinsically grid-sensitive — a psi_accuracy scan (2e-3→2.5e-4) swings dpm[1,1] by ~50% - # (6.2–9.9) and a finer ldp mpsi=512 grid gives ≈8.5 — so treat the diagonal as an - # order-of-magnitude/sign diagnostic, and expect to re-pin whenever the grid generator - # changes (as here). See the "pinning grid-sensitive Δ′ robustly" open problem in + # surfaces now *bracketed* rather than pinned on — see `bracket_mandatory_nodes`). + # + # The q=2 pin is now a CONVERGED value, which closes the "replace this with a converged-Δ′ + # pin on a fixed ldp grid" item the previous pin carried. Three independent grid families + # on this fixture agree to ~1%: + # ldp mpsi 512 / 1024 / 2048 : 8.478 / 9.142 / 9.244 + # uniform mpsi 1024 / 2048 : 9.146 / 9.143 + # auto τ 2e-3/1e-3/5e-4/2.5e-4 : 9.198 / 9.239 / 9.206 / 9.157 + # i.e. dpm[1,1] → ≈9.2. Before this grid fix the τ-scan swung ~50% (6.2–9.9) and the pin + # had to be a mid-spread snapshot; it is now flat in τ to 0.9%, so the default + # psi_accuracy lands on the converged value rather than near it by luck. + # See the "pinning grid-sensitive Δ′ robustly" open problem in # docs/src/developer_notes.md for the plateau criterion meant to replace single-point pins. # (et[1], NTV torque, and ‖resonant flux‖ stay grid-robust to <1% — the sensitivity is # local to the singular-layer matching, not the global response.) Only real parts are @@ -559,15 +565,15 @@ using TOML # FP/platform-sensitive. Near-separatrix surfaces q=5,6 keep only the finiteness/non-zero # checks above. Values use this testset's mode range (mpert=27, vs full-pipeline mpert=35). # - # q=2 additionally carries a platform spread that the other surfaces do not: the same grid - # and inputs give 6.1438 (linux/julia 1.11), 6.6867 (linux/julia 1.x), and 7.7036 - # (macOS/aarch64) — ~25% end to end. This is the surface with the known Δ' plateau problem, - # so its extraction sits on the steepest part of the grid-refinement curve and small - # platform differences in knot placement move it. The q=2 pin is therefore centered on the - # midpoint of the observed spread with rtol=1.5e-1 to span it; q=3 and q=4 are reproducible - # across the same platforms and keep rtol=1e-1. Replacing this with a converged-Δ' pin on a - # fixed ldp grid is tracked as an open issue. - @test isapprox(real(dpm[1, 1]), +6.923700e+00; rtol=1.5e-1) # q=2 (platform spread 6.14–7.70) + # The q=2 platform spread also collapsed with this fix: because the matching stencil no + # longer straddles a knot, q‴ is single-valued across it and the surface no longer sits on + # the steepest part of the refinement curve. Same grid and inputs now give 9.239 + # (macOS/aarch64), 9.333 (linux/julia 1.11) and 9.356 (linux/julia 1.x) — 1.3% end to end, + # against ~25% before (6.14 / 6.69 / 7.70). rtol=3e-2 spans that with margin and is a 5× + # TIGHTENING of the old 1.5e-1. q=3 and q=4 are unchanged: the same ladder puts them at + # −5.716 and −16.033 (ldp 2048), inside their existing tolerances, and they reproduce + # across platforms. + @test isapprox(real(dpm[1, 1]), +9.300000e+00; rtol=3e-2) # q=2, converged (ldp2048 9.244) @test isapprox(real(dpm[2, 2]), -5.344199e+00; rtol=1e-1) # q=3 @test isapprox(real(dpm[3, 3]), -1.590034e+01; rtol=1e-1) # q=4 end From b2ec04c3dc32387acadefbe585aa45897b6df168 Mon Sep 17 00:00:00 2001 From: d-burg Date: Thu, 13 Aug 2026 13:33:11 -0400 Subject: [PATCH 3/4] COILS - BUGFIX - Default the coil control surface to psihigh, not psi_N=1 (Fortran parity) MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit sample_boundary_grid (and the two pass-throughs compute_coil_forcing_modes! / convert_forcing_normalization!) defaulted to psi_N=1, which extrapolates the geometry splines past their last knot. Fortran GPEC does not do this: gpec/gpec.f:431 CALL field_bs_psi(psilim, coilmn(:,j), ...) coil/field.F:167-183 field_bs_psi = this function's analog, same splines, same rfac = SQRT(crzphi_f(1)) coil/field.F:133 Fortran's own name for the mesh: "control surface mesh" dcon/sing.f:170 psilim = psihigh, moved only inward by sas_flag/qhigh/psiedge equil/inverse.f:142 psihigh is the last knot -> Fortran never extrapolates docs/index.rst:60 external field is on the boundary "defined by the psihigh variable in equil.in" Both Fortran-comparison benchmarks already pass Fortran's own psilim explicitly (benchmark_against_fortran_run.jl:548, benchmark_coil_ForcingTerms_against_fortran.jl:236), so the psi_N=1 default was never exercised on a Fortran-validated path. Numerically, on the DIII-D-like example (psihigh=0.995) rzphi_rsquared is +0.29 at psi_N=0.995 but -22.6 at psi_N=1 (negative on 15/96 theta), so sqrt threw a DomainError. The extrapolated distance is fixed at 1-psihigh while the final spline interval shrinks with edge packing, so refining the pass-1 auto grid in this PR exposed it: runtests_coils passes 8/8 at PASS1_INTERVALS=128 and errors at 512/1024. Fixing the root cause keeps PASS1_INTERVALS=1024 and leaves the Delta' convergence untouched (dpm[1,1] = 9.246/9.291/9.239 at 256/512/1024). Impact: coupling surface moves inward by (1-psihigh) — 0.05% of flux at the default psihigh=0.9995. Callers passing psi explicitly are unaffected. Known gap: strictly faithful is psilim (integration limit), not psihigh (spline limit); they differ when dmlim/qhigh/psiedge fire. psihigh is all a bare PlasmaEquilibrium knows — callers with an FFS solve should pass intr.psilim. Full test suite passes locally. Co-Authored-By: Claude Opus 5 --- src/ForcingTerms/CoilFourier.jl | 49 +++++++++++++++++++++++++++++---- 1 file changed, 43 insertions(+), 6 deletions(-) diff --git a/src/ForcingTerms/CoilFourier.jl b/src/ForcingTerms/CoilFourier.jl index 19366a9d1..97382a5e5 100644 --- a/src/ForcingTerms/CoilFourier.jl +++ b/src/ForcingTerms/CoilFourier.jl @@ -47,13 +47,50 @@ struct BoundaryGrid end """ - sample_boundary_grid(equil, mtheta, nzeta; psi=1.0) -> BoundaryGrid + sample_boundary_grid(equil, mtheta, nzeta; psi=equil.rzphi_xs[end]) -> BoundaryGrid Evaluate plasma geometry at a uniform (mtheta × nzeta) grid on the flux surface `psi`. Uses `equil.rzphi_rsquared` and `equil.rzphi_offset` splines at the given `psi`. -Defaults to `psi=1.0` (true plasma boundary). Use a smaller value (e.g. the -Fortran `psilim`) to evaluate on an interior truncation surface. + +Defaults to the **outermost computed surface** `equil.rzphi_xs[end]` (i.e. `psihigh`), NOT +ψ_N = 1, matching Fortran GPEC. + +**Fortran reference** (PrincetonUniversity/GPEC v1.5.7): the coil field is evaluated at +`psilim`, never at ψ_N = 1 — `CALL field_bs_psi(psilim, coilmn(:,j), ...)` (`gpec/gpec.f:431`). +`field_bs_psi` is the direct analog of this function, using the same two splines and the same +square root: `rfac = SQRT(crzphi_f(1))` (`coil/field.F:170`, where `crzphi_f(1)` is +`rzphi_rsquared` and `crzphi_f(2)` is `rzphi_offset`); `coil/field.F:133` calls that mesh the +"control surface mesh". `psilim = psihigh` (`dcon/sing.f:170`) and is only ever moved *inward* +by `sas_flag`/`qhigh`/`psiedge` truncation. Since `psihigh` is the last knot of the radial +grid these splines are built on (`equil/inverse.f:142`), Fortran evaluates exactly ON the last +knot and never extrapolates. The docs state it directly: the external field is specified "on +the surface of the GPEC plasma boundary defined by the psihigh variable in equil.in" +(`docs/index.rst:60`). The previous ψ_N = 1 default was a port divergence. + +Both Fortran-comparison benchmarks already pass the correct surface explicitly +(`benchmark_against_fortran_run.jl:548` and `benchmark_coil_ForcingTerms_against_fortran.jl:236`, +both using Fortran's own `psilim` attribute), so the ψ_N = 1 default was never exercised on any +Fortran-validated path — only by callers that omitted the keyword. + +Why it mattered numerically: the geometry splines are defined only out to `psihigh`, so ψ_N = 1 +extrapolates a cubic past its last knot with no guarantee of remaining physical. On the +DIII-D-like example (psihigh = 0.995) `rzphi_rsquared` is +0.29 at ψ_N = 0.995 but −22.6 at +ψ_N = 1, negative over 15 of 96 θ points, so the `sqrt` below throws a DomainError. The +extrapolated distance is a fixed 1 − psihigh while the final spline interval shrinks with edge +packing, so the error grows with resolution — which is why this surfaced only once the pass-1 +auto grid was refined. + +NOTE: strictly faithful is `psilim`, the *integration* limit, not `psihigh`, the *equilibrium +spline* limit. They are equal unless `dmlim`/`qhigh`/`psiedge` truncation fires, in which case +`psilim < psihigh`. PPPL shipped a fix for exactly this confusion (`docs/releases.rst:281`: +"Fixes inappropriate uses of psihigh, which may not be the end of integration psilim if +sas_flag, qhigh, or peak_flag are used"). `psihigh` is used as the default here because it is +all a bare `PlasmaEquilibrium` knows; callers holding a ForceFreeStates solve should pass +`psi = intr.psilim` explicitly, as the benchmarks do. + +The coupling surface therefore moves inward by (1 − psihigh): 0.05 % of flux at the default +psihigh = 0.9995, 0.5 % on the example, which lowers psihigh to 0.995 to capture q=6. Poloidal derivatives dR/dθ_norm and dZ/dθ_norm (unit-norm θ_norm ∈ [0,1]) are computed via periodic cubic splines on the resulting R(θ_norm), Z(θ_norm) data. @@ -64,7 +101,7 @@ For DIII-D (Bt < 0, Ip > 0 → helicity = -1): phi increases with j (standard di For positive-helicity machines (Bt > 0, Ip > 0 → helicity = +1): phi decreases with j. """ function sample_boundary_grid(equil::Equilibrium.PlasmaEquilibrium, mtheta::Int, nzeta::Int; - psi::Float64=1.0) + psi::Float64=equil.rzphi_xs[end]) # Build uniform theta grid (same convention as equil.rzphi_ys, but potentially finer) theta_grid = range(0; length=mtheta, step=1.0/mtheta) @@ -226,7 +263,7 @@ function compute_coil_forcing_modes!( n::Int, m_low::Int, m_high::Int; - psi::Float64=1.0, + psi::Float64=equil.rzphi_xs[end], # outermost computed surface (psihigh), not ψ_N=1 — see sample_boundary_grid verbose::Bool=false ) nzeta = cfg.nzeta_coil > 0 ? cfg.nzeta_coil : NZETA_POINTS_PER_PERIOD * max(1, abs(n)) @@ -307,7 +344,7 @@ function convert_forcing_normalization!( n::Int, m_low::Int, m_high::Int; - psi::Float64=1.0, + psi::Float64=equil.rzphi_xs[end], # outermost computed surface (psihigh), not ψ_N=1 — see sample_boundary_grid mtheta::Int=256, nzeta::Int=64 ) From 8dab0ed9339a91496dc19a9f709877e9dfb9b121 Mon Sep 17 00:00:00 2001 From: d-burg Date: Thu, 13 Aug 2026 14:24:09 -0400 Subject: [PATCH 4/4] COILS - BUGFIX - Use psilim on the external-forcing path too, not just the coil path MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit The coil branch of PerturbedEquilibrium already passed psi=ffs_intr.psilim (PerturbedEquilibrium.jl:98), but the other branch — external forcing data read from file — called convert_forcing_normalization! with no psi, so it fell back to the default surface. That path therefore normalized on the wrong surface silently, with no error, whenever dmlim/qhigh/psiedge moved psilim inward of psihigh (e.g. 0.998902 vs 0.9995 on a DIII-D case: a 6e-4 gap, same order as the psi_N=1 extrapolation this PR already fixed). Both production paths now use psilim, matching Fortran (gpec/gpec.f:431 field_bs_psi(psilim,...)). psilim is the ODE's last grid point to within eps=1e-10 (the final chunk ends at psilim*(1-eps), EulerLagrange.jl:504) and is snapped to psi_store[end] exactly in the dW-peak truncation path (Riccati.jl:1375,1870, EulerLagrange.jl:230). The psihigh default on the three signatures is retained and re-documented as the bare-PlasmaEquilibrium fallback — reachable only by direct callers with no ForceFreeStates solve, for whom psilim does not exist yet — not as the intended control surface. Full test suite passes locally. Co-Authored-By: Claude Opus 5 --- src/ForcingTerms/CoilFourier.jl | 18 +++++++++++------- .../PerturbedEquilibrium.jl | 6 +++++- 2 files changed, 16 insertions(+), 8 deletions(-) diff --git a/src/ForcingTerms/CoilFourier.jl b/src/ForcingTerms/CoilFourier.jl index 97382a5e5..62996bf59 100644 --- a/src/ForcingTerms/CoilFourier.jl +++ b/src/ForcingTerms/CoilFourier.jl @@ -81,13 +81,17 @@ extrapolated distance is a fixed 1 − psihigh while the final spline interval s packing, so the error grows with resolution — which is why this surfaced only once the pass-1 auto grid was refined. -NOTE: strictly faithful is `psilim`, the *integration* limit, not `psihigh`, the *equilibrium -spline* limit. They are equal unless `dmlim`/`qhigh`/`psiedge` truncation fires, in which case -`psilim < psihigh`. PPPL shipped a fix for exactly this confusion (`docs/releases.rst:281`: -"Fixes inappropriate uses of psihigh, which may not be the end of integration psilim if -sas_flag, qhigh, or peak_flag are used"). `psihigh` is used as the default here because it is -all a bare `PlasmaEquilibrium` knows; callers holding a ForceFreeStates solve should pass -`psi = intr.psilim` explicitly, as the benchmarks do. +NOTE ON THE DEFAULT: the physically correct control surface is `psilim`, the *integration* +limit, not `psihigh`, the *equilibrium spline* limit. They are equal unless +`dmlim`/`qhigh`/`psiedge` truncation fires, in which case `psilim < psihigh`. PPPL shipped a fix +for exactly this confusion (`docs/releases.rst:281`: "Fixes inappropriate uses of psihigh, which +may not be the end of integration psilim if sas_flag, qhigh, or peak_flag are used"). + +**Both production paths pass `psi = ffs_intr.psilim` explicitly** (PerturbedEquilibrium.jl, coil +and forcing-file branches), as do both Fortran-comparison benchmarks. This default of `psihigh` +therefore applies only to direct callers holding a bare `PlasmaEquilibrium` with no +ForceFreeStates solve — for which it is the outermost surface that exists, and always a strict +improvement on extrapolating to ψ_N = 1. Pass `psi` explicitly whenever `psilim` is known. The coupling surface therefore moves inward by (1 − psihigh): 0.05 % of flux at the default psihigh = 0.9995, 0.5 % on the example, which lowers psihigh to 0.995 to capture q=6. diff --git a/src/PerturbedEquilibrium/PerturbedEquilibrium.jl b/src/PerturbedEquilibrium/PerturbedEquilibrium.jl index f94eabbc8..659f979ba 100644 --- a/src/PerturbedEquilibrium/PerturbedEquilibrium.jl +++ b/src/PerturbedEquilibrium/PerturbedEquilibrium.jl @@ -107,8 +107,12 @@ function compute_perturbed_equilibrium( modes_n = filter(m -> m.n == n, intr.forcing_modes) isempty(modes_n) && continue m_vals = [m.m for m in modes_n] + # Same control surface as the coil branch above: psilim, the integration + # limit (Fortran: gpec/gpec.f:431 `field_bs_psi(psilim, ...)`). Without it + # the normalization was taken on the equilibrium-spline limit, which differs + # whenever dmlim/qhigh/psiedge truncation moves psilim inward. convert_forcing_normalization!(modes_n, norm_tag, equil, n, - minimum(m_vals), maximum(m_vals)) + minimum(m_vals), maximum(m_vals); psi=ffs_intr.psilim) end end end