Skip to content

KineticForces - BUGFIX! - Fit the parallel-velocity spline periodically #390

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Open
logan-nc wants to merge 1 commit into
base: develop
Choose a base branch
from
+33 −30
Open

KineticForces - BUGFIX! - Fit the parallel-velocity spline periodically #390

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
1 commit
Select commit
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
38 changes: 19 additions & 19 deletions src/KineticForces/BounceAveraging.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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,
Expand All @@ -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]
Expand All @@ -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,
Expand Down Expand Up @@ -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)
Expand Down Expand Up @@ -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))


"""
Expand All @@ -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)
Expand All @@ -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
Expand All @@ -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
Expand All @@ -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]
Expand Down Expand Up @@ -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,
Expand Down Expand Up @@ -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.
Expand Down
25 changes: 14 additions & 11 deletions src/KineticForces/Torque.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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)
Expand Down Expand Up @@ -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,
Expand Down Expand Up @@ -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,
Expand All @@ -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)

Expand Down Expand Up @@ -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)
Expand Down Expand Up @@ -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,
Expand Down Expand Up @@ -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)

Expand Down
Loading