From 7b1386b45b9472748db6eacfa5888572d538f98a Mon Sep 17 00:00:00 2001 From: logan-nc Date: Fri, 4 Sep 2026 13:52:13 -0400 Subject: [PATCH] EQUIL - BUGFIX - Sample knot midpoints so the round-trip check can see spline ringing The efit_by_inversion round-trip residual (psi,theta) -> (R,Z) -> psi was sampled only at psi knots -- where the rzphi cubic splines interpolate exactly. Inter-knot ringing is invisible to that measurement by construction, so the check reported a clean edge while the reconstruction was unusable. Measured on the DIII-D-like deck at fixed mpsi=512, pushing psihigh toward 1: psihigh on-knot midpoint ratio et[1] 0.995 4.10e-06 5.02e-06 1.2 +0.80 0.999 3.57e-06 4.19e-06 1.2 +0.81 0.9999 1.14e-05 1.29e-03 113.1 +4.70 0.99999 2.65e-05 2.18e-02 823.2 -47707 The on-knot residual is flat across the whole range while the stability energies collapse by six orders of magnitude. The midpoint residual tracks the degradation, and their ratio is what actually separates the two regimes: dimensionless, so it needs no per-grid calibration, and good/bad sit two orders either side of 10. An absolute tolerance alone does not substitute. At psihigh=0.9999 the midpoint residual is 1.29e-3, under the existing 2e-3, so only the ratio catches that rung. Conversely the ratio alone would false-fire on a clean grid where both residuals are at the rounding floor, so it is consulted only once the midpoint residual clears ROUNDTRIP_TOL/20. Sampling widened from two fixed knots to every knot and midpoint in the outer ROUNDTRIP_EDGE_FRAC of the grid -- the region the check is named for and where the failure lives. The wider window matters: with only the two previously sampled knots and their neighbouring midpoints, the ratio at psihigh=0.99999 reads 1.0 and misses the worst rung entirely, because the ringing sits in an interval those samples do not touch. Verification: runtests_equil.jl passes 286/286, including a new false-positive control that a cleanly traced deck still reports at Info. Co-Authored-By: Claude Opus 5 (1M context) Claude-Session: https://claude.ai/code/session_01PSrf6JCViFfVzqzkQ66o6b --- .../DirectEquilibriumByInversion.jl | 63 ++++++++++++++----- test/runtests.jl | 1 + test/runtests_equil.jl | 33 ++++++++++ 3 files changed, 82 insertions(+), 15 deletions(-) diff --git a/src/Equilibrium/DirectEquilibriumByInversion.jl b/src/Equilibrium/DirectEquilibriumByInversion.jl index 136d8e6ac..95f2e7ed7 100644 --- a/src/Equilibrium/DirectEquilibriumByInversion.jl +++ b/src/Equilibrium/DirectEquilibriumByInversion.jl @@ -10,6 +10,20 @@ Select via `eq_type = "efit_by_inversion"` in `gpec.toml`. import Contour as Ctr +# Round-trip (ψ,θ)→(R,Z)→ψ residual above which the traced edge geometry is not trusted. +const ROUNDTRIP_TOL = 2e-3 +# Fraction of the radial grid treated as "the edge" for the round-trip check. +const ROUNDTRIP_EDGE_FRAC = 0.9 +# How much larger the midpoint residual may be than the on-knot one before the rzphi splines are +# judged to be ringing between knots. Grid-invariant, being a ratio of two residuals measured the +# same way. Calibrated on a DIII-D-like psihigh ladder, where it reads 1.2 while the ForceFreeStates +# energies are sound and 10^2-10^3 once they are not; a separation that wide leaves the exact +# threshold uncritical. +const ROUNDTRIP_RATIO_TOL = 10.0 +# Midpoint residual below which the ratio is not consulted, since both residuals are then at the +# rounding floor and their ratio measures noise rather than ringing. +const ROUNDTRIP_RINGING_FLOOR = ROUNDTRIP_TOL / 20 + """ _is_closed_curve(curve) @@ -692,27 +706,46 @@ function equilibrium_solver_by_inversion( pe = equilibrium_solver(inv_input; override_psi_nodes) - # Round-trip validation: (ψ,θ) → (R,Z) → ψ_spline − ψ_target. - # Checks 4 angles at the outermost surface and at 75% of the radial grid. - # A large error indicates that Contour.jl resolution was insufficient near the - # x-point and the traced surface positions are inaccurate. - max_rt_err = 0.0 - for ψ_check in (pe.rzphi_xs[end], pe.rzphi_xs[max(1, length(pe.rzphi_xs) * 3 ÷ 4)]) - for θ_check in (0.0, 0.25, 0.5, 0.75) + # Round-trip validation: (ψ,θ) → (R,Z) → ψ_spline − ψ_target, over 4 angles on the outer + # ROUNDTRIP_EDGE_FRAC of the radial grid. A large residual means Contour.jl resolution was + # insufficient near the x-point and the traced surface positions are inaccurate. + # + # Measured twice, at the ψ knots and at the knot midpoints, because the rzphi splines + # interpolate exactly at their own knots: an on-knot residual cannot see inter-knot ringing + # at all, and stays flat at ~1e-6 on a psihigh ladder whose stability energies degrade by + # six orders of magnitude. The midpoint residual tracks that degradation, and their RATIO is + # the usable signal — dimensionless, so it needs no per-grid calibration. + rt_residual(ψ_samples) = maximum( + begin r2 = pe.rzphi_rsquared((ψ_check, θ_check)) off = pe.rzphi_offset((ψ_check, θ_check)) rfac = sqrt(max(r2, 0.0)) η = 2π * (θ_check + off) - R = pe.ro + rfac * cos(η) - Z = pe.zo + rfac * sin(η) - ψ_rt = 1.0 - raw_profile.psi_in((R, Z)) / psio - max_rt_err = max(max_rt_err, abs(ψ_rt - ψ_check)) + abs((1.0 - raw_profile.psi_in((pe.ro + rfac * cos(η), pe.zo + rfac * sin(η))) / psio) - ψ_check) end - end - if max_rt_err > 2e-3 - @warn "efit_by_inversion: round-trip error at edge = $(@sprintf("%.2e", max_rt_err)) > 2e-3; accuracy near psihigh may be limited. Consider reducing psihigh or increasing resolution_factor." + for ψ_check in ψ_samples, θ_check in (0.0, 0.25, 0.5, 0.75)) + + xs = pe.rzphi_xs + knots = @view xs[max(1, ceil(Int, length(xs) * ROUNDTRIP_EDGE_FRAC)):end] + max_rt_err = rt_residual(knots) + mids = [(knots[i] + knots[i+1]) / 2 for i in 1:(length(knots)-1)] + rt_mid = isempty(mids) ? NaN : rt_residual(mids) + rt_ratio = rt_mid / max(max_rt_err, eps()) + + # The ratio only means something once the midpoint residual is above the noise it would + # otherwise be measuring: on a clean grid both residuals sit near 1e-6, where their ratio is + # dominated by rounding rather than by ringing. + ringing = rt_mid > ROUNDTRIP_RINGING_FLOOR && rt_ratio > ROUNDTRIP_RATIO_TOL + too_large = max_rt_err > ROUNDTRIP_TOL || rt_mid > ROUNDTRIP_TOL + msg = + "efit_by_inversion: round-trip error at edge = $(@sprintf("%.2e", max_rt_err)) on knots, " * + "$(@sprintf("%.2e", rt_mid)) at knot midpoints (ratio $(@sprintf("%.1f", rt_ratio)))" + if too_large + @warn "$msg; exceeds $(ROUNDTRIP_TOL), so accuracy near psihigh may be limited. Consider reducing psihigh or increasing resolution_factor." + elseif ringing + @warn "$msg; the midpoint residual is $(@sprintf("%.0f", rt_ratio))x the on-knot one, so the rzphi splines are ringing between knots even though both residuals are small. Consider reducing psihigh or increasing resolution_factor." else - @info "efit_by_inversion: round-trip error at edge = $(@sprintf("%.2e", max_rt_err))" + @info msg end return pe diff --git a/test/runtests.jl b/test/runtests.jl index 7845eea60..f71414232 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -1,4 +1,5 @@ using Test +using Logging using Pkg using GeneralizedPerturbedEquilibrium.Vacuum using GeneralizedPerturbedEquilibrium.Equilibrium diff --git a/test/runtests_equil.jl b/test/runtests_equil.jl index af3cd5743..2bb25f40c 100644 --- a/test/runtests_equil.jl +++ b/test/runtests_equil.jl @@ -68,6 +68,39 @@ @test all(>(0), B_nodes) end + @testset "Round-trip check sees inter-knot ringing" begin + EQ = GeneralizedPerturbedEquilibrium.Equilibrium + # The rzphi splines interpolate exactly at their own knots, so an on-knot residual is + # blind to ringing between them. The check therefore samples the knot midpoints too and + # compares the two, and the ringing verdict is not consulted while the midpoint residual + # is still at the rounding floor, where the ratio would measure noise. + @test EQ.ROUNDTRIP_RINGING_FLOOR < EQ.ROUNDTRIP_TOL + @test EQ.ROUNDTRIP_RATIO_TOL > 1 + @test 0 < EQ.ROUNDTRIP_EDGE_FRAC < 1 + + # False-positive control: a deck that traces cleanly must report both residuals and stay + # at Info. The DIII-D-like geqdsk is used because the CHEASE deck above is the case that + # already trips the pre-existing absolute tolerance on this path. + cfg = EQ.EquilibriumConfig(; + eq_filename=joinpath(@__DIR__, "..", "examples", "DIIID-like_ideal_example", "TkMkr_D3Dlike_Hmode.geqdsk"), + eq_type="efit_by_inversion", + jac_type="hamada", + grid_type="log_asymptotic", + psilow=1e-4, + psihigh=0.995, + mpsi=128, + mtheta=128 + ) + logs, _ = Test.collect_test_logs(; min_level=Logging.Info) do + EQ.setup_equilibrium(cfg) + end + rt = filter(l -> occursin("round-trip error at edge", string(l.message)), logs) + @test length(rt) == 1 + @test rt[1].level == Logging.Info + @test occursin("at knot midpoints", string(rt[1].message)) + @test occursin("ratio", string(rt[1].message)) + end + @testset "Resolved psihigh" begin # The config holds the user's request and is never written to; the value the # equilibrium is actually formed on rides on params.psihigh_resolved.