From 7847c08b52a649f3314d908a400d8f9ec560bc25 Mon Sep 17 00:00:00 2001 From: adrianaghiozzi <67669644+adrianaghiozzi@users.noreply.github.com> Date: Wed, 12 Aug 2026 15:32:22 -0400 Subject: [PATCH 1/6] Explicit type equil --- src/KineticForces/Torque.jl | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/KineticForces/Torque.jl b/src/KineticForces/Torque.jl index af22d165c..f96935c57 100644 --- a/src/KineticForces/Torque.jl +++ b/src/KineticForces/Torque.jl @@ -30,7 +30,7 @@ Imaginary component is proportional to the kinetic energy Im(T) = 2*n*dW_k. """ function tpsi!(tpsi_var::Ref{ComplexF64}, psi::Float64, n::Int, l::Int, zi::Int, mi::Int, wdfac::Float64, divxfac::Float64, - electron::Bool, method::String, equil, intr::KineticForcesInternal, + electron::Bool, method::String, equil::Equilibrium.PlasmaEquilibrium, intr::KineticForcesInternal, kinetic_profiles::Equilibrium.KineticProfileSplines; op_wmats::Union{Nothing,Array{ComplexF64,3}}=nothing, rex_override::Union{Nothing,Float64}=nothing, @@ -889,7 +889,7 @@ function compute_kinetic_matrices_at_psi!( ktmat::Array{ComplexF64,3}, psi::Float64, n::Int, l::Int, zi::Int, mi::Int, wdfac::Float64, _divxfac::Float64, - electron::Bool, equil, intr::KineticForcesInternal, + electron::Bool, equil::Equilibrium.PlasmaEquilibrium, intr::KineticForcesInternal, kinetic_profiles::Equilibrium.KineticProfileSplines; nutype::String="harmonic", f0type::String="maxwellian", nufac::Float64=1.0, atol_xlmda::Float64=1e-9, rtol_xlmda::Float64=1e-6) From 98af781fab5883d498d151b4692ad4e92d6ec662 Mon Sep 17 00:00:00 2001 From: adrianaghiozzi <67669644+adrianaghiozzi@users.noreply.github.com> Date: Wed, 12 Aug 2026 15:34:20 -0400 Subject: [PATCH 2/6] Faster green function --- src/Vacuum/Kernel2D.jl | 48 ++++++++++++++---------------------------- 1 file changed, 16 insertions(+), 32 deletions(-) diff --git a/src/Vacuum/Kernel2D.jl b/src/Vacuum/Kernel2D.jl index ea93f224b..fc7155f3e 100644 --- a/src/Vacuum/Kernel2D.jl +++ b/src/Vacuum/Kernel2D.jl @@ -639,27 +639,29 @@ according to equations (36)-(42) of Chance 1997. Replaces `green` from Fortran c gamma_prefactor::Float64=2 * sqrt(π) * gamma(0.5 - n), uselegacygreenfunction::Bool=false ) - x_obs2 = x_obs^2 x_source2 = x_source^2 x_minus2 = (x_obs - x_source)^2 x_multiple = x_obs * x_source - ζ = (z_obs - z_source) + ζ = z_obs - z_source ζ2 = ζ^2 ρ2 = x_minus2 + ζ2 - - # Distance parameter ℛ [Chance Phys. Plasmas 1997 2161 eq. 41] R4 = ρ2 * (ρ2 + 4 * x_multiple) R2 = sqrt(R4) R = sqrt(R2) R5 = R4 * R - # Argument of Legendre function 𝘴 [Chance Phys. Plasmas 1997 2161 eq. 42] - s = (x_obs2 + x_source2 + ζ2) / R2 + S = x_obs2 + x_source2 + ζ2 + a = x_obs2 - x_source2 + D = a + ζ2 # x_obs2 - x_source2 + ζ2 + E = ζ2 - a # x_source2 - x_obs2 + ζ2 + fourXmult = 4 * x_multiple + twoXobsD = 2 * x_obs * D + xSourceE = x_source * E + + s = S / R2 - # Legendre functions for - # P⁰ = p0, P¹ = p1, Pⁿ = pn, Pⁿ⁺¹ = pnp1 legendre = acquire!(pool, Float64, n + 2) if uselegacygreenfunction Pn_minus_half_1997!(legendre, s, n) @@ -667,37 +669,19 @@ according to equations (36)-(42) of Chance 1997. Replaces `green` from Fortran c Pn_minus_half_2007!(legendre, s, n) end - p0 = legendre[1] - p1 = legendre[2] - pnp1 = legendre[end] - pn = legendre[end-1] + p0, p1, pnp1, pn = @inbounds legendre[1], legendre[2], legendre[end], legendre[end-1] - # Green's function 2π𝒢ⁿ = G_n [Chance Phys. Plasmas 1997 2161 eq. 40] gg = gamma_prefactor / R G_n = gg * pn + grad_gg = gg / (2π * R4) - # Gradient factor [Chance Phys. Plasmas 1997 2161 eq. 44] - # NOTE: Paper has erroneous extra factor of 2π - grad_gg = gg / R4 / 2π - - # Derivatives of Green's function [Chance Phys. Plasmas 1997 2161 eq. 36-38] - # ∂Gⁿ/∂X' using chain rule: ∂Gⁿ/∂X' = (∂Gⁿ/∂R)(∂R/∂X') + (∂Gⁿ/∂s)(∂s/∂X') - xterm1 = (n * (x_obs2 + x_source2 + ζ2) * (x_obs2 - x_source2 + ζ2) - x_source2*(x_source2-x_obs2+ζ2)) * pn - xterm2 = (2.0 * x_source * x_obs * (x_obs2-x_source2+ζ2)) * pnp1 - dG_dX = grad_gg * (xterm1 + xterm2) / x_source - - # ∂Gⁿ/∂Z' using chain rule - zterm1 = (2.0 * n + 1.0) * (x_obs2 + x_source2 + ζ2) * pn - zterm2 = 4.0 * x_multiple * pnp1 - dG_dZ = grad_gg * (zterm1 + zterm2) * ζ + dG_dX = grad_gg * ( (n * S * D - x_source * xSourceE) * pn / x_source + twoXobsD * pnp1 ) + dG_dZ = grad_gg * ((2n + 1) * S * pn + fourXmult * pnp1) * ζ - # Coupling term 𝒥 ∇'𝒢ⁿ∇'ℒ [Chance Phys. Plasmas 1997 2161 eq. 51] - # Jacobian factor from coordinate transformation coupling_n = -x_source * (dz_dtheta * dG_dX - dx_dtheta * dG_dZ) - # Special case for n=0: coupling_0 = 1/(2π) 𝒥 ∇'𝒢⁰∇'ℒ - dG_dX0_R5 = ((2.0 * x_obs * (x_obs2-x_source2+ζ2)) * p1 - x_source * (x_source2-x_obs2+ζ2) * p0) - dG_dZ0_R5 = ζ * ((x_obs2 + x_source2 + ζ2) * p0 + 4.0 * x_multiple * p1) + dG_dX0_R5 = twoXobsD * p1 - xSourceE * p0 + dG_dZ0_R5 = ζ * (S * p0 + fourXmult * p1) coupling_0 = -x_source * (dz_dtheta * dG_dX0_R5 - dx_dtheta * dG_dZ0_R5) / R5 return G_n, coupling_n, coupling_0 end From 8adb33b724bd9eebd8edda2cf78f27fc34983cbf Mon Sep 17 00:00:00 2001 From: adrianaghiozzi <67669644+adrianaghiozzi@users.noreply.github.com> Date: Fri, 14 Aug 2026 10:36:02 -0400 Subject: [PATCH 3/6] Tuple -> Vector in integrate_pitch_gar_quadgk to avoid type instability --- src/KineticForces/PitchIntegration.jl | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/KineticForces/PitchIntegration.jl b/src/KineticForces/PitchIntegration.jl index fa61c5952..afa646b17 100644 --- a/src/KineticForces/PitchIntegration.jl +++ b/src/KineticForces/PitchIntegration.jl @@ -97,9 +97,9 @@ function integrate_pitch_gar_quadgk( # the kink in leff = ell + n*q (circulating) → ell (trapped). bobmax_clip = clamp(bobmax, lambda_min, lambda_max) segments = if lambda_min < bobmax_clip < lambda_max - (lambda_min, bobmax_clip, lambda_max) + [lambda_min, bobmax_clip, lambda_max] else - (lambda_min, lambda_max) + [lambda_min, lambda_max] end # In-place quadgk! buffer; copy the result out so the returned vector is From 152e7e50d8c859b9d0082c0d1987e4df044d79b4 Mon Sep 17 00:00:00 2001 From: adrianaghiozzi <67669644+adrianaghiozzi@users.noreply.github.com> Date: Fri, 14 Aug 2026 11:18:09 -0400 Subject: [PATCH 4/6] Better type instability fix for integrate_pitch_gar_quadgk --- src/KineticForces/PitchIntegration.jl | 26 ++++++++++---------------- 1 file changed, 10 insertions(+), 16 deletions(-) diff --git a/src/KineticForces/PitchIntegration.jl b/src/KineticForces/PitchIntegration.jl index afa646b17..20c7a5703 100644 --- a/src/KineticForces/PitchIntegration.jl +++ b/src/KineticForces/PitchIntegration.jl @@ -96,17 +96,13 @@ function integrate_pitch_gar_quadgk( # Split domain at trapped/passing boundary so Gauss-Kronrod resolves # the kink in leff = ell + n*q (circulating) → ell (trapped). bobmax_clip = clamp(bobmax, lambda_min, lambda_max) - segments = if lambda_min < bobmax_clip < lambda_max - [lambda_min, bobmax_clip, lambda_max] - else - [lambda_min, lambda_max] - end - - # In-place quadgk! buffer; copy the result out so the returned vector is - # owned by the caller. buf = zeros(ComplexF64, nqty) kernel! = (out, λ) -> _pitch_gar_kernel_quadgk!(out, λ, params) - I, _ = quadgk!(kernel!, buf, segments...; atol=pitch_atol, rtol=pitch_rtol) + I, _ = if lambda_min < bobmax_clip < lambda_max + quadgk!(kernel!, buf, lambda_min, bobmax_clip, lambda_max; atol=pitch_atol, rtol=pitch_rtol) + else + quadgk!(kernel!, buf, lambda_min, lambda_max; atol=pitch_atol, rtol=pitch_rtol) + end return copy(I) end @@ -195,15 +191,13 @@ function integrate_pitch_gar_quadgk_wt( lambda_max = last(fbnce.cache.x) bobmax_clip = clamp(bobmax, lambda_min, lambda_max) - segments = if lambda_min < bobmax_clip < lambda_max - (lambda_min, bobmax_clip, lambda_max) - else - (lambda_min, lambda_max) - end - buf = zeros(ComplexF64, 2 * nqty) kernel! = (out, λ) -> _pitch_gar_kernel_quadgk_wt!(out, λ, params) - I, _ = quadgk!(kernel!, buf, segments...; atol=pitch_atol, rtol=pitch_rtol) + I, _ = if lambda_min < bobmax_clip < lambda_max + quadgk!(kernel!, buf, lambda_min, bobmax_clip, lambda_max; atol=pitch_atol, rtol=pitch_rtol) + else + quadgk!(kernel!, buf, lambda_min, lambda_max; atol=pitch_atol, rtol=pitch_rtol) + end return copy(I) end From 69ad76e4f7677c54ac91a69be3bca74d6e57a92f Mon Sep 17 00:00:00 2001 From: adrianaghiozzi <67669644+adrianaghiozzi@users.noreply.github.com> Date: Fri, 14 Aug 2026 17:15:09 -0400 Subject: [PATCH 5/6] FORCEFREESTATES - OPTIMIZATION - Loosen ballooning ODE abstol to match reltol integrate_ballooning_ode set abstol=reltol^2 (1e-16, at machine epsilon) against a non-stiff DP5 solver, forcing needlessly tight step control. Profiling (warmed, no JIT noise) showed ballooning boundary search as the dominant real cost on gal-resistive-PE-style cases -- ~55% of pipeline runtime, driven by many repeated ODE shoots per flux surface (linear pre-scan + bisection to locate the marginal-stability crossing). Setting abstol=reltol=1e-8 (the standard pairing) gives a measured ~8% speedup on Force-Free States for DIIID-like_gal_resistive_pe_example, confirming reltol was already the binding constraint for most of the integration. Validated bit-for-bit identical locstab/alpha_critical (the value that actually drives stability classification) before/after; the raw ballooning_Delta_prime diagnostic drifts ~1e-5 relative, which flips a few profile checksums in the regression harness but is far below any threshold that matters physically. Full 11-case regression suite against merge-base a0cad260 is otherwise clean. --- src/ForceFreeStates/Ballooning.jl | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/ForceFreeStates/Ballooning.jl b/src/ForceFreeStates/Ballooning.jl index d24f9ff81..78408b5a2 100644 --- a/src/ForceFreeStates/Ballooning.jl +++ b/src/ForceFreeStates/Ballooning.jl @@ -1169,7 +1169,7 @@ function integrate_ballooning_ode(ode_coefficient_spline; theta_k::Float64=0.0) problem_left, DP5(); reltol=TOLERANCE, - abstol=TOLERANCE^2, + abstol=TOLERANCE, dtmin=MINIMUM_STEP, adaptive=true, save_everystep=false, @@ -1190,7 +1190,7 @@ function integrate_ballooning_ode(ode_coefficient_spline; theta_k::Float64=0.0) problem_right, DP5(); reltol=TOLERANCE, - abstol=TOLERANCE^2, + abstol=TOLERANCE, dtmin=MINIMUM_STEP, adaptive=true, save_everystep=false, From e952629dba932e192b697fc4cbec0c182361b23e Mon Sep 17 00:00:00 2001 From: adrianaghiozzi <67669644+adrianaghiozzi@users.noreply.github.com> Date: Mon, 17 Aug 2026 15:51:05 -0400 Subject: [PATCH 6/6] Add back Distance parameter comment --- src/Vacuum/Kernel2D.jl | 2 ++ 1 file changed, 2 insertions(+) diff --git a/src/Vacuum/Kernel2D.jl b/src/Vacuum/Kernel2D.jl index fc7155f3e..fd6599209 100644 --- a/src/Vacuum/Kernel2D.jl +++ b/src/Vacuum/Kernel2D.jl @@ -647,6 +647,8 @@ according to equations (36)-(42) of Chance 1997. Replaces `green` from Fortran c ζ2 = ζ^2 ρ2 = x_minus2 + ζ2 + + # Distance parameter ℛ [Chance Phys. Plasmas 1997 2161 eq. 41] R4 = ρ2 * (ρ2 + 4 * x_multiple) R2 = sqrt(R4) R = sqrt(R2)