Switch to a standardized spline package #70
Description
Summarizing discussion in #68:
@logan-nc suggested we spend out time using an existing Julia spline package rather than converting the custom splines from fortran to Julia and doing lots of modifications to make them more readible and thread stable.
@bclyons12 listed the most well-known ones to consider:
@logan-nc identified a winner:
Gemini and ChatGPT both recommend Interpolations.jl.
Gemini says
Yes, you can use Interpolations.jl in parallelized code, such as with Julia's built-in multithreading and distributed computing capabilities. An Interpolations.jl object is read-only after construction, making it safe to use across multiple threads for interpolation lookups.
ChatGPT warns agains Direckx.jl
Avoid Dierckx.jl for threaded evaluation of derivatives. There’s an open issue showing wildly incorrect results when calling Dierckx.derivative from multiple threads; single-thread is fine, but multithreaded derivatives are not thread-safe in practice.
while recommending Interpolations.jl as
Pure Julia; very fast construction & throughput on both B-spline and gridded methods; handles multi-D regular and irregularly spaced axes; well-documented boundary conditions/extrapolation; works with modern Julia versions. The docs include benchmark “shootouts” vs Dierckx and others where Interpolations wins across cases. Maintenance stance is “maintained, semver, fixes prioritized,” which signals low surprise from version to version.
The main complaint for Interpolations.jl seems to be the assumption of regular grids, but there is an option to scale the grid that I think would work for us.... except that I don't think we even need it!
- Our irregular equilibrium psi grids are all just known functional packings of the points so we can spline in regular variable
xwherex=psin^2and then easily just eval likeqspl(psin^2)or even add explicitx=psin^2; qspl(x)type stuff in all the functions that happen at a psin value - GPEC splines the solutions on the dynamically determined force-free solution grid (
ugrid) but only ever does calculations at nodes - so that is just wasted of overhead really. - Theta grids are uniform.
Conclusion
Someone needs to spend a day attempting to place Interpolations.jl splines in for the old fortran ones we currently have wrappers for.
My recommendation would be to just replace something easy like sq to see if it works and/or impacts speed. Then move on to others if it looks promising.