FastInterpolations bump and fix

Project.toml

@@ -40,7 +40,7 @@ DiffEqCallbacks = "4.9.0"
 Documenter = "1.14.1"
 FFTW = "1.9.0"
 FastGaussQuadrature = "1.1.0"
-FastInterpolations = "0.4"
+FastInterpolations = "0.4.10"
 HDF5 = "0.17.2"
 IMASdd = "8"
 JLD2 = "0.6.3"

src/Analysis/PerturbedEquilibriumModes.jl

@@ -71,7 +71,7 @@ function modes_to_theta(h5_file::String, variable::String;
             nu_vals  = read(f, "splines/rzphi/nu")
             nu_spline = cubic_interp(
                 (rzphi_xs, rzphi_ys), nu_vals;
-                bc=(CubicFit(), PeriodicBC(; check=false)),
+                bc=(CubicFit(), PeriodicBC()),
                 extrap=(ExtendExtrap(), WrapExtrap())
             )
 

src/Equilibrium/AnalyticEquilibrium.jl

@@ -217,7 +217,9 @@ function lar_run(equil_input::EquilibriumConfig, lar_input::LargeAspectRatioConf
     rz_in_xs = r_nodes
     rz_in_ys = collect(rzphi_y_nodes)
 
-    itp_2d_opts = (bc=(CubicFit(), PeriodicBC(; check=false)), extrap=(ExtendExtrap(), WrapExtrap()))
+    @views rzphi_fs_nodes[:, end, :] .= rzphi_fs_nodes[:, 1, :]
+
+    itp_2d_opts = (bc=(CubicFit(), PeriodicBC()), extrap=(ExtendExtrap(), WrapExtrap()))
     rz_in_R = cubic_interp((rz_in_xs, rz_in_ys), rzphi_fs_nodes[:, :, 1]; itp_2d_opts...)
     rz_in_Z = cubic_interp((rz_in_xs, rz_in_ys), rzphi_fs_nodes[:, :, 2]; itp_2d_opts...)
 

src/Equilibrium/DirectEquilibrium.jl

@@ -570,7 +570,7 @@ robustness.
 
         ff_fs_nodes[end, :] .= ff_fs_nodes[1, :]  # enforce periodic endpoint
 
-        ff_interp = cubic_interp(ff_x_nodes, Series(ff_fs_nodes); bc=PeriodicBC(; check=false))
+        ff_interp = cubic_interp(ff_x_nodes, Series(ff_fs_nodes); bc=PeriodicBC())
         ff_deriv = deriv1(ff_interp)
 
         # Resample ff onto uniform theta grid
@@ -632,7 +632,10 @@ robustness.
     rzphi_ys = collect(theta_nodes)
 
     grid2d = (rzphi_xs, theta_nodes)
-    opts2d = (bc=(CubicFit(), PeriodicBC(; check=false)), extrap=(ExtendExtrap(), WrapExtrap()))
+    opts2d = (bc=(CubicFit(), PeriodicBC()), extrap=(ExtendExtrap(), WrapExtrap()))
+
+    # Snap periodic endpoint: ff_interp evaluation at theta_nodes[end] may drift by machine eps from theta_nodes[1]
+    @views rzphi_nodes[:, end, :] .= rzphi_nodes[:, 1, :]
 
     rzphi_rsquared = cubic_interp(grid2d, rzphi_nodes[:, :, 1]; opts2d...)
     rzphi_offset = cubic_interp(grid2d, rzphi_nodes[:, :, 2]; opts2d...)
@@ -697,6 +700,8 @@ robustness.
         end
     end
 
+    @views eqfun_fs_nodes[:, end, :] .= eqfun_fs_nodes[:, 1, :]
+
     eqfun_B = cubic_interp(grid2d, eqfun_fs_nodes[:, :, 1]; opts2d...)
     eqfun_metric1 = cubic_interp(grid2d, eqfun_fs_nodes[:, :, 2]; opts2d...)
     eqfun_metric2 = cubic_interp(grid2d, eqfun_fs_nodes[:, :, 3]; opts2d...)

src/Equilibrium/DirectEquilibriumByInversion.jl

@@ -150,7 +150,7 @@ function resample_contour_to_theta_grid!(
     η_ext[n_ext] = η_ext[1] + 2π
     logρ_ext[n_ext] = logρ_ext[1]
 
-    logρ_spl = cubic_interp(η_ext, logρ_ext; bc=PeriodicBC(; check=false))
+    logρ_spl = cubic_interp(η_ext, logρ_ext; bc=PeriodicBC())
 
     # theta_grid is in turns [0, 1]; sample at 2π*θ radians and convert polar → Cartesian.
     # Monotonically increasing η → shared hint gives O(1) lookups per step.
@@ -658,7 +658,10 @@ function equilibrium_solver_by_inversion(
     # Build InverseRunInput — same type consumed by equilibrium_solver(::InverseRunInput)
     rz_in_xs = psi_nodes
     rz_in_ys = theta_grid
-    itp_opts2d = (bc=(CubicFit(), PeriodicBC(; check=false)), extrap=(ExtendExtrap(), WrapExtrap()))
+    itp_opts2d = (bc=(CubicFit(), PeriodicBC()), extrap=(ExtendExtrap(), WrapExtrap()))
+
+    @views R_table[:, end] .= R_table[:, 1]
+    @views Z_table[:, end] .= Z_table[:, 1]
 
     rz_in_R = cubic_interp((rz_in_xs, rz_in_ys), R_table; itp_opts2d...)
     rz_in_Z = cubic_interp((rz_in_xs, rz_in_ys), Z_table; itp_opts2d...)

src/Equilibrium/Equilibrium.jl

@@ -457,7 +457,8 @@ function equilibrium_gse!(equil::PlasmaEquilibrium)
         end
     end
     # Create flux interpolants for Grad-Shafranov diagnostics
-    flux_opts = (bc=(CubicFit(), PeriodicBC(; check=false)), extrap=(ExtendExtrap(), WrapExtrap()))
+    @views flux_fs[:, end, :] .= flux_fs[:, 1, :]
+    flux_opts = (bc=(CubicFit(), PeriodicBC()), extrap=(ExtendExtrap(), WrapExtrap()))
     flux1 = cubic_interp((equil.rzphi_xs, equil.rzphi_ys), flux_fs[:, :, 1]; flux_opts...)
     flux2 = cubic_interp((equil.rzphi_xs, equil.rzphi_ys), flux_fs[:, :, 2]; flux_opts...)
 

src/Equilibrium/InverseEquilibrium.jl

@@ -119,7 +119,7 @@ function equilibrium_solver(input::InverseRunInput)
     @views r2[:, end] .= r2[:, 1]
     @views deta[:, end] .= deta[:, 1]
 
-    itp_opts2d = (bc=(CubicFit(), PeriodicBC(; check=false)), extrap=(ExtendExtrap(), WrapExtrap()))
+    itp_opts2d = (bc=(CubicFit(), PeriodicBC()), extrap=(ExtendExtrap(), WrapExtrap()))
 
     # Create 2D interpolants for r² and dη
     rz_rsq = cubic_interp((rz_in_xs, rz_in_ys), r2; itp_opts2d...)
@@ -253,7 +253,7 @@ function equilibrium_solver(input::InverseRunInput)
         # (Numerical operations may have broken exact periodicity)
         @views spl_fs[end, :] .= spl_fs[1, :]
 
-        spl = cubic_interp(spl_xs, Series(spl_fs); bc=PeriodicBC(; check=false))
+        spl = cubic_interp(spl_xs, Series(spl_fs); bc=PeriodicBC())
         spl_fsi = FastInterpolations.cumulative_integrate(spl)
 
         spl_xs .= spl_fsi[:, 5] ./ spl_fsi[mtheta+1, 5]
@@ -268,7 +268,7 @@ function equilibrium_solver(input::InverseRunInput)
         # then evaluate at the uniform SFL theta grid (rzphi_ys). This correctly
         # propagates the SFL coordinate transformation into the rzphi splines.
         # (Using spl.y directly would give pre-transformation values — wrong for eqfun.)
-        spl_post = cubic_interp(spl_xs, Series(spl_fs); bc=PeriodicBC(; check=false))
+        spl_post = cubic_interp(spl_xs, Series(spl_fs); bc=PeriodicBC())
         hint_post = Ref(1)
         for itheta in 0:mtheta
             spl_post(spl_post_buf, rzphi_ys[itheta+1]; hint=hint_post)
@@ -312,6 +312,8 @@ function equilibrium_solver(input::InverseRunInput)
     # Create 2D interpolants for geometric quantities (rzphi)
     rzphi_grid2d = (rzphi_xs, theta_range)
 
+    @views rzphi_fs[:, end, :] .= rzphi_fs[:, 1, :]
+
     rzphi_rsquared = cubic_interp(rzphi_grid2d, rzphi_fs[:, :, 1]; itp_opts2d...)
     rzphi_offset = cubic_interp(rzphi_grid2d, rzphi_fs[:, :, 2]; itp_opts2d...)
     rzphi_nu = cubic_interp(rzphi_grid2d, rzphi_fs[:, :, 3]; itp_opts2d...)
@@ -365,6 +367,7 @@ function equilibrium_solver(input::InverseRunInput)
     end
     # Create 2D interpolants for physics quantities (eqfun)
     eqfun_grid2d = (eqfun_xs, theta_range)
+    @views eqfun_fs[:, end, :] .= eqfun_fs[:, 1, :]
     eqfun_B = cubic_interp(eqfun_grid2d, eqfun_fs[:, :, 1]; itp_opts2d...)
     eqfun_metric1 = cubic_interp(eqfun_grid2d, eqfun_fs[:, :, 2]; itp_opts2d...)
     eqfun_metric2 = cubic_interp(eqfun_grid2d, eqfun_fs[:, :, 3]; itp_opts2d...)

src/Equilibrium/ReadEquilibrium.jl

@@ -210,8 +210,8 @@ function read_chease_binary(config::EquilibriumConfig)
         # Create separate interpolants for R and Z coordinates
         rz_in_xs = xs
         rz_in_ys = range(0, 1; length=mtau) |> collect
-        rz_in_R = cubic_interp((rz_in_xs, rz_in_ys), fs_2d[:, :, 1]; bc=(CubicFit(), PeriodicBC(; check=false)), extrap=(ExtendExtrap(), WrapExtrap()))
-        rz_in_Z = cubic_interp((rz_in_xs, rz_in_ys), fs_2d[:, :, 2]; bc=(CubicFit(), PeriodicBC(; check=false)), extrap=(ExtendExtrap(), WrapExtrap()))
+        rz_in_R = cubic_interp((rz_in_xs, rz_in_ys), fs_2d[:, :, 1]; bc=(CubicFit(), PeriodicBC()), extrap=(ExtendExtrap(), WrapExtrap()))
+        rz_in_Z = cubic_interp((rz_in_xs, rz_in_ys), fs_2d[:, :, 2]; bc=(CubicFit(), PeriodicBC()), extrap=(ExtendExtrap(), WrapExtrap()))
 
         @info "Finished reading CHEASE equilibrium (Binary)"
         return InverseRunInput(config, sq_in, rz_in_xs, rz_in_ys, rz_in_R, rz_in_Z, ro, zo, psio)
@@ -373,7 +373,10 @@ function read_chease_ascii(config::EquilibriumConfig)
     # Create separate interpolants for R and Z coordinates
     rz_in_xs = xs
 
-    opts2d = (bc=(CubicFit(), PeriodicBC(; check=false)), extrap=(ExtendExtrap(), WrapExtrap()))
+    @views R_data[:, end] .= R_data[:, 1]
+    @views Z_data[:, end] .= Z_data[:, 1]
+
+    opts2d = (bc=(CubicFit(), PeriodicBC()), extrap=(ExtendExtrap(), WrapExtrap()))
     rz_in_R = cubic_interp((rz_in_xs, rz_in_ys), R_data; opts2d...)
     rz_in_Z = cubic_interp((rz_in_xs, rz_in_ys), Z_data; opts2d...)
     @info "Finished reading CHEASE equilibrium. Magnetic axis at (ro=$(@sprintf("%.3f", ro)), zo=$(@sprintf("%.3f", zo))), psio=$(@sprintf("%.3e", psio))"
@@ -430,6 +433,7 @@ function read_imas(config::EquilibriumConfig, dd)
     # Extract 1D profiles, converting psi from COCOS 11 to internal
     psi_1d = eqt.profiles_1d.psi .* cocos_factor
     f_1d = eqt.profiles_1d.f          # F(ψ) = R·Bt [T·m], COCOS-independent
+    bt_sign = Int(sign(f_1d[end]))    # sign of toroidal field (before abs is applied below)
     p_1d = eqt.profiles_1d.pressure   # plasma pressure P(ψ) [Pa], COCOS-independent
     q_1d = eqt.profiles_1d.q          # safety factor, COCOS-independent
 
@@ -479,5 +483,5 @@ function read_imas(config::EquilibriumConfig, dd)
           "\n    R ∈ [$(round(rmin; sigdigits=4)), $(round(rmax; sigdigits=4))] m" *
           "\n    Z ∈ [$(round(zmin; sigdigits=4)), $(round(zmax; sigdigits=4))] m"
 
-    return DirectRunInput(config, sq_in, psi_in, psi_in_xs, psi_in_ys, rmin, rmax, zmin, zmax, psio)
+    return DirectRunInput(config, sq_in, psi_in, psi_in_xs, psi_in_ys, rmin, rmax, zmin, zmax, psio, bt_sign)
 end

src/ForceFreeStates/Bal.jl

@@ -289,7 +289,7 @@ function prepare_ballooning_coefficients(ipsi::Int, plasma_eq::Equilibrium.Plasm
     end
 
     # compute curvature terms using native cubic_interp with PeriodicBC
-    spl0_interp = cubic_interp(theta_grid, Series(hcat(1 ./ bsq, jac .* b1 ./ bsq)); bc=PeriodicBC(; check=false))
+    spl0_interp = cubic_interp(theta_grid, Series(hcat(1 ./ bsq, jac .* b1 ./ bsq)); bc=PeriodicBC())
     spl0_d1 = deriv1(spl0_interp)
     # Evaluate derivatives at all theta points (returns Vector of Vectors, stack to matrix)
     spl0_fs1 = stack(spl0_d1(theta_grid))
@@ -306,11 +306,11 @@ function prepare_ballooning_coefficients(ipsi::Int, plasma_eq::Equilibrium.Plasm
     bf_fs[:, 3] = (dbdb2 - dbdb1 .^ 2 ./ (4 .* dbdb0)) ./ (bsq .* jacfac)
     bf_fs[:, 4] = 2 .* kappan .* pressure_gradient ./ chi_prime .* jacfac
     bf_fs[:, 5] = -2 .* kappas .* pressure_gradient ./ chi_prime .* q_derivative .* jacfac
-    ode_coeff_interp = cubic_interp(theta_grid, Series(bf_fs); bc=PeriodicBC(; check=false))
+    ode_coeff_interp = cubic_interp(theta_grid, Series(bf_fs); bc=PeriodicBC())
 
     # initialize bg values
     spl0_bg_fs = -pressure_gradient .* q_derivative .* two_pi_f ./ (bsq .* chi_prime^2)
-    itp_bg = cubic_interp(theta_grid, spl0_bg_fs; bc=PeriodicBC(; check=false))
+    itp_bg = cubic_interp(theta_grid, spl0_bg_fs; bc=PeriodicBC())
     spl0_bg_fsi = FastInterpolations.cumulative_integrate(itp_bg)
     bg_fs = zeros(mtheta + 1, 5)
     bg_fs[:, 5] = spl0_bg_fsi .- spl0_bg_fsi[end]
@@ -329,7 +329,7 @@ function prepare_ballooning_coefficients(ipsi::Int, plasma_eq::Equilibrium.Plasm
         spl1_fs[itheta, 4] = -spl1_fs[itheta, 1]
     end
 
-    itp_spl1 = cubic_interp(theta_grid, Series(spl1_fs); bc=PeriodicBC(; check=false))
+    itp_spl1 = cubic_interp(theta_grid, Series(spl1_fs); bc=PeriodicBC())
     spl1_totals = FastInterpolations.integrate(itp_spl1)
 
     d0bar = [spl1_totals[1] spl1_totals[2]; spl1_totals[3] spl1_totals[4]]
@@ -338,7 +338,7 @@ function prepare_ballooning_coefficients(ipsi::Int, plasma_eq::Equilibrium.Plasm
 
     # if stable by Mercier, no more work to do - return empty asymptotic data
     if di > 0
-        asymptotic_interp = cubic_interp(theta_grid, Series(bg_fs); bc=PeriodicBC(; check=false))
+        asymptotic_interp = cubic_interp(theta_grid, Series(bg_fs); bc=PeriodicBC())
         return di, NaN, ode_coeff_interp, asymptotic_interp, zeros(2, 2), 0.0
     end
 
@@ -363,13 +363,13 @@ function prepare_ballooning_coefficients(ipsi::Int, plasma_eq::Equilibrium.Plasm
         spl2_fs[itheta, 4] = spl1_fs[itheta, 3] * v0[1, 2] + (spl1_fs[itheta, 4] + alpha) * v0[2, 2]
     end
 
-    itp_spl2 = cubic_interp(theta_grid, Series(spl2_fs); bc=PeriodicBC(; check=false))
+    itp_spl2 = cubic_interp(theta_grid, Series(spl2_fs); bc=PeriodicBC())
     spl2_fsi = FastInterpolations.cumulative_integrate(itp_spl2)
     # CRITICAL: Use integrated values as the new fs (matching Fortran's spl2%fs=spl2%fsi)
     spl2_fs_new = copy(spl2_fsi)
 
     # Compute derivatives of spl1 for second-order terms
-    spl1_interp = cubic_interp(theta_grid, Series(spl1_fs); bc=PeriodicBC(; check=false))
+    spl1_interp = cubic_interp(theta_grid, Series(spl1_fs); bc=PeriodicBC())
     spl1_d1 = deriv1(spl1_interp)
     # Evaluate derivatives at all theta points (returns Vector of Vectors, stack to matrix)
     spl1_fs1 = stack(spl1_d1(theta_grid))
@@ -403,7 +403,7 @@ function prepare_ballooning_coefficients(ipsi::Int, plasma_eq::Equilibrium.Plasm
                              d1_21 * v0[1, 2] + d1_22 * v0[2, 2]
     end
 
-    itp_spl3 = cubic_interp(theta_grid, Series(spl3_fs); bc=PeriodicBC(; check=false))
+    itp_spl3 = cubic_interp(theta_grid, Series(spl3_fs); bc=PeriodicBC())
     spl3_totals = FastInterpolations.integrate(itp_spl3)
 
     # Compute first-order constants for both eigenfunctions
@@ -432,7 +432,7 @@ function prepare_ballooning_coefficients(ipsi::Int, plasma_eq::Equilibrium.Plasm
     bg_fs[:, 3] = spl2_fs_new[:, 3] .+ v10[1, 2]
     bg_fs[:, 4] = spl2_fs_new[:, 4] .+ v10[2, 2]
 
-    asymptotic_interp = cubic_interp(theta_grid, Series(bg_fs); bc=PeriodicBC(; check=false))
+    asymptotic_interp = cubic_interp(theta_grid, Series(bg_fs); bc=PeriodicBC())
 
     reference_angle = 0.0 # Central poloidal angle (typically 0)
     return di, alpha, ode_coeff_interp, asymptotic_interp, v0, reference_angle

src/ForceFreeStates/Mercier.jl

@@ -61,7 +61,7 @@ function mercier_scan!(locstab_fs::Matrix{Float64}, plasma_eq::Equilibrium.Plasm
         end
 
         # Integrate quantities with respect to theta using FastInterpolations
-        itp = cubic_interp(plasma_eq.rzphi_ys, Series(ff_fs); bc=PeriodicBC(; check=false))
+        itp = cubic_interp(plasma_eq.rzphi_ys, Series(ff_fs); bc=PeriodicBC())
         avg = FastInterpolations.integrate(itp)
 
         # Evaluate Mercier criterion and related quantities

src/Vacuum/DataTypes.jl

@@ -226,9 +226,9 @@ function PlasmaGeometry(inputs::VacuumInput)
     θ_out = range(; start=0, length=inputs.mtheta, step=2π/inputs.mtheta) # VACUUM uses [0, 2π) grid
 
     # Use one-shot API with PeriodicBC
-    x = cubic_interp(θ_in, inputs.x, θ_out; bc=PeriodicBC(; check=false)) # no endpoint handling needed!
-    z = cubic_interp(θ_in, inputs.z, θ_out; bc=PeriodicBC(; check=false))
-    ν = cubic_interp(θ_in, inputs.ν, θ_out; bc=PeriodicBC(; check=false))
+    x = cubic_interp(θ_in, inputs.x, θ_out; bc=PeriodicBC()) # no endpoint handling needed!
+    z = cubic_interp(θ_in, inputs.z, θ_out; bc=PeriodicBC())
+    ν = cubic_interp(θ_in, inputs.ν, θ_out; bc=PeriodicBC())
 
     return PlasmaGeometry(x, z, ν)
 end
@@ -346,7 +346,7 @@ function PlasmaGeometry3D(inputs::VacuumInput)
         ζ_flat = repeat(collect(ζ_grid); inner=mtheta)
         grid_points = (θ_flat, ζ_flat)
         for (k, data) in enumerate((inputs.x, inputs.y, inputs.z))
-            itp = cubic_interp((θ_in, ζ_in), reshape(data, inputs.mtheta_in, inputs.nzeta_in); bc=(PeriodicBC(; check=false), PeriodicBC(; check=false)))
+            itp = cubic_interp((θ_in, ζ_in), reshape(data, inputs.mtheta_in, inputs.nzeta_in); bc=(PeriodicBC(), PeriodicBC()))
             r[:, k] = itp(grid_points)
         end
     end