From de0db4f551b4774df959109e146dfffd6a5532b3 Mon Sep 17 00:00:00 2001 From: logan-nc Date: Sun, 16 Aug 2026 18:28:59 -0400 Subject: [PATCH] EQUIL - IMPROVEMENT - Sample every flux surface at the same straight-fieldline angles Each surface was traced independently and then splined on that surface's OWN solver-chosen abscissae before being resampled onto the common theta grid, so the resample error was uncorrelated between neighbouring surfaces -- white noise in psi that grid refinement amplifies rather than reduces. Neither reltol nor abstol touched it, because it is remap interpolation error rather than integration error. The trace now returns its dense solution, and equilibrium_solver root-solves (Brent, bracketed by the monotone jac-weighted flux integral) for the angle at which the normalised straight-fieldline angle reaches each target node, evaluating there. Every surface is sampled at identical abscissae and the resample error at the output nodes is zero. The arclength tracer returns nothing for the solution and keeps the previous path. Measured on DIII-D stripped decks at eulerlagrange_tolerance 1e-10, accepted Euler-Lagrange steps fall 2309/3977/7638 -> 1881/2768/4403 for mpsi 256/512/1024, a 42% reduction at mpsi=1024, and the per-doubling growth drops from 1.92x to 1.59x. Surface geometry residuals now converge with refinement instead of sitting on a floor, and every knot-to-knot correlation turns positive. Co-Authored-By: Claude Opus 5 (1M context) --- src/Equilibrium/DirectEquilibrium.jl | 78 +++++++++++++++---- src/Equilibrium/DirectEquilibriumArcLength.jl | 4 +- 2 files changed, 67 insertions(+), 15 deletions(-) diff --git a/src/Equilibrium/DirectEquilibrium.jl b/src/Equilibrium/DirectEquilibrium.jl index 961a0f7f5..586d81001 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..ec4c661d3 100644 --- a/src/Equilibrium/DirectEquilibriumArcLength.jl +++ b/src/Equilibrium/DirectEquilibriumArcLength.jl @@ -75,7 +75,7 @@ outboard midplane (Z = zo, R > ro) after a minimum arc-length guard. """ @with_pool pool function arclength_fieldline_int( psifac::Float64, raw_profile::DirectRunInput, ro::Float64, zo::Float64, rs2::Float64 -)::Tuple{Matrix{Float64},DirectBField} +)::Tuple{Matrix{Float64},DirectBField,Nothing} psi0_guess = raw_profile.psio * (1.0 - psifac) r = ro + sqrt(psifac) * (rs2 - ro) @@ -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