Use FastInterpolations Splines with no adaptations

docs/make.jl

@@ -24,12 +24,11 @@ makedocs(;
         "Home" => "index.md",
         "Setup" => "set_up.md",
         "API Reference" => [
-            "Splines" => "splines.md",
             "Vacuum" => "vacuum.md",
             "Equilibrium" => "equilibrium.md",
             "Utilities" => "utilities.md",
             "Perturbed Equilibrium" => "perturbed_equilibrium.md"
-        ],
+        ]
     ],
     checkdocs=:exports
 )

src/Equilibrium/AnalyticEquilibrium.jl

@@ -208,10 +208,10 @@ function lar_run(equil_input::EquilibriumConfig, lar_input::LargeAspectRatioConf
     # Create separate interpolants for R and Z coordinates
     rz_in_xs = r_nodes
     rz_in_ys = collect(rzphi_y_nodes)
-    rz_in_R = cubic_interp((rz_in_xs, rz_in_ys), rzphi_fs_nodes[:, :, 1];
-        bc=(CubicFit(), Spl.PeriodicBC()), extrap=(:extension, :wrap))
-    rz_in_Z = cubic_interp((rz_in_xs, rz_in_ys), rzphi_fs_nodes[:, :, 2];
-        bc=(CubicFit(), Spl.PeriodicBC()), extrap=(:extension, :wrap))
+    rz_in_R = cubic_interp((rz_in_xs, rz_in_ys), rzphi_fs_nodes[:, :, 1]; search=LinearBinary(),
+        bc=(CubicFit(), PeriodicBC()), extrap=(:extension, :wrap))
+    rz_in_Z = cubic_interp((rz_in_xs, rz_in_ys), rzphi_fs_nodes[:, :, 2]; search=LinearBinary(),
+        bc=(CubicFit(), PeriodicBC()), extrap=(:extension, :wrap))
 
     # LAR equilibrium has midplane symmetry, so magnetic axis is at Z = 0
     return InverseRunInput(equil_input, sq_in, rz_in_xs, rz_in_ys, rz_in_R, rz_in_Z, lar_r0, 0.0, psio)
@@ -288,7 +288,7 @@ function sol_run(equil_inputs::EquilibriumConfig, sol_inputs::SolovevConfig)
     end
     psi_in_xs = r
     psi_in_ys = z
-    psi_in = cubic_interp((psi_in_xs, psi_in_ys), psifs;
+    psi_in = cubic_interp((psi_in_xs, psi_in_ys), psifs; search=LinearBinary(),
         bc=(CubicFit(), CubicFit()), extrap=(:extension, :extension))
 
     # Print out equilibrium info

src/Equilibrium/DirectEquilibrium.jl

@@ -191,7 +191,7 @@ function direct_position!(raw_profile::DirectRunInput)
     # Access nodal values from psi_in interpolant: partials[1,:,:] = function values
     new_psi_fs = raw_profile.psi_in.nodal_derivs.partials[1, :, :] .* raw_profile.psio / bfield.psi
     # Because DirectRunInput is a mutable struct, we can update the spline here
-    raw_profile.psi_in = cubic_interp((x_coords, y_coords), new_psi_fs;
+    raw_profile.psi_in = cubic_interp((x_coords, y_coords), new_psi_fs; search=LinearBinary(),
         bc=(CubicFit(), CubicFit()), extrap=(:extension, :extension))
 
     # Helper function for robust Newton-Raphson search with restarts
@@ -467,7 +467,7 @@ function equilibrium_solver(raw_profile::DirectRunInput)
         ff_fs_nodes[end, :] .= ff_fs_nodes[1, :]
 
         # Create series interpolant for all columns
-        ff_interp = cubic_interp(ff_x_nodes, ff_fs_nodes; bc=Spl.PeriodicBC())
+        ff_interp = cubic_interp(ff_x_nodes, ff_fs_nodes; bc=PeriodicBC())
         ff_deriv = deriv1(ff_interp)
 
         # Interpolate `ff` onto the uniform `theta` grid for `rzphi`
@@ -524,14 +524,14 @@ function equilibrium_solver(raw_profile::DirectRunInput)
     # theta_nodes includes both 0 and 1 (closed periodic grid).
     rzphi_xs = psi_nodes
     rzphi_ys = collect(theta_nodes)
-    rzphi_rsquared = cubic_interp((rzphi_xs, rzphi_ys), rzphi_nodes[:, :, 1];
-        bc=(CubicFit(), Spl.PeriodicBC()), extrap=(:extension, :wrap))
-    rzphi_offset = cubic_interp((rzphi_xs, rzphi_ys), rzphi_nodes[:, :, 2];
-        bc=(CubicFit(), Spl.PeriodicBC()), extrap=(:extension, :wrap))
-    rzphi_nu = cubic_interp((rzphi_xs, rzphi_ys), rzphi_nodes[:, :, 3];
-        bc=(CubicFit(), Spl.PeriodicBC()), extrap=(:extension, :wrap))
-    rzphi_jac = cubic_interp((rzphi_xs, rzphi_ys), rzphi_nodes[:, :, 4];
-        bc=(CubicFit(), Spl.PeriodicBC()), extrap=(:extension, :wrap))
+    rzphi_rsquared = cubic_interp((rzphi_xs, rzphi_ys), rzphi_nodes[:, :, 1]; search=LinearBinary(),
+        bc=(CubicFit(), PeriodicBC()), extrap=(:extension, :wrap))
+    rzphi_offset = cubic_interp((rzphi_xs, rzphi_ys), rzphi_nodes[:, :, 2]; search=LinearBinary(),
+        bc=(CubicFit(), PeriodicBC()), extrap=(:extension, :wrap))
+    rzphi_nu = cubic_interp((rzphi_xs, rzphi_ys), rzphi_nodes[:, :, 3]; search=LinearBinary(),
+        bc=(CubicFit(), PeriodicBC()), extrap=(:extension, :wrap))
+    rzphi_jac = cubic_interp((rzphi_xs, rzphi_ys), rzphi_nodes[:, :, 4]; search=LinearBinary(),
+        bc=(CubicFit(), PeriodicBC()), extrap=(:extension, :wrap))
 
     # Calculate physics quantities (B-field, metric components, etc.) in 2D spline `eqfun`
     # for use in stability and transport codes
@@ -596,12 +596,12 @@ function equilibrium_solver(raw_profile::DirectRunInput)
         end
     end
     # Create 2D interpolants for physics quantities (eqfun)
-    eqfun_B = cubic_interp((rzphi_xs, rzphi_ys), eqfun_fs_nodes[:, :, 1];
-        bc=(CubicFit(), Spl.PeriodicBC()), extrap=(:extension, :wrap))
-    eqfun_metric1 = cubic_interp((rzphi_xs, rzphi_ys), eqfun_fs_nodes[:, :, 2];
-        bc=(CubicFit(), Spl.PeriodicBC()), extrap=(:extension, :wrap))
-    eqfun_metric2 = cubic_interp((rzphi_xs, rzphi_ys), eqfun_fs_nodes[:, :, 3];
-        bc=(CubicFit(), Spl.PeriodicBC()), extrap=(:extension, :wrap))
+    eqfun_B = cubic_interp((rzphi_xs, rzphi_ys), eqfun_fs_nodes[:, :, 1]; search=LinearBinary(),
+        bc=(CubicFit(), PeriodicBC()), extrap=(:extension, :wrap))
+    eqfun_metric1 = cubic_interp((rzphi_xs, rzphi_ys), eqfun_fs_nodes[:, :, 2]; search=LinearBinary(),
+        bc=(CubicFit(), PeriodicBC()), extrap=(:extension, :wrap))
+    eqfun_metric2 = cubic_interp((rzphi_xs, rzphi_ys), eqfun_fs_nodes[:, :, 3]; search=LinearBinary(),
+        bc=(CubicFit(), PeriodicBC()), extrap=(:extension, :wrap))
 
     return PlasmaEquilibrium(raw_profile.config, EquilibriumParameters(), profiles,
         rzphi_xs, rzphi_ys,

src/Equilibrium/Equilibrium.jl

@@ -1,12 +1,11 @@
 module Equilibrium
 
 # --- Module-level Dependencies ---
-import ..Spl
 
 using Printf, OrdinaryDiffEq, DiffEqCallbacks, LinearAlgebra, HDF5
 using TOML
 import FastInterpolations
-using FastInterpolations: cubic_interp, deriv1, deriv2, deriv3, LinearBinary, CubicFit
+using FastInterpolations: cubic_interp, deriv1, deriv2, deriv3, LinearBinary, CubicFit, PeriodicBC
 import StaticArrays: @MMatrix, SVector
 
 # --- Internal Module Structure ---
@@ -107,11 +106,12 @@ function equilibrium_separatrix_find!(pe::PlasmaEquilibrium)
 
     for iside in 1:2
         it = 0
+        hint2d = (Ref(1), Ref(1))  # Shared 2D hint for Newton iteration
         while true
             it += 1
             # Evaluate offset and its derivative
-            offset = pe.rzphi_offset((psi_edge, theta))
-            offset_y = pe.rzphi_offset((psi_edge, theta); deriv=(0,1))
+            offset = pe.rzphi_offset((psi_edge, theta); hint=hint2d)
+            offset_y = pe.rzphi_offset((psi_edge, theta); deriv=(0, 1), hint=hint2d)
 
             eta = theta + offset - eta0
             eta_theta = 1 + offset_y
@@ -143,15 +143,16 @@ function equilibrium_separatrix_find!(pe::PlasmaEquilibrium)
         z = 0.0
         max_iter = 1000
         iter = 0
+        hint2d = (Ref(1), Ref(1))  # Shared 2D hint for Newton iteration
         while iter < max_iter
             iter += 1
             # Evaluate rcoord and offset with derivatives
-            r2 = pe.rzphi_rsquared((psi_edge, theta))
-            r2y = pe.rzphi_rsquared((psi_edge, theta); deriv=(0,1))
-            r2yy = pe.rzphi_rsquared((psi_edge, theta); deriv=(0,2))
-            eta = pe.rzphi_offset((psi_edge, theta))
-            eta1 = pe.rzphi_offset((psi_edge, theta); deriv=(0,1))
-            eta2 = pe.rzphi_offset((psi_edge, theta); deriv=(0,2))
+            r2 = pe.rzphi_rsquared((psi_edge, theta); hint=hint2d)
+            r2y = pe.rzphi_rsquared((psi_edge, theta); deriv=(0, 1), hint=hint2d)
+            r2yy = pe.rzphi_rsquared((psi_edge, theta); deriv=(0, 2), hint=hint2d)
+            eta = pe.rzphi_offset((psi_edge, theta); hint=hint2d)
+            eta1 = pe.rzphi_offset((psi_edge, theta); deriv=(0, 1), hint=hint2d)
+            eta2 = pe.rzphi_offset((psi_edge, theta); deriv=(0, 2), hint=hint2d)
 
             rfac = sqrt(r2)
             rfac1 = r2y / (2 * rfac)
@@ -441,17 +442,18 @@ function equilibrium_gse!(equil::PlasmaEquilibrium)
         end
     end
     # Create flux interpolants for Grad-Shafranov diagnostics
-    flux1 = cubic_interp((equil.rzphi_xs, equil.rzphi_ys), flux_fs[:, :, 1];
-        bc=(CubicFit(), Spl.PeriodicBC()), extrap=(:extension, :wrap))
-    flux2 = cubic_interp((equil.rzphi_xs, equil.rzphi_ys), flux_fs[:, :, 2];
-        bc=(CubicFit(), Spl.PeriodicBC()), extrap=(:extension, :wrap))
+    flux1 = cubic_interp((equil.rzphi_xs, equil.rzphi_ys), flux_fs[:, :, 1]; search=LinearBinary(),
+        bc=(CubicFit(), PeriodicBC()), extrap=(:extension, :wrap))
+    flux2 = cubic_interp((equil.rzphi_xs, equil.rzphi_ys), flux_fs[:, :, 2]; search=LinearBinary(),
+        bc=(CubicFit(), PeriodicBC()), extrap=(:extension, :wrap))
     # Compute flux derivatives at all grid points for diagnostics
+    hint2d = (Ref(1), Ref(1))  # Shared 2D hint for hot loop optimization
     for ipsi in 0:mpsi
         for itheta in 0:mtheta
-            flux_fsx[ipsi+1, itheta+1, 1] = flux1((equil.rzphi_xs[ipsi+1], equil.rzphi_ys[itheta+1]); deriv=(1,0))
-            flux_fsx[ipsi+1, itheta+1, 2] = flux2((equil.rzphi_xs[ipsi+1], equil.rzphi_ys[itheta+1]); deriv=(1,0))
-            flux_fsy[ipsi+1, itheta+1, 1] = flux1((equil.rzphi_xs[ipsi+1], equil.rzphi_ys[itheta+1]); deriv=(0,1))
-            flux_fsy[ipsi+1, itheta+1, 2] = flux2((equil.rzphi_xs[ipsi+1], equil.rzphi_ys[itheta+1]); deriv=(0,1))
+            flux_fsx[ipsi+1, itheta+1, 1] = flux1((equil.rzphi_xs[ipsi+1], equil.rzphi_ys[itheta+1]); deriv=(1, 0), hint=hint2d)
+            flux_fsx[ipsi+1, itheta+1, 2] = flux2((equil.rzphi_xs[ipsi+1], equil.rzphi_ys[itheta+1]); deriv=(1, 0), hint=hint2d)
+            flux_fsy[ipsi+1, itheta+1, 1] = flux1((equil.rzphi_xs[ipsi+1], equil.rzphi_ys[itheta+1]); deriv=(0, 1), hint=hint2d)
+            flux_fsy[ipsi+1, itheta+1, 2] = flux2((equil.rzphi_xs[ipsi+1], equil.rzphi_ys[itheta+1]); deriv=(0, 1), hint=hint2d)
         end
     end
 
@@ -492,8 +494,9 @@ function equilibrium_gse!(equil::PlasmaEquilibrium)
         fs_matrix[:, 1] = flux_fsx[ipsi, :, 1]
         fs_matrix[:, 2] = source[ipsi, :]
 
-        # Compute total integral using exact spline integration (only final value needed)
-        term[ipsi, :] .= Spl.total_integral(equil.rzphi_ys, fs_matrix; bc=Spl.PeriodicBC())
+        # Compute total integral using FastInterpolations native integration
+        itp = cubic_interp(equil.rzphi_ys, fs_matrix; bc=PeriodicBC())
+        term[ipsi, :] .= FastInterpolations.integrate(itp)
     end
 
     totali = sum(term; dims=2)