Equilibrium - BUGFIX! - Sample all flux surfaces on common straight-fieldline abscissae and cap the coefficient-spline core knots

docs/src/kinetic_forces.md

@@ -138,6 +138,35 @@ terms respectively, and do not modify the stored kinetic profile splines.
        when computing the `toroidal_rotation_factor` back-solve. Julia uses a
        clean reimplementation with consistent pre-scaling derivatives throughout.
 
+## Regularization: why kinetic runs set `reg_spot = 0`
+
+`[PerturbedEquilibrium] reg_spot` smooths the displacements before they drive the NTV
+integrand, multiplying ``\xi^{\psi\prime}`` and ``\xi^\alpha`` by
+``Q^2/(Q^2 + \mathrm{reg\_spot}^2)`` with ``Q = m - nq``. It exists because **ideal** MHD is
+singular at the rationals: ``\bar F_\mathrm{ideal} = Q F Q`` has ``\det \bar F = 0`` there, so
+those two components diverge as ``1/Q`` and the torque integral does not converge.
+
+The **self-consistent kinetic** workflow has no such singularity. Park & Logan
+([Phys. Plasmas 24, 032505 (2017)](https://doi.org/10.1063/1.4978562), §III D) decompose the
+kinetic composite matrix as ``F_k = Q \bar F_k Q - P_l^\dagger Q - Q P_u + R_1`` where
+``R_1 \neq 0`` at ``Q = 0``; with finite torque ``\det \bar F`` is complex, its zeros leave the
+real ``\psi`` axis, and the singularity is removed from both the solution and the torque
+integral. (Torque-free kinetic energy principles instead *shift and split* the zeros to
+``\psi_r \mp r_{L,R}``, where the singularity is logarithmic and integrable — still not a case
+for smoothing.) Logan, *Electromagnetic Torque in Tokamaks with Toroidal Asymmetries* (PhD thesis,
+Princeton, 2015), Ch. 7 Eq. 7.46 states the same resonance shifting/splitting result in GPEC's own
+notation.
+
+GPEC therefore **forces `reg_spot = 0` whenever `kinetic_factor > 0`**, logging the override.
+Leaving it on suppresses a finite physical response and does so inconsistently — ``\xi^\psi``
+is never regularized, so damping the other two breaks their near-resonance cancellation in
+``\delta B/B`` and leaves a spurious residue driving the NTV integrand.
+
+The check that this is right is that the NTV torque and the Euler–Lagrange solution's own
+dissipation ``-2n\,\mathrm{Im}\langle \xi, u_2\rangle/4\mu_0`` — two independent calculations
+of the same quantity — agree to a fraction of a percent with the regularization off, and disagree
+by tens of percent with it on.
+
 ## HDF5 outputs: complex torque convention and the EnergyIntegrals layout
 
 The method level of `KineticForces/<method>/` reports the two physical scalars a user

examples/DIIID-like_ideal_example/gpec.toml

@@ -94,5 +94,6 @@ f0type = "maxwellian"                   # Equilibrium distribution
 moment = "pressure"                     # Pressure-moment NTV torque
 atol_xlmda = 1e-9                       # Absolute tolerance for inner pitch + energy integrations
 rtol_xlmda = 1e-5                       # Relative tolerance for inner pitch + energy integrations
+axis_validity_suppression = true        # Suppress kinetic terms inside the near-axis orbit-width validity limit
 write_outputs_to_HDF5 = true            # Write outputs to the HDF5 file
 verbose = true                          # Enable verbose logging

examples/Solovev_kinetic_NTV_example/gpec.toml

@@ -87,3 +87,4 @@ f0fac = 1      # Scale toroidal field at constant pressure (β, q change; Φ, p,
 
 [KineticForces]
 kinetic_file = "kinetic.dat"            # Kinetic profile file: psi_n, n_i, n_e, T_i, T_e, omega_E columns
+axis_validity_suppression = true        # Suppress kinetic terms inside the near-axis orbit-width validity limit

examples/a10_kinetic_example/gpec.toml

@@ -59,3 +59,4 @@ nutype = "harmonic"            # Collision operator (zero, small, krook, harmoni
 f0type = "maxwellian"          # Distribution function (maxwellian, jkp, cgl)
 atol_xlmda = 1e-9              # Absolute tolerance for inner pitch + energy integrations
 rtol_xlmda = 1e-5              # Relative tolerance for inner pitch + energy integrations
+axis_validity_suppression = true   # Suppress kinetic terms inside the near-axis orbit-width validity limit

src/Equilibrium/DirectEquilibrium.jl

@@ -228,7 +228,51 @@ function direct_position!(raw_profile::DirectRunInput)
 end
 
 """
-    direct_fieldline_int(psifac, raw_profile, ro, zo, rs2)
+    sfl_sample_time(sol, ts, sfl, sfl_index, x) -> (t, lo)
+
+Solver time `t` at which the normalised straight-fieldline angle reaches `x` ∈ [0, 1], and the index
+`lo` of the saved step that starts its bracket. `sfl` is the SFL integral ∫jac·dl/Bp at the saved
+times `ts`; it is monotone, so it brackets the root to one step, and Brent finishes on component
+`sfl_index` of the dense solution.
+"""
+function sfl_sample_time(sol, ts::Vector{Float64}, sfl::AbstractVector{Float64}, sfl_index::Int, x::Float64)
+    n = length(ts)
+    x <= 0 && return ts[1], 1
+    x >= 1 && return ts[n], n - 1
+    target = x * sfl[n]
+    hi = clamp(searchsortedfirst(sfl, target), 2, n)
+    lo = hi - 1
+    flo, fhi = sfl[lo] - target, sfl[hi] - target
+    flo == 0 && return ts[lo], lo
+    fhi == 0 && return ts[hi], lo
+    flo < 0 < fhi || error("sfl_sample_time: SFL integral is not monotone across a solver step (x = $x)")
+    return find_zero(t -> sol(t)[sfl_index] - target, (ts[lo], ts[hi]), Roots.Brent()), lo
+end
+
+"""
+    sample_trace_at_sfl_angles(sol, y_out, sfl_index, theta_nodes, state_row) -> Matrix{Float64}
+
+Sample a traced flux surface at the straight-fieldline angles `theta_nodes` (0 to 1). Row `i` is
+`state_row(sol(t), t, lo)` at the time `t` where the normalised SFL angle reaches `theta_nodes[i]`
+(see `sfl_sample_time`), in the 5-column `y_out` layout. The last row is the saved closed endpoint,
+so the surface totals `y_out[end, :]` are exact. Every surface sharing the same abscissae keeps the
+resample error from varying between neighbouring surfaces. Every field-line tracer must return its
+surface through this function.
+"""
+function sample_trace_at_sfl_angles(sol, y_out::Matrix{Float64}, sfl_index::Int, theta_nodes, state_row::F) where {F}
+    ts = sol.t::Vector{Float64}
+    sfl = view(y_out, :, 5)
+    samples = Matrix{Float64}(undef, length(theta_nodes), 5)
+    for i in 1:(length(theta_nodes)-1)
+        t, lo = sfl_sample_time(sol, ts, sfl, sfl_index, Float64(theta_nodes[i]))
+        samples[i, :] .= state_row(sol(t), t, lo)
+    end
+    samples[end, :] .= @view y_out[end, :]
+    return samples
+end
+
+"""
+    direct_fieldline_int(psifac, raw_profile, ro, zo, rs2, theta_nodes) -> (y_out, bfield)
 
 Performs the field-line integration for a single flux surface. This is a Julia adaptation
 of the Fortran `direct_fl_int` subroutine. Note that the array `y_out` is now indexed
@@ -240,10 +284,12 @@ from 1:5 rather than 0:4 as in Fortran.
   - `raw_profile`: `DirectRunInput` object containing splines and parameters.
   - `ro`, `zo`: Coordinates of the magnetic axis [m].
   - `rs2`: R-coordinate of the outboard separatrix [m].
+  - `theta_nodes`: straight-fieldline angles (0 to 1) at which to sample the surface.
 
 ## Returns:
 
-    - `y_out`: A matrix containing the integrated quantities vs. the geometric angle `η`.
+    - `y_out`: A matrix containing the integrated quantities at each of `theta_nodes`
+      (see `sample_trace_at_sfl_angles`).
         - `y_out[:, 1]`: η (geometric poloidal angle)
         - `y_out[:, 2]`: ∫(dl/Bp)
         - `y_out[:, 3]`: rfac (radial distance from magnetic axis)
@@ -252,7 +298,8 @@ from 1:5 rather than 0:4 as in Fortran.
 
   - `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,
+    theta_nodes)::Tuple{Matrix{Float64},DirectBField}
 
     # Find the starting point on the flux surface (outboard midplane)
     psi0_guess = raw_profile.psio * (1.0 - psifac)
@@ -277,10 +324,11 @@ function direct_fieldline_int(psifac::Float64, raw_profile::DirectRunInput, ro::
     u0 = zeros(Float64, 4)
     u0[2] = sqrt((r - ro)^2 + (z - zo)^2)
 
-    bfield = DirectBField()
+    # The ODE right-hand side overwrites its own field buffer, so keep `bfield` at the start point (F, P at ψ₀).
+    bfield_ode = DirectBField()
     equil_config = raw_profile.config
     params = FieldLineDerivParams(ro, zo, raw_profile.psi_in, raw_profile.sq_in, sq_in_deriv, raw_profile.psio,
-        equil_config.power_bp, equil_config.power_b, equil_config.power_r, equil_config.power_rc, bfield)
+        equil_config.power_bp, equil_config.power_b, equil_config.power_r, equil_config.power_rc, bfield_ode)
 
     # Use a callback to refine the solution at each step to stay on the flux surface
     function refine_affect!(integrator)
@@ -291,7 +339,8 @@ 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 so the surface can be sampled at prescribed SFL angles between solver steps.
+    sol = solve(prob, Vern9(); callback=callback, reltol=equil_config.etol, abstol=1e-8, dt=2π / 200, adaptive=true, dense=true)
 
     # A failed solve returns a truncated solution instead of throwing; check both the retcode and
     # that the field line reached η = 2π, since a callback can end it early and still report Success.
@@ -312,7 +361,9 @@ function direct_fieldline_int(psifac::Float64, raw_profile::DirectRunInput, ro::
     end
 
     sol_matrix = reduce(hcat, sol.u::Vector{Vector{Float64}})'
-    return hcat(sol.t::Vector{Float64}, sol_matrix), bfield
+    y_out = hcat(sol.t::Vector{Float64}, sol_matrix)
+    # State is [∫dl/Bp, rfac, ∫dl/(R²Bp), ∫jac·dl/Bp] and the integration variable is η itself.
+    return sample_trace_at_sfl_angles(sol, y_out, 4, theta_nodes, (u, eta, _) -> (eta, u[1], u[2], u[3], u[4])), bfield
 end
 
 """
@@ -485,6 +536,9 @@ robustness.
 
   - `raw_profile`: A `DirectRunInput` object containing the initial splines (`psi_in`, `sq_in`)
     and run parameters (`equil_input`).
+  - `fieldline_int`: the flux-surface tracer, called as `fieldline_int(psifac, raw_profile, ro, zo,
+    rs2, theta_nodes) -> (y_out, bfield)` with `y_out` already sampled at `theta_nodes` through
+    `sample_trace_at_sfl_angles`.
 
 ## Returns:
 
@@ -510,20 +564,18 @@ robustness.
     mpsi = length(psi_nodes) - 1
     theta_nodes = range(0.0, 1.0; length=mtheta + 1)
 
+    ff_x_nodes = collect(theta_nodes)
     sq_nodes = zeros!(pool, Float64, mpsi + 1, 4)
     rzphi_nodes = zeros!(pool, Float64, mpsi + 1, mtheta + 1, 4)
     ff_val = zeros!(pool, Float64, 4)
     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)
+        # Every tracer returns the surface sampled at theta_nodes, so all surfaces share SFL abscissae.
+        y_out, bfield = fieldline_int(psi_nodes[ipsi], raw_profile, ro, zo, rs2, theta_nodes)
         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 = acquire!(pool, Float64, mtheta + 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])

src/Equilibrium/DirectEquilibriumArcLength.jl

@@ -59,10 +59,11 @@ so the solver never restricts step size for them.
 end
 
 """
-    arclength_fieldline_int(psifac, raw_profile, ro, zo, rs2)
+    arclength_fieldline_int(psifac, raw_profile, ro, zo, rs2, theta_nodes) -> (y_out, bfield)
 
 Arc-length-parameterized flux surface integration. Drop-in replacement for
-`direct_fieldline_int` with identical return format:
+`direct_fieldline_int` with identical return format, sampled at the straight-fieldline angles
+`theta_nodes` (see `sample_trace_at_sfl_angles`):
 
 - `y_out[:, 1]`: geometric angle η ∈ 0 to 2π (CCW from outboard midplane)
 - `y_out[:, 2]`: accumulated ∫dl/Bp
@@ -74,7 +75,7 @@ The ODE is terminated by a `ContinuousCallback` that detects the return to the
 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
+    psifac::Float64, raw_profile::DirectRunInput, ro::Float64, zo::Float64, rs2::Float64, theta_nodes
 )::Tuple{Matrix{Float64},DirectBField}
 
     psi0_guess = raw_profile.psio * (1.0 - psifac)
@@ -114,13 +115,14 @@ outboard midplane (Z = zo, R > ro) after a minimum arc-length guard.
     reltol_vec = [equil_config.etol, equil_config.etol, 1e20, 1e20, 1e20]
     abstol_vec = [1e-8, 1e-8, 1e20, 1e20, 1e20]
     sol = solve(prob, BS5(); callback=callback, reltol=reltol_vec, abstol=abstol_vec,
-        dt=2π / 200, adaptive=true, dense=false)
+        dt=2π / 200, adaptive=true, dense=true)
 
     n = length(sol.u)
     y_out = Matrix{Float64}(undef, n, 5)
 
     # Compute geometric angles and unwrap to monotone [0, 2π]
-    η_raw = [atan(sol.u[i][2] - zo, sol.u[i][1] - ro) for i in 1:n]
+    η_atan = [atan(sol.u[i][2] - zo, sol.u[i][1] - ro) for i in 1:n]
+    η_raw = copy(η_atan)
     for i in 2:n
         Δ = η_raw[i] - η_raw[i-1]
         if Δ > π
@@ -141,7 +143,13 @@ outboard midplane (Z = zo, R > ro) after a minimum arc-length guard.
         y_out[i, 5] = sol.u[i][5]
     end
 
+    # The integration variable is arc length, so η at a sample is unwrapped against its bracketing step.
+    function state_row(u, _s, lo)
+        η = y_out[lo, 1] + rem2pi(atan(u[2] - zo, u[1] - ro) - η_atan[lo], RoundNearest)
+        return (η, u[3], sqrt((u[1] - ro)^2 + (u[2] - zo)^2), u[4], u[5])
+    end
+
     # bfield at the starting point carries F and P for the surface-averaged quantities
-    return y_out, bfield
+    return sample_trace_at_sfl_angles(sol, y_out, 5, theta_nodes, state_row), bfield
 end
 

src/Equilibrium/GridRefinement.jl

@@ -419,6 +419,7 @@ 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;
+`pinned` lists ψ values inserted as plain knots without a cleared zone (kinetic-resonance surfaces);
 `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
@@ -428,6 +429,7 @@ function refined_psi_grid(equil::PlasmaEquilibrium;
     tau::Float64,
     kin::Union{Nothing,KineticProfileSplines}=nothing,
     mandatory::Vector{Float64}=Float64[],
+    pinned::Vector{Float64}=Float64[],
     singfac_min::Float64=1e-4,
     n_min::Int=1,
     bracket_coef::Float64=BRACKET_COEF,
@@ -449,6 +451,9 @@ function refined_psi_grid(equil::PlasmaEquilibrium;
     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 = enforce_min_spacing(_equidistribute(xs, rho, N), min_spacing)
+    # Pinned knots use knot-at-node semantics (no cleared zone); inserted before rational bracketing
+    # so a rational's bracket still wins locally where the two overlap.
+    grid = isempty(pinned) ? grid : merge_mandatory_nodes(grid, pinned)
     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)