diff --git a/src/Equilibrium/DirectEquilibrium.jl b/src/Equilibrium/DirectEquilibrium.jl
index 02f18410..24b78f21 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 03935640..cf08d8a4 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 9478d96b..e58eb81d 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 01256e8c..725391b1 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 2e10a1bd..54d47426 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 31c80311..35c350c9 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