EQUIL + GPEC - Two-pass measured-curvature auto psi grid with rational-surface knots

benchmarks/benchmark_q_vs_iota_edge.jl

@@ -0,0 +1,108 @@
+# Quantifies whether splining iota = 1/q instead of q would reduce the knot count needed
+# for a given edge accuracy (GitHub discussion around the auto-grid redesign; PR #179
+# lineage owns any production iota work). For a grid ending at psihigh < 1, interpolation
+# error transforms as err_q ≈ err_iota · q², cancelling the flattening of the profile, so
+# no leading-order gain is expected — this script documents that with numbers on the
+# DIII-D-like example equilibrium.
+#
+# Usage: julia --project=. benchmarks/benchmark_q_vs_iota_edge.jl
+# Outputs (not committed): benchmarks/q_vs_iota_edge.csv, benchmarks/q_vs_iota_edge.png
+
+using Pkg;
+Pkg.activate(joinpath(@__DIR__, ".."))
+using GeneralizedPerturbedEquilibrium
+using FastInterpolations
+using Printf
+using Plots
+
+const GPE = GeneralizedPerturbedEquilibrium
+const EXAMPLE_DIR = joinpath(@__DIR__, "..", "examples", "DIIID-like_ideal_example")
+
+# 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
+    equil = GPE.Equilibrium.setup_equilibrium(eq_config, additional_input)
+    return equil, eq_config
+end
+
+# Rational surface ψ_s(q = m/n) and q'(ψ_s) recovered from a fitted spline by bisection
+function rational_locations(q_itp, dq_itp, psis, q_targets)
+    out = Dict{Float64,Tuple{Float64,Float64}}()
+    for qt in q_targets
+        lo, hi = psis[1], psis[end]
+        (q_itp(lo) - qt) * (q_itp(hi) - qt) > 0 && continue
+        for _ in 1:200
+            mid = 0.5 * (lo + hi)
+            (q_itp(lo) - qt) * (q_itp(mid) - qt) <= 0 ? (hi = mid) : (lo = mid)
+        end
+        ps = 0.5 * (lo + hi)
+        out[qt] = (ps, dq_itp(ps))
+    end
+    return out
+end
+
+function main()
+    equil, eq_config = reference_q()
+    xs_ref = equil.profiles.xs
+    q_ref_itp = equil.profiles.q_spline
+    dq_ref_itp = equil.profiles.q_deriv
+    psilow, psihigh = xs_ref[1], xs_ref[end]
+
+    q_lo = ceil(q_ref_itp(psilow))
+    q_hi = floor(q_ref_itp(psihigh))
+    q_targets = collect(q_lo:q_hi)  # n=1 rationals
+    ref_rats = rational_locations(q_ref_itp, dq_ref_itp, xs_ref, q_targets)
+
+    psi_eval = collect(range(psilow, psihigh; length=4001))
+    edge_band = psi_eval .> 0.9
+
+    rows = String[]
+    push!(rows, "N,repr,max_rel_q_err,edge_rel_q_err,edge_rel_qprime_err,max_psi_s_err,max_rel_q1_err")
+    results = Dict{String,Vector{NTuple{2,Float64}}}("q" => [], "iota" => [])
+
+    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
+        knots = GPE.Equilibrium._build_psi_grid(cfg, psilow, psihigh)
+        q_nodes = [q_ref_itp(p) for p in knots]
+
+        for repr in ("q", "iota")
+            itp = repr == "q" ? cubic_interp(knots, q_nodes; extrap=ExtendExtrap()) :
+                  cubic_interp(knots, 1.0 ./ q_nodes; extrap=ExtendExtrap())
+            ditp = deriv1(itp)
+            qf = repr == "q" ? (p -> itp(p)) : (p -> 1.0 / itp(p))
+            dqf = repr == "q" ? (p -> ditp(p)) : (p -> -ditp(p) / itp(p)^2)  # q' = -ι'/ι²
+
+            q_err = [abs(qf(p) - q_ref_itp(p)) / abs(q_ref_itp(p)) for p in psi_eval]
+            dq_err = [abs(dqf(p) - dq_ref_itp(p)) / max(abs(dq_ref_itp(p)), 1e-10) for p in psi_eval]
+            rats = rational_locations(qf, dqf, knots, q_targets)
+            psi_s_err = maximum(abs(rats[qt][1] - ref_rats[qt][1]) for qt in keys(ref_rats); init=0.0)
+            q1_err = maximum(abs(rats[qt][2] - ref_rats[qt][2]) / abs(ref_rats[qt][2]) for qt in keys(ref_rats); init=0.0)
+
+            push!(rows, @sprintf("%d,%s,%.3e,%.3e,%.3e,%.3e,%.3e",
+                N, repr, maximum(q_err), maximum(q_err[edge_band]), maximum(dq_err[edge_band]), psi_s_err, q1_err))
+            push!(results[repr], (Float64(N), maximum(q_err[edge_band])))
+        end
+    end
+
+    csv_path = joinpath(@__DIR__, "q_vs_iota_edge.csv")
+    open(csv_path, "w") do io
+        foreach(r -> println(io, r), rows)
+    end
+    println("Wrote ", abspath(csv_path))
+    foreach(println, rows)
+
+    p = plot(; xscale=:log10, yscale=:log10, xlabel="knots N", ylabel="max relative q error, ψ > 0.9",
+        title="q-spline vs ι-spline edge accuracy", legend=:topright, left_margin=12Plots.mm, bottom_margin=4Plots.mm)
+    for (repr, pts) in results
+        plot!(p, first.(pts), last.(pts); marker=:circle, lw=2, label="$repr-spline")
+    end
+    png_path = joinpath(@__DIR__, "q_vs_iota_edge.png")
+    savefig(p, png_path)
+    println("Wrote ", abspath(png_path))
+end
+
+main()

benchmarks/plot_grid_knot_placement.jl

@@ -0,0 +1,67 @@
+# Visualizes where the two-pass auto grid places its radial knots compared to the fixed
+# ldp (sin²) grids, on the DIII-D-like example. Two panels: cumulative knot fraction
+# index/mpsi vs ψ_N (packing profile), and local knot spacing Δψ vs ψ_N (log scale),
+# with the n=1 rational surfaces marked. Documents the auto-grid PR.
+#
+# Usage: julia --project=. benchmarks/plot_grid_knot_placement.jl
+# Output (not committed): benchmarks/grid_knot_placement.png
+
+using Pkg;
+Pkg.activate(joinpath(@__DIR__, ".."))
+using GeneralizedPerturbedEquilibrium
+using Plots
+
+const GPE = GeneralizedPerturbedEquilibrium
+const EXAMPLE_DIR = joinpath(@__DIR__, "..", "examples", "DIIID-like_ideal_example")
+
+function build_grids()
+    _, eq_config, additional_input = GPE.build_inputs_from_toml(EXAMPLE_DIR)
+    tau = eq_config.psi_accuracy
+
+    # Pass 1 (coarse layout) and the two-pass refined grid, as the driver builds them
+    equil1 = GPE.Equilibrium.setup_equilibrium(eq_config, additional_input)
+    pass1 = copy(equil1.profiles.xs)
+    rationals = GPE.ForceFreeStates.rational_psi_nodes(equil1; nlow=1, nhigh=1)
+    auto = GPE.Equilibrium.refined_psi_grid(equil1; tau, mandatory=rationals)
+
+    ldp(mpsi) = [eq_config.psilow + (eq_config.psihigh - eq_config.psilow) * sin((i / mpsi) * (π / 2))^2 for i in 0:mpsi]
+
+    grids = [
+        ("auto two-pass, τ=$(tau), mpsi=$(length(auto) - 1)", auto),
+        ("pass-1 layout, mpsi=$(length(pass1) - 1)", pass1),
+        ("ldp, mpsi=256", ldp(256)),
+        ("ldp, mpsi=512", ldp(512))
+    ]
+    return grids, rationals
+end
+
+function main()
+    grids, rationals = build_grids()
+    # Okabe-Ito colorblind-safe hues, fixed order; pass-1 recessive gray
+    colors = ["#0072B2", "#999999", "#E69F00", "#009E73"]
+    styles = [:solid, :dash, :solid, :solid]
+
+    p1 = plot(; xlabel="ψ_N", ylabel="knot index / mpsi", title="Radial knot packing — DIIID-like example (n=1)",
+        legend=:bottomright, left_margin=12Plots.mm, bottom_margin=4Plots.mm, size=(900, 420))
+    p2 = plot(; xlabel="ψ_N", ylabel="local knot spacing Δψ_N", yscale=:log10,
+        legend=:bottomleft, left_margin=12Plots.mm, bottom_margin=6Plots.mm, size=(900, 420))
+
+    for (i, (label, g)) in enumerate(grids)
+        n = length(g) - 1
+        plot!(p1, g, (0:n) ./ n; lw=2, color=colors[i], ls=styles[i], label=label)
+        mids = 0.5 .* (g[1:(end-1)] .+ g[2:end])
+        plot!(p2, mids, diff(g); lw=2, color=colors[i], ls=styles[i], label=label)
+    end
+    for (j, r) in enumerate(rationals)
+        for p in (p1, p2)
+            vline!(p, [r]; color="#CC79A7", ls=:dot, lw=1.5, label=(j == 1 && p === p1) ? "rational surfaces q=m/1" : "")
+        end
+    end
+
+    fig = plot(p1, p2; layout=(2, 1), size=(900, 840))
+    out = joinpath(@__DIR__, "grid_knot_placement.png")
+    savefig(fig, out)
+    println("Wrote ", abspath(out))
+end
+
+main()

docs/src/equilibrium.md

@@ -57,21 +57,66 @@ tuning is needed.
 
 ## Radial grid packing
 
-The default `grid_type = "log_asymptotic"` uses a three-region grid that respects
-the asymptotic behavior of q near both the magnetic axis and the separatrix:
-
-- **Core** (ψ < 0.15): geometric spacing in log(ψ) — handles the axis where profiles
-  behave as ψⁿ
-- **Middle** (0.15 ≤ ψ ≤ 0.95): uniform spacing — preserves resolution through the
-  pedestal region
-- **Edge** (ψ > 0.95): geometric spacing in log(1−ψ) — tracks the logarithmic
-  divergence q ~ −A·ln(1−ψ) near a diverted separatrix
-
-With `mpsi = 0` (the default), the number of radial knots is chosen automatically
-from the `psi_accuracy` parameter (target absolute error in q). Two probe field-line
-integrations near psihigh estimate the local log-slope A, and the knot count is set so
-the cubic spline error stays below `psi_accuracy` throughout the domain. The legacy
-`grid_type = "ldp"` (sin²-spaced) and explicit `mpsi` are still supported.
+With `grid_type = "auto"` and `mpsi = 0` (the defaults; `"log_asymptotic"` is a legacy alias), the radial grid is
+built by a **two-pass measured-curvature refinement** driven by the single accuracy knob
+`psi_accuracy` (τ):
+
+1. **Pass 1** forms the equilibrium on a coarse three-region layout (geometric in
+   log(ψ) at the core, uniform in the middle, geometric in log(1−ψ) at the edge).
+2. A knot density is derived from **one sizing rule** — the cubic-spline derivative
+   error model err(f′) ≈ h³|f''''|/24 ≤ τ·|f| — applied to three sources of
+   curvature |f''''|:
+   - **Measured** (mid-radius and edge): divided differences on the pass-1 nodes of the
+     1D profiles (F, P, dV/dψ, q), the 2D geometry channels (r², η-offset, ν, Jacobian)
+     along sampled θ-lines (the bicubic ψ-axis shares the same knots), and the kinetic
+     profiles (n, T, ω_E) when loaded, so steep pedestal gradients attract knots. Nodal
+     values come from independent field-line integrations, so the estimate is immune to
+     inter-knot spline ringing. Curvature is measured against ρ = √ψ, in which the
+     equilibrium is regular at the magnetic axis — in ψ itself the geometry channels
+     behave as ψ^(k/2) (R−R₀ ~ √ψ), so their ψ-curvature diverges under refinement while
+     their ρ-curvature converges; the ρ-spacing maps back through dψ/dρ = 2√ψ, which by
+     itself packs √ψ-tight toward the axis.
+   - **Separatrix model** (edge floor, ψ ≥ 0.9): the same rule on q ≈ −A·ln(1−ψ) gives
+     geometric-in-log(1−ψ) packing with uniform relative q′ error, independent of A.
+     Normally inactive — the pass-1 layout is already log-packed at the edge, so the
+     measured density dominates — it only guards a pass-1 that under-sampled the
+     divergence.
+   - **Axis model** (core, ψ ≤ 0.03): the same rule on the power-law axis form gives
+     geometric-in-log(ψ) packing (constant-ratio spacing ∝ ψ down to `psilow`). Here the
+     model *replaces* measurement: nodal data from the smallest flux surfaces is
+     dominated by integration and axis-extrapolation noise, which grid refinement
+     amplifies rather than resolves.
+   - **Rational surfaces**: local packing around every q = m/n in the requested n range
+     (pinned as mandatory knots by the main driver), at any ψ including inside the
+     modeled core region. This local resolution of the ψ-splined Euler-Lagrange
+     coefficient matrices is what converges the Δ′ boundary-value problem.
+3. The density integral is equidistributed into the knot vector, mandatory
+   rational-surface knots are inserted with a minimum-spacing snap guard, and
+   **pass 2** re-forms the equilibrium on the refined grid from the in-memory input
+   (no file re-read).
+
+Every region's knot count scales as τ^(-1/3), so tightening `psi_accuracy` refines
+the core, pedestal, and edge proportionally. The legacy `grid_type = "ldp"`
+(sin²-spaced), `"pow1"`, `"uniform"`, and explicit `mpsi > 0` (single-pass, fixed
+layout) are still supported. Library users calling `setup_equilibrium` directly with
+`mpsi = 0` receive the coarse pass-1 grid; use `refined_psi_grid` and the
+`override_psi_nodes` keyword to apply the refinement manually.
+
+The packing on the DIII-D-like example (n=1) compared to fixed `ldp` grids — note the
+coarse spacing across the smooth mid-radius, the spacing dips at each rational surface,
+and the core/pedestal/edge packing (`benchmarks/plot_grid_knot_placement.jl` regenerates
+this figure):
+
+![Radial knot packing: auto two-pass vs ldp](assets/grid_knot_placement.png)
+
+Decomposing the density by source on the same example shows the pedestal band
+(ψ_N ≈ 0.85–0.98) is driven by *measured* curvature, not the edge floor: the pressure,
+q, and dV/dψ profiles contribute comparably, and the rzphi geometry channels — the
+Grad-Shafranov response to the same pedestal p′ (Shafranov-shift / angle-offset
+steepening) — contribute the most at every node in the band. Because the density is
+measured from the formed solution, the packing follows the pedestal wherever it sits:
+
+![Knot-density decomposition by source](assets/density_decomposition.png)
 
 ## API Reference
 

examples/DIIID-like_gal_resistive_example/gpec.toml

@@ -8,8 +8,8 @@ power_r = 0                    # Major radius power exponent for Jacobian
 grid_type = "ldp"              # Radial grid packing type (ldp = linear-derivative packing toward rationals)
 psilow = 1e-4                  # Lower limit of normalized flux coordinate
 psihigh = 0.993                # Upper limit of normalized flux coordinate (0.993 stays clear of the separatrix; truncating at ≳0.998 is numerically unreasonable here)
-mpsi = 128                     # Number of radial grid points (0 = auto-compute from psi_accuracy)
-psi_accuracy = 0.001           # Target absolute error in q for auto-mpsi
+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)
 etol = 1e-10                   # Error tolerance for equilibrium solver

examples/DIIID-like_gal_resistive_pe_example/gpec.toml

@@ -8,8 +8,8 @@ power_r = 0                    # Major radius power exponent for Jacobian
 grid_type = "ldp"              # Radial grid packing type (ldp = linear-derivative packing toward rationals)
 psilow = 1e-4                  # Lower limit of normalized flux coordinate
 psihigh = 0.993                # Upper limit of normalized flux coordinate (0.993 stays clear of the separatrix; truncating at ≳0.998 is numerically unreasonable here)
-mpsi = 128                     # Number of radial grid points (0 = auto-compute from psi_accuracy)
-psi_accuracy = 0.001           # Target absolute error in q for auto-mpsi
+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)
 etol = 1e-10                   # Error tolerance for equilibrium solver

examples/DIIID-like_ideal_example/gpec.toml

@@ -7,11 +7,11 @@
 eq_filename = "TkMkr_D3Dlike_Hmode.geqdsk" # Path to equilibrium file
 eq_type = "efit"               # Type of the input 2D equilibrium file
 jac_type = "hamada"            # Coordinate system (hamada, pest, boozer, equal_arc, park, other)
-grid_type = "log_asymptotic"   # Radial grid packing type
+grid_type = "auto"   # Radial grid packing type ("auto" = two-pass measured-curvature grid)
 psilow = 1e-4                  # Lower limit of normalized poloidal flux
 psihigh = 0.9995               # Upper limit of normalized poloidal flux
-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
+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)
 etol = 1e-10                   # Error tolerance for equilibrium solver

examples/DIIID-like_ideal_example_IMAS/gpec.toml

@@ -13,7 +13,7 @@ eq_filename  = "from_dd"       # Path to equilibrium file (placeholder: equilibr
 jac_type     = "boozer"        # Coordinate system (hamada, pest, boozer, equal_arc, park, other)
 psilow       = 0.01            # Lower limit of normalized poloidal flux
 psihigh      = 0.9995          # Upper limit of normalized poloidal flux
-mpsi         = 128             # Number of radial grid points (0 = auto-compute from psi_accuracy)
+mpsi         = 128             # Number of radial grid intervals (0 = two-pass auto grid from psi_accuracy)
 mtheta       = 256             # Number of poloidal grid points
 imas_cocos   = 11              # Expected COCOS convention of the input IMAS dd (11 = IMAS standard)
 

examples/LAR_beta_scan/gpec.toml

@@ -13,7 +13,7 @@ jac_type  = "hamada"           # Coordinate system (hamada, pest, boozer, equal_
 grid_type = "ldp"              # Radial grid packing type
 psilow    = 0.01               # Lower limit of normalized poloidal flux
 psihigh   = 0.995              # Upper limit of normalized poloidal flux
-mpsi      = 128                # Number of radial grid points (0 = auto-compute from psi_accuracy)
+mpsi      = 128                # Number of radial grid intervals (0 = two-pass auto grid from psi_accuracy)
 mtheta    = 512                # Number of poloidal grid points
 
 # TJ-analytic equilibrium parameters (cf. R. Fitzpatrick, TJ code, https://github.com/rfitzp/TJ).

examples/LAR_epsilon_scan/gpec.toml

@@ -13,7 +13,7 @@ jac_type  = "hamada"           # Coordinate system (hamada, pest, boozer, equal_
 grid_type = "ldp"              # Radial grid packing type
 psilow    = 0.01               # Lower limit of normalized poloidal flux
 psihigh   = 0.995              # Upper limit of normalized poloidal flux
-mpsi      = 128                # Number of radial grid points (0 = auto-compute from psi_accuracy)
+mpsi      = 128                # Number of radial grid intervals (0 = two-pass auto grid from psi_accuracy)
 mtheta    = 512                # Number of poloidal grid points
 
 # TJ-analytic equilibrium parameters (cf. R. Fitzpatrick, TJ code, https://github.com/rfitzp/TJ).

examples/Solovev_ideal_example/gpec.toml

@@ -9,7 +9,7 @@ jac_type = "pest"                       # Coordinate system (hamada, pest, booze
 grid_type = "ldp"                       # Radial grid packing type
 psilow = 1e-4                           # Lower limit of normalized poloidal flux
 psihigh = 0.9995                        # Upper limit of normalized poloidal flux
-mpsi = 128                              # Number of radial grid points (0 = auto-compute from psi_accuracy)
+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)
 etol = 1e-7                             # Error tolerance for equilibrium solver

examples/Solovev_ideal_example_3D/gpec.toml

@@ -8,7 +8,7 @@ jac_type = "pest"                       # Coordinate system (hamada, pest, booze
 grid_type = "ldp"                       # Radial grid packing type
 psilow = 1e-4                           # Lower limit of normalized poloidal flux
 psihigh = 0.9995                        # Upper limit of normalized poloidal flux
-mpsi = 128                              # Number of radial grid points (0 = auto-compute from psi_accuracy)
+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)
 etol = 1e-7                             # Error tolerance for equilibrium solver