FFS - EXPERIMENT - Implement route (a) and measure: -42% steps, ratio 1.92 -> 1.59

Dense output plus a bracketed Brent root-solve for the eta where the normalised SFL
angle hits each target node, so ff_fs_nodes is built directly on the uniform theta
grid and every surface is sampled at the same abscissae. Arclength tracer falls back.

Accepted steps 2309/3977/7638 -> 1881/2768/4403. Geometry residuals now converge and
every r1 flips positive; mid-plasma C goes 1.21e-5 -> 1.86e-6 with r1 -0.401 ->
+0.856. Delta' diagonal and singular surfaces preserved.

Two honest negatives recorded: the etol self-consistency prediction failed (still
flat at 2857/2768/2817), and the crude dtmax probe still reaches 3614 vs 4403 despite
matching route (a) on every available diagnostic -- an 18% gap with no measurable
difference in the integrand, flagged as the open question.

Implementation saved as route_a.patch; src reverted.

Co-Authored-By: Claude Opus 5 (1M context) <noreply@anthropic.com>

Author: logan-nc

handoff/issue376/RESULTS.md

@@ -602,6 +602,56 @@ still leaves finite resample error. The Phase 2 fix evaluates *exactly* at the t
 the resample error at the output nodes is zero rather than merely small, and it does not force
 extra steps.
 
+## 16. Route (a) implemented and measured (`route_a.patch`)
+
+Implementation: `direct_fieldline_int` solves with `dense=true` and returns the solution;
+`equilibrium_solver` then root-solves (Brent, bracketed by the monotone `y_out[:,5]`) for the η
+where the normalised SFL angle hits each target node and evaluates there, so `ff_fs_nodes` is built
+**directly on the uniform θ grid** and every surface is sampled at the same abscissae. The
+arclength tracer returns `nothing` and falls back to the old path. Patch saved to
+`route_a.patch`; `src/` is reverted (Phase 4 applies it on a branch off `develop`).
+
+### Primary — accepted steps
+
+| mpsi | stock | route (a) | change |
+|---|---|---|---|
+| 256  | 2309 | 1881 | −19% |
+| 512  | 3977 | 2768 | −30% |
+| 1024 | 7638 | **4403** | **−42%** |
+
+Per-doubling ratio **1.72 / 1.92 → 1.47 / 1.59**. This clears the plan's "strong win" threshold
+(≤5350 at mpsi=1024) and is well past the §11 ceiling that the partial repairs implied.
+
+### Diagnostics
+
+- **Geometry**: near-axis `offset` residual 6.39e-4 → **3.48e-6** at mpsi=256, and every `nu` /
+  `offset` / `rcoords` r1 flips positive (e.g. near-axis `nu` −0.345 → **+0.938**). Across the
+  ladder the residuals now **converge**: mid-plasma `nu` 2.65e-7 → 6.36e-8 → 6.57e-9.
+- **Matrices**: mid-plasma at mpsi=1024, C 1.21e-5 → **1.86e-6** with r1 −0.401 → **+0.856**;
+  G 1.70e-5 → 1.29e-6 with r1 −0.176 → **+0.947**. A is unchanged, as expected.
+- **Physics preserved**: Δ′ diagonal at mpsi=1024 [8.58, −5.53, −16.09, −2368.85, 55.11] →
+  [8.59, −5.51, −16.08, −2367.99, 55.62]; et[1] 0.802840 → 0.799037 (4.7e-3 relative, an accuracy
+  change from sampling the geometry correctly). Singular surfaces and `psilim` identical.
+
+### Two results that say the job is not finished
+
+1. **The etol self-consistency prediction failed.** The plan predicted that once the remap error
+   was gone, ε would become integration-limited and the `etol` sweep would stop being null. At
+   mpsi=512 after route (a): **2857 / 2768 / 2817** accepted steps for etol 1e-8 / 1e-10 / 1e-12 —
+   still flat. A residual ε remains that is not controlled by integration tolerance.
+2. **The crude `dtmax` probe still beats it** — 3614 vs 4403 at mpsi=1024 (ratio 1.43 vs 1.59) —
+   even though every diagnostic available says the two runs are the same: geometry residuals equal
+   to 3-4 significant figures (mid `nu` 6.567e-9 vs 6.499e-9), matrix residuals equal (C 1.86e-6
+   vs 1.84e-6), profiles equal, singular surfaces and `psilim` identical. Elementwise the two
+   geometries differ by only ~1e-6 relative, against the ~1e-4 by which both differ from stock.
+
+**Open question for the next session**: an 18% step difference with no measurable difference in the
+integrand. Whatever it is sits outside the node values, the coefficient matrices, the profiles and
+the surface positions — the leading suspect is the dense-output interpolant's behaviour *between*
+solver steps (route (a) evaluates it; `dtmax` makes the steps small enough that it barely matters),
+which the current diagnostics, all built on node values, cannot see. Closing that gap is worth
+roughly another 18% on top of route (a)'s 42%.
+
 ## Implications / ranked follow-ups
 
 1. **Use the two-pass auto grid** (`mpsi=0`, `psi_accuracy`) — already the example default; it

handoff/issue376/route_a.patch

@@ -0,0 +1,128 @@
+diff --git a/src/Equilibrium/DirectEquilibrium.jl b/src/Equilibrium/DirectEquilibrium.jl
+index ed6582131..c3c9ccfdd 100644
+--- a/src/Equilibrium/DirectEquilibrium.jl
++++ b/src/Equilibrium/DirectEquilibrium.jl
+@@ -227,6 +227,32 @@ function direct_position!(raw_profile::DirectRunInput)
+     return ro, zo, rs1, rs2
+ end
+ 
++"""
++    eta_at_sfl_angle(sol, y_out, x, total_x) -> Float64
++
++Integration angle η at which the normalised straight-fieldline angle ∫jac·dl/Bp reaches `x`.
++
++`y_out[:, 5]` is monotone in η, so it brackets the root to one solver step and Brent converges in
++a handful of dense-output evaluations. Used to sample every flux surface at the *same* SFL angles
++instead of resampling each surface's own solver steps (issue #376).
++"""
++function eta_at_sfl_angle(sol, y_out::Matrix{Float64}, x::Float64, total_x::Float64)
++    x <= 0 && return y_out[1, 1]
++    x >= 1 && return y_out[end, 1]
++    target = x * total_x
++    hi = searchsortedfirst(view(y_out, :, 5), target)
++    hi = clamp(hi, 2, size(y_out, 1))
++    lo = hi - 1
++    eta_lo, eta_hi = y_out[lo, 1], y_out[hi, 1]
++    f(eta) = sol(eta)[4] - target
++    flo, fhi = f(eta_lo), f(eta_hi)
++    # Degenerate bracket (repeated η, or the root sitting exactly on a step) needs no solve.
++    flo == 0 && return eta_lo
++    fhi == 0 && return eta_hi
++    (flo * fhi > 0 || eta_hi <= eta_lo) && return eta_lo + (eta_hi - eta_lo) * (target - y_out[lo, 5]) / max(y_out[hi, 5] - y_out[lo, 5], eps())
++    return find_zero(f, (eta_lo, eta_hi), Roots.Brent())
++end
++
+ """
+     direct_fieldline_int(psifac, raw_profile, ro, zo, rs2)
+ 
+@@ -250,9 +276,12 @@ from 1:5 rather than 0:4 as in Fortran.
+         - `y_out[:, 4]`: ∫(dl/(R²Bp))
+         - `y_out[:, 5]`: ∫(jac*dl/Bp)
+ 
++  - `sol`: the dense ODE solution, so callers can evaluate the trace at prescribed SFL angles
++    rather than resampling this surface's own solver steps (`nothing` for tracers without it).
++
+   - `bfield`: A `DirectBField` object with values at the integration start point.
+ """
+-function direct_fieldline_int(psifac::Float64, raw_profile::DirectRunInput, ro::Float64, zo::Float64, rs2::Float64)::Tuple{Matrix{Float64},DirectBField}
++function direct_fieldline_int(psifac::Float64, raw_profile::DirectRunInput, ro::Float64, zo::Float64, rs2::Float64)
+ 
+     # Find the starting point on the flux surface (outboard midplane)
+     psi0_guess = raw_profile.psio * (1.0 - psifac)
+@@ -291,10 +320,12 @@ function direct_fieldline_int(psifac::Float64, raw_profile::DirectRunInput, ro::
+     callback = DiscreteCallback((u, t, i) -> true, refine_affect!; save_positions=(true, false))
+ 
+     prob = ODEProblem{true}(direct_fieldline_der!, u0, (0.0, 2π), params)
+-    sol = solve(prob, Vern9(); callback=callback, reltol=equil_config.etol, abstol=1e-8, dt=2π / 200, adaptive=true, dense=false)
++    # Dense output lets the caller evaluate the trace at the SFL angles it actually wants, instead
++    # of splining this surface's solver-chosen steps and resampling (issue #376).
++    sol = solve(prob, Vern9(); callback=callback, reltol=equil_config.etol, abstol=1e-8, dt=2π / 200, adaptive=true, dense=true)
+ 
+     sol_matrix = reduce(hcat, sol.u::Vector{Vector{Float64}})'
+-    return hcat(sol.t::Vector{Float64}, sol_matrix), bfield
++    return hcat(sol.t::Vector{Float64}, sol_matrix), bfield, sol
+ end
+ 
+ """
+@@ -495,18 +526,39 @@ robustness.
+     ff_deriv_val = zeros!(pool, Float64, 4)
+ 
+     for ipsi in (mpsi+1):-1:1  # outermost to innermost
+-        y_out, bfield = fieldline_int(psi_nodes[ipsi], raw_profile, ro, zo, rs2)
++        y_out, bfield, sol = fieldline_int(psi_nodes[ipsi], raw_profile, ro, zo, rs2)
+         checkpoint!(pool, Float64)
+ 
+-        # Fit data into temporary straight fieldline poloidal angle splines
+-        ff_x_nodes = acquire!(pool, Float64, size(y_out, 1))
+-        @. ff_x_nodes = @view(y_out[:, 5]) / y_out[end, 5]
+-
+-        ff_fs_nodes = acquire!(pool, Float64, size(y_out, 1), 4)
+-        @. ff_fs_nodes[:, 1] = @view(y_out[:, 3])^2
+-        @. ff_fs_nodes[:, 2] = @view(y_out[:, 1]) / (2π) - ff_x_nodes
+-        @. ff_fs_nodes[:, 3] = bfield.f * (@view(y_out[:, 4]) - ff_x_nodes * y_out[end, 4])
+-        @. ff_fs_nodes[:, 4] = @view(y_out[:, 2]) / y_out[end, 2] - ff_x_nodes
++        # Straight-fieldline angle x = normalised ∫jac·dl/Bp, monotone in the integration angle η.
++        #
++        # Sampling x at this surface's own solver steps and resampling onto theta_nodes leaves a
++        # resample error that is uncorrelated between neighbouring surfaces, i.e. white noise in ψ
++        # that grid refinement then amplifies (issue #376). With dense output we instead solve for
++        # the η where x hits each target node and evaluate there, so every surface is sampled at
++        # the same abscissae and the resample error at the output nodes is zero.
++        nff = sol === nothing ? size(y_out, 1) : mtheta + 1
++        ff_x_nodes = acquire!(pool, Float64, nff)
++        ff_fs_nodes = acquire!(pool, Float64, nff, 4)
++
++        if sol === nothing
++            @. ff_x_nodes = @view(y_out[:, 5]) / y_out[end, 5]
++            @. ff_fs_nodes[:, 1] = @view(y_out[:, 3])^2
++            @. ff_fs_nodes[:, 2] = @view(y_out[:, 1]) / (2π) - ff_x_nodes
++            @. ff_fs_nodes[:, 3] = bfield.f * (@view(y_out[:, 4]) - ff_x_nodes * y_out[end, 4])
++            @. ff_fs_nodes[:, 4] = @view(y_out[:, 2]) / y_out[end, 2] - ff_x_nodes
++        else
++            total_x = y_out[end, 5]
++            for itheta in 1:(mtheta+1)
++                x = theta_nodes[itheta]
++                eta = eta_at_sfl_angle(sol, y_out, x, total_x)
++                u = sol(eta)
++                ff_x_nodes[itheta] = x
++                ff_fs_nodes[itheta, 1] = u[2]^2
++                ff_fs_nodes[itheta, 2] = eta / (2π) - x
++                ff_fs_nodes[itheta, 3] = bfield.f * (u[3] - x * y_out[end, 4])
++                ff_fs_nodes[itheta, 4] = u[1] / y_out[end, 2] - x
++            end
++        end
+ 
+         ff_fs_nodes[end, :] .= ff_fs_nodes[1, :]  # enforce periodic endpoint
+ 
+diff --git a/src/Equilibrium/DirectEquilibriumArcLength.jl b/src/Equilibrium/DirectEquilibriumArcLength.jl
+index 32b93a0ac..6774e5601 100644
+--- a/src/Equilibrium/DirectEquilibriumArcLength.jl
++++ b/src/Equilibrium/DirectEquilibriumArcLength.jl
+@@ -142,6 +142,6 @@ outboard midplane (Z = zo, R > ro) after a minimum arc-length guard.
+     end
+ 
+     # bfield at the starting point carries F and P for the surface-averaged quantities
+-    return y_out, bfield
++    return y_out, bfield, nothing
+ end
+