From 2afc6c3ebe732b2f99083253c12a48bb9aef6c3f Mon Sep 17 00:00:00 2001 From: logan-nc Date: Sat, 15 Aug 2026 13:58:02 -0400 Subject: [PATCH] KineticForces - BUG FIX - fit the parallel-velocity spline periodically The bounce-point spline samples 1-(lmda/bo)*B(theta), which is periodic on the closed poloidal interval, but was fitted with non-periodic endpoint conditions while the adjacent equilibrium B spline (tspl) was already periodic. The resulting fit is only C0 at the theta=0/1 seam, so it can manufacture false near-seam extrema and bounce-root pairs. On a synthetic shaped field the seam derivative mismatch drops ~5x (1.4e-4 -> 2.6e-5) when the fit is made periodic. Fixed at both construction sites (the tpsi! quadrature path and the _setup_surface_state kinetic-matrix path). Renamed B_extrap -> B_vpar and _vpar_from_extrap -> _vpar_from_spline, since after this change the old names describe a fit that no longer exists. Julia analog of Fortran GPEC PR #284, which fixed the same spline_fit(vspl, "extrap") -> "periodic" defect in pentrc/torque.F90 after it produced non-finite omega_b/omega_D and an LSODE failure on an ITER case. Co-Authored-By: Claude Fable 5 --- src/KineticForces/BounceAveraging.jl | 38 ++++++++++++++-------------- src/KineticForces/Torque.jl | 25 ++++++++++-------- 2 files changed, 33 insertions(+), 30 deletions(-) diff --git a/src/KineticForces/BounceAveraging.jl b/src/KineticForces/BounceAveraging.jl index 63a59aedc..9db0033fe 100644 --- a/src/KineticForces/BounceAveraging.jl +++ b/src/KineticForces/BounceAveraging.jl @@ -147,7 +147,7 @@ end """ compute_bounce_data(psi, n, l, q, bo, bmax, bmin, theta_bmax, - tspl, B_extrap, mfac, chi1, ro, dbob_m_f, divx_m_f, + tspl, B_vpar, mfac, chi1, ro, dbob_m_f, divx_m_f, divxfac, wdfac, mass, chrg, T_s, method; nlmda=128, ntheta=128, smat=nothing, tmat=nothing, xmat=nothing, @@ -168,8 +168,8 @@ Ports Fortran torque.F90 lines 530-816 (GAR branch). - `bmax, bmin`: Max/min of B(θ) at this ψ - `theta_bmax`: θ location of Bmax (nodal knot; the passing-transit start) - `tspl`: Periodic poloidal interpolant: tspl(θ) → [B, dB/dψ, dB/dθ, J, dJ/dψ] -- `B_extrap`: Endpoint-fit (non-periodic) cubic of B(θ) used for v_par and the - bounce-point roots (the Fortran `vspl` equivalent) +- `B_vpar`: Periodic cubic of B(θ) used for v_par and the bounce-point roots + (the Fortran `vspl` equivalent) - `mfac`: Poloidal mode numbers [mlow:mhigh] - `chi1`: 2π·ψ₀ flux normalization - `ro`: Major radius [m] @@ -190,7 +190,7 @@ function compute_bounce_data( psi::Float64, n::Int, l::Int, q::Float64, bo::Float64, bmax::Float64, bmin::Float64, theta_bmax::Float64, - tspl, B_extrap, mfac::Vector{Int}, chi1::Float64, ro::Float64, + tspl, B_vpar, mfac::Vector{Int}, chi1::Float64, ro::Float64, dbob_m_f::Vector{ComplexF64}, divx_m_f::Vector{ComplexF64}, divxfac::Float64, wdfac::Float64, mass::Float64, chrg::Float64, @@ -241,12 +241,12 @@ function compute_bounce_data( # Find bounce points and build θ sub-grid _, _, tdt_pts, tdt_wts = _find_bounce_points_and_grid( - lmda, bo, sigma, B_extrap, theta_bmax, psi, ntheta) + lmda, bo, sigma, B_vpar, theta_bmax, psi, ntheta) # Bounce integrals over θ (Fortran lines 674-735) wbbar, wdbar, dJdJ_val, wmats_lmda = _bounce_integrate( tdt_pts, tdt_wts, lmda, lnq, sigma, n, q, bo, - tspl, B_extrap, chi1, ro, mfac, dbob_m_f, divx_m_f, divxfac, wdfac, + tspl, B_vpar, chi1, ro, mfac, dbob_m_f, divx_m_f, divxfac, wdfac, do_matrices, mpert, smat, tmat, xmat, ymat, zmat, scr) # Physical frequencies (Fortran lines 744-745) @@ -409,12 +409,12 @@ end """ -Parallel-velocity factor `v_par = 1 − (λ/bo)·B(θ)` from the endpoint-fit cubic of B -(`B_extrap`, built where the surface interpolants are constructed), keeping v_par +Parallel-velocity factor `v_par = 1 − (λ/bo)·B(θ)` from the periodic cubic of B +(`B_vpar`, built where the surface interpolants are constructed), keeping v_par consistent with the bounce-point roots as in Fortran's `vspl`. """ -@inline _vpar_from_extrap(B_extrap, lmda::Float64, bo::Float64, θ::Float64) = - 1.0 - (lmda / bo) * B_extrap(mod(θ, 1.0)) +@inline _vpar_from_spline(B_vpar, lmda::Float64, bo::Float64, θ::Float64) = + 1.0 - (lmda / bo) * B_vpar(mod(θ, 1.0)) """ @@ -423,14 +423,14 @@ Returns (t1, t2, theta_points, theta_weights). """ function _find_bounce_points_and_grid( lmda::Float64, bo::Float64, sigma::Int, - B_extrap, theta_bmax::Float64, psi::Float64, + B_vpar, theta_bmax::Float64, psi::Float64, ntheta::Int ) if sigma == 0 # trapped - # Bounce points: all roots of v_par(θ) = 1 − (λ/bo)·B_extrap(θ) in (0,1), + # Bounce points: all roots of v_par(θ) = 1 − (λ/bo)·B_vpar(θ) in (0,1), # sorted descending — the same order as Fortran spline_roots, which the # marginally-trapped and deepest-well wrap logic below assume. - vpar_fn = θ -> _vpar_from_extrap(B_extrap, lmda, bo, θ) + vpar_fn = θ -> _vpar_from_spline(B_vpar, lmda, bo, θ) bpts = sort!(Roots.find_zeros(vpar_fn, 0.0, 1.0); rev=true) nbpts = length(bpts) @@ -443,7 +443,7 @@ function _find_bounce_points_and_grid( t1 = bpts[1] t2 = bpts[1] + 1.0 else - t1, t2 = _find_deepest_well(bpts, B_extrap, lmda, bo) + t1, t2 = _find_deepest_well(bpts, B_vpar, lmda, bo) end # Power-law grid refined near bounce points @@ -463,7 +463,7 @@ end Find the deepest potential well (largest midpoint v_par) among bounce-point pairs, handling pairs that wrap through θ = 0/1. """ -function _find_deepest_well(bpts::Vector{Float64}, B_extrap, lmda::Float64, bo::Float64) +function _find_deepest_well(bpts::Vector{Float64}, B_vpar, lmda::Float64, bo::Float64) nbpts = length(bpts) best_vpar = 0.0 best_t1 = 0.0 @@ -477,7 +477,7 @@ function _find_deepest_well(bpts::Vector{Float64}, B_extrap, lmda::Float64, bo:: else θmid = 0.5 * (bpts[i] + bpts[j]) end - vpar_mid = _vpar_from_extrap(B_extrap, lmda, bo, θmid) + vpar_mid = _vpar_from_spline(B_vpar, lmda, bo, θmid) if vpar_mid > best_vpar best_t1 = bpts[i] best_t2 = bpts[j] @@ -506,7 +506,7 @@ Ports Fortran torque.F90 lines 674-793. function _bounce_integrate( tdt_pts::Vector{Float64}, tdt_wts::Vector{Float64}, lmda::Float64, lnq::Float64, sigma::Int, n::Int, q::Float64, bo::Float64, - tspl, B_extrap, chi1::Float64, ro::Float64, + tspl, B_vpar, chi1::Float64, ro::Float64, mfac::Vector{Int}, dbob_m_f::Vector{ComplexF64}, divx_m_f::Vector{ComplexF64}, divxfac::Float64, wdfac::Float64, do_matrices::Bool, mpert::Int, @@ -550,9 +550,9 @@ function _bounce_integrate( jac = tspl_f[4] djdpsi = tspl_f[5] - # v_par from the endpoint-fit cubic (consistent with the bounce points); + # v_par from the periodic cubic (consistent with the bounce points); # the periodic tspl B_val remains the numerator field in the integrands. - vpar = 1.0 - (lmda / bo) * B_extrap(θmod) + vpar = 1.0 - (lmda / bo) * B_vpar(θmod) if vpar <= 0 # Negative v_par near a bounce point: same fill rules as the Fortran bounce loop. diff --git a/src/KineticForces/Torque.jl b/src/KineticForces/Torque.jl index f9af86485..e4bbd3c35 100644 --- a/src/KineticForces/Torque.jl +++ b/src/KineticForces/Torque.jl @@ -94,10 +94,12 @@ function tpsi!(tpsi_var::Ref{ComplexF64}, psi::Float64, n::Int, l::Int, # Create periodic interpolant for poloidal quantities tspl = cubic_interp(xs, Series(hcat(B_vals, dBdpsi_vals, dBdtheta_vals, jac_vals, djdpsi_vals)); bc=PeriodicBC()) - # v_par and the bounce points use a separate endpoint-fit (non-periodic) cubic - # of B, like Fortran's vspl. The endpoint fit of 1−(λ/bo)B equals 1−(λ/bo) - # times the fit of B, so one B_extrap per surface serves every λ. - B_extrap = cubic_interp(xs, B_vals; bc=CubicFit()) + # v_par and the bounce points use a separate scalar cubic of B (Fortran's vspl). + # 1−(λ/bo)B is periodic on the closed θ interval, so the fit must be periodic too: + # a non-periodic endpoint fit is only C⁰ at the θ=0/1 seam and manufactures false + # near-seam extrema and root pairs there. The fit of 1−(λ/bo)B equals 1−(λ/bo) + # times the fit of B, so one B_vpar per surface serves every λ. + B_vpar = cubic_interp(xs, B_vals; bc=PeriodicBC()) bmax = maximum(B_vals) ibmax = argmax(B_vals) @@ -226,7 +228,7 @@ function tpsi!(tpsi_var::Ref{ComplexF64}, psi::Float64, n::Int, l::Int, method, op_wmats; chi1=intr.chi1, ro=intr.ro, mfac=intr.mfac, mpert=intr.mpert, theta_bmax=theta_bmax, - B_extrap=B_extrap, + B_vpar=B_vpar, smat=smat_f, tmat=tmat_f, xmat=xmat_f, ymat=ymat_f, zmat=zmat_f, energy_atol=atol_xlmda, energy_rtol=rtol_xlmda, @@ -401,7 +403,7 @@ function calculate_gar(psi, n, l, q, epsr, wdian, wdiat, welec, nuk, bo, bmax, bmin, n_s::Float64, T_s::Float64, mass, chrg, tspl, dbob_m_f, divx_m_f, divxfac, wdfac, method, op_wmats; chi1::Float64, ro::Float64, mfac::Vector{Int}, mpert::Int, - theta_bmax::Float64, B_extrap, + theta_bmax::Float64, B_vpar, smat=nothing, tmat=nothing, xmat=nothing, ymat=nothing, zmat=nothing, nlmda::Int=128, ntheta::Int=128, @@ -415,7 +417,7 @@ function calculate_gar(psi, n, l, q, epsr, wdian, wdiat, welec, nuk, bo, # Bounce-averaged quantities per pitch angle bounce = compute_bounce_data( psi, n, l, q, bo, bmax, bmin, theta_bmax, - tspl, B_extrap, mfac, chi1, ro, dbob_m_f, divx_m_f, divxfac, wdfac, + tspl, B_vpar, mfac, chi1, ro, dbob_m_f, divx_m_f, divxfac, wdfac, mass, chrg, T_s, method; nlmda, ntheta, smat, tmat, xmat, ymat, zmat) @@ -615,8 +617,9 @@ function _setup_surface_state( end tspl = cubic_interp(xs, Series(hcat(B_vals, dBdpsi_vals, dBdtheta_vals, jac_vals, djdpsi_vals)); bc=PeriodicBC()) - # Endpoint-fit (non-periodic) cubic of B for v_par and bounce points (Fortran vspl equivalent). - B_extrap = cubic_interp(xs, B_vals; bc=CubicFit()) + # Periodic cubic of B for v_par and bounce points (Fortran vspl equivalent); 1−(λ/bo)B + # is periodic on the closed θ interval, so a non-periodic fit would break at the seam. + B_vpar = cubic_interp(xs, B_vals; bc=PeriodicBC()) bmax = maximum(B_vals) ibmax = argmax(B_vals) @@ -685,7 +688,7 @@ function _setup_surface_state( return (; chrg, mass, - tspl, B_extrap, bmax, bmin, theta_bmax, + tspl, B_vpar, bmax, bmin, theta_bmax, q, n_s, T_s, welec, wdian, wdiat, wtran, wgyro, nuk, epsr, @@ -750,7 +753,7 @@ function kinetic_energy_matrices_for_euler_lagrange!( bounce = compute_bounce_data( psi, n, l, state.q, bo, state.bmax, state.bmin, state.theta_bmax, - state.tspl, state.B_extrap, mfac, chi1, ro, dbob_m_f, divx_m_f, 1.0, wdfac, + state.tspl, state.B_vpar, mfac, chi1, ro, dbob_m_f, divx_m_f, 1.0, wdfac, state.mass, state.chrg, state.T_s, "fwmm"; nlmda, ntheta, smat=smat_f, tmat=tmat_f, xmat=xmat_f, ymat=ymat_f, zmat=zmat_f)