From 720d6856532259f7b4ee6bfa8d711f5521a920ff Mon Sep 17 00:00:00 2001 From: logan-nc <6198372+logan-nc@users.noreply.github.com> Date: Fri, 13 Feb 2026 13:52:09 -0500 Subject: [PATCH 1/8] SPLINES - IMPROVEMENT - Add shared 2D hints for hot loop optimization Add shared (Ref{Int}, Ref{Int}) hints to all hot loops with 2D bicubic spline evaluations for O(1) interval search optimization: - Add rzphi_hint field to OdeState for future use in ODE integration - Free.jl: compute_vacuum_inputs loop - VacuumFromEquilibrium.jl: extract_plasma_surface_at_psi loop - SingularCoupling.jl: compute_current_density and compute_surface_area loops - Equilibrium.jl: flux derivative computation loop and equilibrium_separatrix_find Newton iteration while loops Implements FastInterpolations v0.2.10 shared hint optimization for 2D splines. Co-Authored-By: Claude Sonnet 4.5 --- src/Equilibrium/Equilibrium.jl | 27 +++-- src/ForceFreeStates/ForceFreeStatesStructs.jl | 3 + src/ForceFreeStates/Free.jl | 7 +- src/PerturbedEquilibrium/SingularCoupling.jl | 114 ++++++++++-------- src/Vacuum/VacuumFromEquilibrium.jl | 7 +- 5 files changed, 90 insertions(+), 68 deletions(-) diff --git a/src/Equilibrium/Equilibrium.jl b/src/Equilibrium/Equilibrium.jl index 40b5789b..25bccdf1 100644 --- a/src/Equilibrium/Equilibrium.jl +++ b/src/Equilibrium/Equilibrium.jl @@ -107,11 +107,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 +144,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) @@ -446,12 +448,13 @@ function equilibrium_gse!(equil::PlasmaEquilibrium) flux2 = cubic_interp((equil.rzphi_xs, equil.rzphi_ys), flux_fs[:, :, 2]; bc=(CubicFit(), Spl.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 diff --git a/src/ForceFreeStates/ForceFreeStatesStructs.jl b/src/ForceFreeStates/ForceFreeStatesStructs.jl index a86b0059..41907fbd 100644 --- a/src/ForceFreeStates/ForceFreeStatesStructs.jl +++ b/src/ForceFreeStates/ForceFreeStatesStructs.jl @@ -466,6 +466,9 @@ and a small set of temporary matrices and factors used to compute singular-layer spline_hint::Base.RefValue{Int} = Ref(1) # Separate hint for wvmat splines (different grid size than equilibrium profiles) wv_hint::Base.RefValue{Int} = Ref(1) + # Shared 2D hint for CubicInterpolantND (rzphi splines) during ODE integration + # Tuple of (psi_hint, theta_hint) for O(1) interval lookups in 2D bicubic splines + rzphi_hint::Tuple{Base.RefValue{Int},Base.RefValue{Int}} = (Ref(1), Ref(1)) end # Initialize function for OdeState with relevant parameters for array initialization diff --git a/src/ForceFreeStates/Free.jl b/src/ForceFreeStates/Free.jl index 9dd14a14..9ebbbdc7 100644 --- a/src/ForceFreeStates/Free.jl +++ b/src/ForceFreeStates/Free.jl @@ -158,11 +158,12 @@ function compute_vacuum_inputs(psifac::Float64, n::Int, equil::Equilibrium.Plasm # Compute r, z, and ν at the plasma boundary qa = equil.profiles.q_spline(psifac) + hint2d = (Ref(1), Ref(1)) # Shared 2D hint for hot loop optimization for itheta in 1:mtheta # Evaluate geometric quantities at (ψ, θ) - r2 = equil.rzphi_rsquared((psifac, theta_norm[itheta])) - offset = equil.rzphi_offset((psifac, theta_norm[itheta])) - nu_val = equil.rzphi_nu((psifac, theta_norm[itheta])) + r2 = equil.rzphi_rsquared((psifac, theta_norm[itheta]); hint=hint2d) + offset = equil.rzphi_offset((psifac, theta_norm[itheta]); hint=hint2d) + nu_val = equil.rzphi_nu((psifac, theta_norm[itheta]); hint=hint2d) rfac[itheta] = sqrt(r2) angle[itheta] = 2π * (theta_norm[itheta] + offset) diff --git a/src/PerturbedEquilibrium/SingularCoupling.jl b/src/PerturbedEquilibrium/SingularCoupling.jl index 99d25c70..9f2a2512 100644 --- a/src/PerturbedEquilibrium/SingularCoupling.jl +++ b/src/PerturbedEquilibrium/SingularCoupling.jl @@ -28,28 +28,32 @@ Calculates the coupling between forcing modes and resonant surfaces, which deter the effectiveness of external perturbations at rational surfaces where q = m/n. Implements GPEC algorithm from gpout.f: -1. For each forcing mode, apply unit field and compute plasma response -2. For each singular surface, evaluate field jump and compute coupling quantities -3. Calculate diagnostic island widths and overlap parameters + + 1. For each forcing mode, apply unit field and compute plasma response + 2. For each singular surface, evaluate field jump and compute coupling quantities + 3. Calculate diagnostic island widths and overlap parameters ## Arguments -- `state`: Output state structure to store results -- `equil`: Equilibrium solution with q-profile and flux surfaces -- `ForceFreeStates_results`: ForceFreeStates stability calculation results -- `vac_data`: Vacuum response data including Green's functions -- `ffs_intr`: ForceFreeStates internal state with singular surface data -- `intr`: PerturbedEquilibrium internal state with permeability matrix -- `ctrl`: Control parameters + + - `state`: Output state structure to store results + - `equil`: Equilibrium solution with q-profile and flux surfaces + - `ForceFreeStates_results`: ForceFreeStates stability calculation results + - `vac_data`: Vacuum response data including Green's functions + - `ffs_intr`: ForceFreeStates internal state with singular surface data + - `intr`: PerturbedEquilibrium internal state with permeability matrix + - `ctrl`: Control parameters ## Modifies + Populates the following fields in `state`: -- `resonant_flux[msing, numpert_total]`: Normalized resonant flux Φ_r/A -- `resonant_current[msing, numpert_total]`: Resonant current density -- `island_width_sq[msing, numpert_total]`: Square of island half-width (w/2)² -- `penetrated_field[msing, numpert_total]`: Normal field at resonant surface -- `delta_prime[msing, numpert_total]`: Tearing stability parameter Δ' -- `island_half_width[msing]`: Dimensional island half-width -- `chirikov_parameter[msing]`: Island overlap metric + + - `resonant_flux[msing, numpert_total]`: Normalized resonant flux Φ_r/A + - `resonant_current[msing, numpert_total]`: Resonant current density + - `island_width_sq[msing, numpert_total]`: Square of island half-width (w/2)² + - `penetrated_field[msing, numpert_total]`: Normal field at resonant surface + - `delta_prime[msing, numpert_total]`: Tearing stability parameter Δ' + - `island_half_width[msing]`: Dimensional island half-width + - `chirikov_parameter[msing]`: Island overlap metric Note: numpert_total = mpert × npert handles all (m,n) mode combinations """ @@ -120,7 +124,7 @@ function compute_singular_coupling_metrics!( nhigh = ffs_intr.nhigh # Perturbed equilibrium calculations always use nowall - wall_settings = Vacuum.WallShapeSettings(shape="nowall") + wall_settings = Vacuum.WallShapeSettings(; shape="nowall") if ctrl.verbose println(" Processing toroidal modes n = $nlow:$nhigh") @@ -311,7 +315,7 @@ Interpolate magnetic field at arbitrary flux surface. Interpolates displacement from ForceFreeStates solution and converts to field using ideal MHD relation from flux coordinates: - b^ψ = i * χ₁ * (m - n*q) * ξ_ψ +b^ψ = i * χ₁ * (m - n*q) * ξ_ψ This matches the field reconstruction in FieldReconstruction.jl. """ @@ -367,10 +371,10 @@ end Compute effective current density coefficient at given flux surface. Implements GPEC's j_c calculation (gpout.f line 560-578): - j_c = χ₁² * q / (μ₀ * integral) +j_c = χ₁² * q / (μ₀ * integral) where the integral is computed via flux surface integration: - integral = ∫ (jac * |∇ψ| * sqreqb / |∇ψ|³) dθ +integral = ∫ (jac * |∇ψ| * sqreqb / |∇ψ|³) dθ ## GPEC Formula @@ -394,9 +398,10 @@ j_c(ising)=1.0/j_c(ising)*chi1**2*sq%f(4)/mu0 Uses trapezoidal rule integration around the flux surface with metric quantities from the equilibrium bicubic spline (rzphi). The integrand includes: -- jac: Jacobian of flux coordinates -- |∇ψ|: Flux gradient magnitude (delpsi) -- sqreqb: Magnetic field quantity (F² + χ₁²|∇ψ|²)/(2πR)² + + - jac: Jacobian of flux coordinates + - |∇ψ|: Flux gradient magnitude (delpsi) + - sqreqb: Magnetic field quantity (F² + χ₁²|∇ψ|²)/(2πR)² """ function compute_current_density( equil::Equilibrium.PlasmaEquilibrium, @@ -427,17 +432,18 @@ function compute_current_density( last_delpsi = 0.0 last_jcfun = 0.0 + hint2d = (Ref(1), Ref(1)) # Shared 2D hint for hot loop optimization for itheta in 0:mthsurf # Theta coordinate normalized to [0, 1] theta = itheta / mthsurf # Evaluate bicubic splines with derivatives at (psi, theta) # New API uses separate interpolants for each component - r2 = equil.rzphi_rsquared((psi, theta)) # rfac² - deta = equil.rzphi_offset((psi, theta)) # angle offset - jac = equil.rzphi_jac((psi, theta)) # Jacobian - r2_y = equil.rzphi_rsquared((psi, theta); deriv=(0,1)) # ∂(rfac²)/∂theta - deta_y = equil.rzphi_offset((psi, theta); deriv=(0,1)) # ∂(deta)/∂theta + r2 = equil.rzphi_rsquared((psi, theta); hint=hint2d) # rfac² + deta = equil.rzphi_offset((psi, theta); hint=hint2d) # angle offset + jac = equil.rzphi_jac((psi, theta); hint=hint2d) # Jacobian + r2_y = equil.rzphi_rsquared((psi, theta); deriv=(0, 1), hint=hint2d) # ∂(rfac²)/∂theta + deta_y = equil.rzphi_offset((psi, theta); deriv=(0, 1), hint=hint2d) # ∂(deta)/∂theta rfac = sqrt(abs(r2)) fy_rfac2 = r2_y @@ -493,11 +499,13 @@ Apply Green's function to Fourier mode coefficients. Simplified version that extracts only plasma surface rows (first mtheta rows). ## Arguments -- `green`: Green's function matrix [2*mtheta, 2*mpert] -- `mode_coeffs`: Complex Fourier coefficients [mpert] + + - `green`: Green's function matrix [2*mtheta, 2*mpert] + - `mode_coeffs`: Complex Fourier coefficients [mpert] ## Returns -- Potential at plasma surface theta points [mtheta] + + - Potential at plasma surface theta points [mtheta] """ function apply_green_function_simple( green::Matrix{Float64}, @@ -509,7 +517,7 @@ function apply_green_function_simple( # Pack complex to real/imag format packed = zeros(Float64, 2 * mpert) for i in 1:mpert - packed[2*i - 1] = real(mode_coeffs[i]) + packed[2*i-1] = real(mode_coeffs[i]) packed[2*i] = imag(mode_coeffs[i]) end @@ -530,21 +538,25 @@ end Compute surface inductance matrix from Green's functions at flux surface. This implements the full GPEC gpvacuum_flxsurf algorithm (gpvacuum.f line 299-331): -1. For each mode, apply Green's functions to unit flux -2. Compute surface current from potential jump: kax = (chi - che) / μ₀ -3. Fourier transform to mode space -4. Return inductance: L_surf = flux * inv(current) + + 1. For each mode, apply Green's functions to unit flux + 2. Compute surface current from potential jump: kax = (chi - che) / μ₀ + 3. Fourier transform to mode space + 4. Return inductance: L_surf = flux * inv(current) ## Arguments -- `grri`: Interior Green's function [2*mtheta, 2*mpert] (stored in sing_surf) -- `grre`: Exterior Green's function [2*mtheta, 2*mpert] (stored in sing_surf) -- `ffs_intr`: ForceFreeStates internal state -- `intr`: PerturbedEquilibrium internal state with mode arrays + + - `grri`: Interior Green's function [2*mtheta, 2*mpert] (stored in sing_surf) + - `grre`: Exterior Green's function [2*mtheta, 2*mpert] (stored in sing_surf) + - `ffs_intr`: ForceFreeStates internal state + - `intr`: PerturbedEquilibrium internal state with mode arrays ## Returns + Surface inductance matrix [numpert_total, numpert_total] ## GPEC Reference + ```fortran DO i=1,mpert vbwp_mn=0; vbwp_mn(i)=1.0 @@ -650,7 +662,7 @@ end Compute singular flux from resonant current using surface inductance. Uses surface inductance matrix at the singular surface: - Φ = L_surf * I +Φ = L_surf * I where L_surf is computed from Green's functions stored in sing_surf. @@ -718,7 +730,7 @@ end Compute flux surface area at given ψ. Implements GPEC's area calculation (gpout.f line 568): - area = ∫ jac * |∇ψ| dθ +area = ∫ jac * |∇ψ| dθ where the integral is computed around the flux surface. @@ -740,8 +752,9 @@ area(ising)=area(ising)-jac*delpsi(mthsurf)/mthsurf ! trapezoidal rule ## Implementation Uses trapezoidal rule integration around the flux surface with: -- jac: Jacobian of flux coordinates from rzphi -- |∇ψ|: Flux gradient magnitude (delpsi) from metric tensor + + - jac: Jacobian of flux coordinates from rzphi + - |∇ψ|: Flux gradient magnitude (delpsi) from metric tensor """ function compute_surface_area( equil::Equilibrium.PlasmaEquilibrium, @@ -764,16 +777,17 @@ function compute_surface_area( last_jac = 0.0 last_delpsi = 0.0 + hint2d = (Ref(1), Ref(1)) # Shared 2D hint for hot loop optimization for itheta in 0:mthsurf # Theta coordinate normalized to [0, 1] theta = itheta / mthsurf # Evaluate bicubic splines with derivatives at (psi, theta) - r2 = equil.rzphi_rsquared((psi, theta)) - jac = equil.rzphi_jac((psi, theta)) - deta = equil.rzphi_offset((psi, theta)) - r2_y = equil.rzphi_rsquared((psi, theta); deriv=(0,1)) - deta_y = equil.rzphi_offset((psi, theta); deriv=(0,1)) + r2 = equil.rzphi_rsquared((psi, theta); hint=hint2d) + jac = equil.rzphi_jac((psi, theta); hint=hint2d) + deta = equil.rzphi_offset((psi, theta); hint=hint2d) + r2_y = equil.rzphi_rsquared((psi, theta); deriv=(0, 1), hint=hint2d) + deta_y = equil.rzphi_offset((psi, theta); deriv=(0, 1), hint=hint2d) # Compute rfac rfac = sqrt(abs(r2)) diff --git a/src/Vacuum/VacuumFromEquilibrium.jl b/src/Vacuum/VacuumFromEquilibrium.jl index 296621a6..fc4d6dbf 100644 --- a/src/Vacuum/VacuumFromEquilibrium.jl +++ b/src/Vacuum/VacuumFromEquilibrium.jl @@ -70,15 +70,16 @@ function extract_plasma_surface_at_psi( # Evaluate equilibrium around the flux surface twopi = 2π + hint2d = (Ref(1), Ref(1)) # Shared 2D hint for hot loop optimization for itheta in 0:mtheta_eq-1 # Theta coordinate normalized to [0, 1) theta = itheta / mtheta_eq # Evaluate bicubic splines at (psi, theta) # New API uses separate interpolants for each component - r2 = equil.rzphi_rsquared((psi, theta)) # r² or rfac² - deta = equil.rzphi_offset((psi, theta)) # angle offset - dphi = equil.rzphi_nu((psi, theta)) # toroidal angle offset (nu) + r2 = equil.rzphi_rsquared((psi, theta); hint=hint2d) # r² or rfac² + deta = equil.rzphi_offset((psi, theta); hint=hint2d) # angle offset + dphi = equil.rzphi_nu((psi, theta); hint=hint2d) # toroidal angle offset (nu) rfac = sqrt(abs(r2)) From 0a47ccd18a95c13ff88d40acbb4d503693e22637 Mon Sep 17 00:00:00 2001 From: logan-nc <6198372+logan-nc@users.noreply.github.com> Date: Fri, 13 Feb 2026 13:53:17 -0500 Subject: [PATCH 2/8] SPLINES - IMPROVEMENT - Convert to PeriodicBC(endpoint=:exclusive) API MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Use FastInterpolations v0.2.9+ PeriodicBC(endpoint=:exclusive) instead of manually appending endpoints with vcat for periodic boundary conditions: - Kernel2D.jl: Remove vcat() endpoint duplication, use period=2π - MathUtils.jl: Infer period from grid and use exclusive endpoint This eliminates redundant endpoint storage and leverages native periodic spline support in FastInterpolations. Co-Authored-By: Claude Sonnet 4.5 --- src/Vacuum/Kernel2D.jl | 10 ++++------ src/Vacuum/MathUtils.jl | 11 ++++++----- 2 files changed, 10 insertions(+), 11 deletions(-) diff --git a/src/Vacuum/Kernel2D.jl b/src/Vacuum/Kernel2D.jl index bce1b7a7..0a97f9e2 100644 --- a/src/Vacuum/Kernel2D.jl +++ b/src/Vacuum/Kernel2D.jl @@ -110,12 +110,10 @@ function kernel!( log_correction_2=4.0*dtheta*(7.0*log(2*dtheta)-11.0/15.0)/45.0 # Used for Z'_θ and X'_θ in eq.(51) - # Close the loop for periodic BC by appending first point at the end - theta_closed = vcat(collect(theta_grid), theta_grid[end] + dtheta) - x_closed = vcat(x_sourcepoints, x_sourcepoints[1]) - z_closed = vcat(z_sourcepoints, z_sourcepoints[1]) - spline_x = cubic_interp(theta_closed, x_closed; bc=PeriodicBC()) - spline_z = cubic_interp(theta_closed, z_closed; bc=PeriodicBC()) + # Use exclusive endpoint PeriodicBC (FastInterpolations v0.2.9+) + theta_vec = collect(theta_grid) + spline_x = cubic_interp(theta_vec, x_sourcepoints; bc=PeriodicBC(; endpoint=:exclusive, period=2π)) + spline_z = cubic_interp(theta_vec, z_sourcepoints; bc=PeriodicBC(; endpoint=:exclusive, period=2π)) # Create derivative views once, reuse for all evaluations (avoids allocation per call) d1_spline_x = deriv1(spline_x) d1_spline_z = deriv1(spline_z) diff --git a/src/Vacuum/MathUtils.jl b/src/Vacuum/MathUtils.jl index 0df35dae..ccf03728 100644 --- a/src/Vacuum/MathUtils.jl +++ b/src/Vacuum/MathUtils.jl @@ -4,12 +4,13 @@ # Returns the array of derivatives at all x points, I think this acts like difspl # in the Fortran but need to check/consolidate spline routines later function periodic_cubic_deriv(theta, vals) - # Close the loop for periodic BC by appending first point at the end - theta_closed = vcat(collect(theta), theta[end] + (theta[2] - theta[1])) - vals_closed = vcat(vals, vals[1]) - spline = cubic_interp(theta_closed, vals_closed; bc=PeriodicBC()) + # Use exclusive endpoint PeriodicBC (FastInterpolations v0.2.9+) + # Infer period from theta grid + theta_vec = collect(theta) + period = theta_vec[end] - theta_vec[1] + (theta_vec[2] - theta_vec[1]) + spline = cubic_interp(theta_vec, vals; bc=PeriodicBC(; endpoint=:exclusive, period=period)) d1 = deriv1(spline) - return d1.(collect(theta)) + return d1.(theta_vec) end """ From bb54d91b1eccaed9a548815fa9d8c5fe606958e0 Mon Sep 17 00:00:00 2001 From: logan-nc <6198372+logan-nc@users.noreply.github.com> Date: Fri, 13 Feb 2026 13:54:51 -0500 Subject: [PATCH 3/8] SPLINES - IMPROVEMENT - Use native FastInterpolations integration functions Replace custom Spl.cumulative_integral and Spl.total_integral with native FastInterpolations.cumulative_integrate() and FastInterpolations.integrate(): - Equilibrium.jl: flux integral computation - ReadEquilibrium.jl: pressure integral normalization - InverseEquilibrium.jl: periodic spline integration - Mercier.jl: theta-averaged quantities - Bal.jl: ballooning mode integrals Native integration functions provide exact spline integration using cubic spline coefficients without requiring custom wrapper functions. Co-Authored-By: Claude Sonnet 4.5 --- src/Equilibrium/Equilibrium.jl | 5 +++-- src/Equilibrium/InverseEquilibrium.jl | 2 +- src/Equilibrium/ReadEquilibrium.jl | 7 ++++--- src/ForceFreeStates/Bal.jl | 12 ++++++++---- src/ForceFreeStates/Mercier.jl | 5 +++-- 5 files changed, 19 insertions(+), 12 deletions(-) diff --git a/src/Equilibrium/Equilibrium.jl b/src/Equilibrium/Equilibrium.jl index 25bccdf1..606dfe2d 100644 --- a/src/Equilibrium/Equilibrium.jl +++ b/src/Equilibrium/Equilibrium.jl @@ -495,8 +495,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=Spl.PeriodicBC()) + term[ipsi, :] .= FastInterpolations.integrate(itp) end totali = sum(term; dims=2) diff --git a/src/Equilibrium/InverseEquilibrium.jl b/src/Equilibrium/InverseEquilibrium.jl index e589fd8d..be9507ba 100644 --- a/src/Equilibrium/InverseEquilibrium.jl +++ b/src/Equilibrium/InverseEquilibrium.jl @@ -207,7 +207,7 @@ function equilibrium_solver(input::InverseRunInput) spl_fs[end, :] .= spl_fs[1, :] spl = cubic_interp(spl_xs, spl_fs; bc=Spl.PeriodicBC()) - spl_fsi = Spl.cumulative_integral(spl_xs, spl_fs; bc=Spl.PeriodicBC()) + spl_fsi = FastInterpolations.cumulative_integrate(spl) spl_xs = spl_fsi[:, 5] ./ spl_fsi[mtheta+1, 5] spl_fs[:, 2] .+= rzphi_ys .- spl_xs diff --git a/src/Equilibrium/ReadEquilibrium.jl b/src/Equilibrium/ReadEquilibrium.jl index 1ce8a974..f676a381 100644 --- a/src/Equilibrium/ReadEquilibrium.jl +++ b/src/Equilibrium/ReadEquilibrium.jl @@ -171,8 +171,9 @@ function read_chease_binary(config::EquilibriumConfig) fs[:, 2] .= zcppr fs[:, 3] .= zq - # Compute cumulative integral of pressure column for normalization - fsi_pressure = Spl.cumulative_integral(xs, fs[:, 2]; bc=CubicFit()) + # Compute cumulative integral of pressure column for normalization using FastInterpolations + itp_pressure = cubic_interp(xs, fs[:, 2]; bc=CubicFit()) + fsi_pressure = FastInterpolations.cumulative_integrate(itp_pressure) # Make a writable copy and normalize pressure integral fs_copy = copy(fs) fs_copy[:, 2] .= (fsi_pressure .- fsi_pressure[ma]) .* psio @@ -379,4 +380,4 @@ function read_chease_ascii(config::EquilibriumConfig) println("--> Finished reading CHEASE equilibrium.") println(" Magnetic axis at (ro=$ro, zo=$zo), psio=$psio") return InverseRunInput(config, sq_in, rz_in_xs, rz_in_ys, rz_in_R, rz_in_Z, ro, zo, psio) -end \ No newline at end of file +end diff --git a/src/ForceFreeStates/Bal.jl b/src/ForceFreeStates/Bal.jl index 1f858319..52f81b9a 100644 --- a/src/ForceFreeStates/Bal.jl +++ b/src/ForceFreeStates/Bal.jl @@ -310,7 +310,8 @@ function prepare_ballooning_coefficients(ipsi::Int, plasma_eq::Equilibrium.Plasm # initialize bg values spl0_bg_fs = -pressure_gradient .* q_derivative .* two_pi_f ./ (bsq .* chi_prime^2) - spl0_bg_fsi = Spl.cumulative_integral(theta_grid, spl0_bg_fs; bc=PeriodicBC()) + 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] @@ -328,7 +329,8 @@ function prepare_ballooning_coefficients(ipsi::Int, plasma_eq::Equilibrium.Plasm spl1_fs[itheta, 4] = -spl1_fs[itheta, 1] end - spl1_totals = Spl.total_integral(theta_grid, spl1_fs; bc=PeriodicBC()) + itp_spl1 = cubic_interp(theta_grid, spl1_fs; bc=PeriodicBC()) + spl1_totals = FastInterpolations.integrate(itp_spl1) d0bar = [spl1_totals[1] spl1_totals[2]; spl1_totals[3] spl1_totals[4]] @@ -361,7 +363,8 @@ 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 - spl2_fsi = Spl.cumulative_integral(theta_grid, spl2_fs; bc=PeriodicBC()) + itp_spl2 = cubic_interp(theta_grid, 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) @@ -400,7 +403,8 @@ function prepare_ballooning_coefficients(ipsi::Int, plasma_eq::Equilibrium.Plasm d1_21 * v0[1, 2] + d1_22 * v0[2, 2] end - spl3_totals = Spl.total_integral(theta_grid, spl3_fs; bc=PeriodicBC()) + itp_spl3 = cubic_interp(theta_grid, spl3_fs; bc=PeriodicBC()) + spl3_totals = FastInterpolations.integrate(itp_spl3) # Compute first-order constants for both eigenfunctions # First eigenfunction (with -alpha correction) diff --git a/src/ForceFreeStates/Mercier.jl b/src/ForceFreeStates/Mercier.jl index 79f6621e..cefd0a81 100644 --- a/src/ForceFreeStates/Mercier.jl +++ b/src/ForceFreeStates/Mercier.jl @@ -60,8 +60,9 @@ function mercier_scan!(locstab_fs::Matrix{Float64}, plasma_eq::Equilibrium.Plasm @views ff_fs[itheta, :] .*= jac / v1 end - # Integrate quantities with respect to theta using exact spline integral - avg = Spl.total_integral(plasma_eq.rzphi_ys, ff_fs; bc=Spl.PeriodicBC()) + # Integrate quantities with respect to theta using FastInterpolations + itp = cubic_interp(plasma_eq.rzphi_ys, ff_fs; bc=Spl.PeriodicBC()) + avg = FastInterpolations.integrate(itp) # Evaluate Mercier criterion and related quantities term = twopif * p1 * v1 / (q1 * chi1^3) * avg[2] From b60981af77dfec84b851873eae6ec999a71f4873 Mon Sep 17 00:00:00 2001 From: logan-nc <6198372+logan-nc@users.noreply.github.com> Date: Fri, 13 Feb 2026 14:24:25 -0500 Subject: [PATCH 4/8] SPLINES - CLEANUP - Remove src/Splines directory, use FastInterpolations directly Remove the Splines module wrapper and use FastInterpolations directly throughout the codebase: - Remove src/Splines/ directory (FastInterpolationsAdaptor.jl, Splines.jl) - Remove import .SplinesMod as Spl from JPEC.jl - Replace all Spl.PeriodicBC() with PeriodicBC() - Replace all Spl.CubicFit() with CubicFit() - Update test/runtests_spline.jl to use native FastInterpolations functions - Remove legacy MultiQuantityProfile and BicubicSpline tests The Splines module was only re-exporting FastInterpolations with thin wrappers, which created unnecessary indirection. Using FastInterpolations directly improves code clarity and leverages all native FastInterpolations v0.2.10 features. Co-Authored-By: Claude Sonnet 4.5 --- src/Equilibrium/AnalyticEquilibrium.jl | 4 +- src/Equilibrium/DirectEquilibrium.jl | 16 +- src/Equilibrium/Equilibrium.jl | 7 +- src/Equilibrium/InverseEquilibrium.jl | 26 +- src/Equilibrium/ReadEquilibrium.jl | 13 +- src/ForceFreeStates/ForceFreeStates.jl | 1 - src/ForceFreeStates/Mercier.jl | 2 +- src/JPEC.jl | 17 +- .../PerturbedEquilibrium.jl | 14 +- src/Splines/FastInterpolationsAdaptor.jl | 223 ------------------ src/Splines/Splines.jl | 18 -- src/Vacuum/Vacuum.jl | 1 - test/runtests_spline.jl | 86 +------ 13 files changed, 60 insertions(+), 368 deletions(-) delete mode 100644 src/Splines/FastInterpolationsAdaptor.jl delete mode 100644 src/Splines/Splines.jl diff --git a/src/Equilibrium/AnalyticEquilibrium.jl b/src/Equilibrium/AnalyticEquilibrium.jl index c4370004..daf49d7f 100644 --- a/src/Equilibrium/AnalyticEquilibrium.jl +++ b/src/Equilibrium/AnalyticEquilibrium.jl @@ -209,9 +209,9 @@ function lar_run(equil_input::EquilibriumConfig, lar_input::LargeAspectRatioConf 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)) + bc=(CubicFit(), 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)) + 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) diff --git a/src/Equilibrium/DirectEquilibrium.jl b/src/Equilibrium/DirectEquilibrium.jl index 8393c89f..c678d59d 100644 --- a/src/Equilibrium/DirectEquilibrium.jl +++ b/src/Equilibrium/DirectEquilibrium.jl @@ -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` @@ -525,13 +525,13 @@ function equilibrium_solver(raw_profile::DirectRunInput) 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)) + bc=(CubicFit(), PeriodicBC()), extrap=(:extension, :wrap)) rzphi_offset = cubic_interp((rzphi_xs, rzphi_ys), rzphi_nodes[:, :, 2]; - bc=(CubicFit(), Spl.PeriodicBC()), extrap=(:extension, :wrap)) + bc=(CubicFit(), PeriodicBC()), extrap=(:extension, :wrap)) rzphi_nu = cubic_interp((rzphi_xs, rzphi_ys), rzphi_nodes[:, :, 3]; - bc=(CubicFit(), Spl.PeriodicBC()), extrap=(:extension, :wrap)) + bc=(CubicFit(), PeriodicBC()), extrap=(:extension, :wrap)) rzphi_jac = cubic_interp((rzphi_xs, rzphi_ys), rzphi_nodes[:, :, 4]; - bc=(CubicFit(), Spl.PeriodicBC()), extrap=(:extension, :wrap)) + 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 @@ -597,11 +597,11 @@ function equilibrium_solver(raw_profile::DirectRunInput) 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)) + bc=(CubicFit(), PeriodicBC()), extrap=(:extension, :wrap)) eqfun_metric1 = cubic_interp((rzphi_xs, rzphi_ys), eqfun_fs_nodes[:, :, 2]; - bc=(CubicFit(), Spl.PeriodicBC()), extrap=(:extension, :wrap)) + bc=(CubicFit(), PeriodicBC()), extrap=(:extension, :wrap)) eqfun_metric2 = cubic_interp((rzphi_xs, rzphi_ys), eqfun_fs_nodes[:, :, 3]; - bc=(CubicFit(), Spl.PeriodicBC()), extrap=(:extension, :wrap)) + bc=(CubicFit(), PeriodicBC()), extrap=(:extension, :wrap)) return PlasmaEquilibrium(raw_profile.config, EquilibriumParameters(), profiles, rzphi_xs, rzphi_ys, diff --git a/src/Equilibrium/Equilibrium.jl b/src/Equilibrium/Equilibrium.jl index 606dfe2d..da017d39 100644 --- a/src/Equilibrium/Equilibrium.jl +++ b/src/Equilibrium/Equilibrium.jl @@ -1,7 +1,6 @@ module Equilibrium # --- Module-level Dependencies --- -import ..Spl using Printf, OrdinaryDiffEq, DiffEqCallbacks, LinearAlgebra, HDF5 using TOML @@ -444,9 +443,9 @@ function equilibrium_gse!(equil::PlasmaEquilibrium) 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)) + bc=(CubicFit(), PeriodicBC()), extrap=(:extension, :wrap)) flux2 = cubic_interp((equil.rzphi_xs, equil.rzphi_ys), flux_fs[:, :, 2]; - bc=(CubicFit(), Spl.PeriodicBC()), extrap=(:extension, :wrap)) + 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 @@ -496,7 +495,7 @@ function equilibrium_gse!(equil::PlasmaEquilibrium) fs_matrix[:, 2] = source[ipsi, :] # Compute total integral using FastInterpolations native integration - itp = cubic_interp(equil.rzphi_ys, fs_matrix; bc=Spl.PeriodicBC()) + itp = cubic_interp(equil.rzphi_ys, fs_matrix; bc=PeriodicBC()) term[ipsi, :] .= FastInterpolations.integrate(itp) end diff --git a/src/Equilibrium/InverseEquilibrium.jl b/src/Equilibrium/InverseEquilibrium.jl index be9507ba..046be5c4 100644 --- a/src/Equilibrium/InverseEquilibrium.jl +++ b/src/Equilibrium/InverseEquilibrium.jl @@ -119,9 +119,9 @@ function equilibrium_solver(input::InverseRunInput) # Create 2D interpolants for r² and dη rz_rsq = cubic_interp((rz_in_xs, rz_in_ys), r2; - bc=(CubicFit(), Spl.PeriodicBC()), extrap=(:extension, :wrap)) + bc=(CubicFit(), PeriodicBC()), extrap=(:extension, :wrap)) rz_deta = cubic_interp((rz_in_xs, rz_in_ys), deta; - bc=(CubicFit(), Spl.PeriodicBC()), extrap=(:extension, :wrap)) + bc=(CubicFit(), PeriodicBC()), extrap=(:extension, :wrap)) # c----------------------------------------------------------------------- # c prepare new spline type for surface quantities. @@ -137,7 +137,7 @@ function equilibrium_solver(input::InverseRunInput) if grid_type == "ldp" sq_xs = psilow .+ (psihigh - psilow) .* (sin.(range(0.0, 1.0; length=mpsi+1) .* (π/2))) .^ 2 sq_fs = zeros(Float64, mpsi+1, 4) - sq = cubic_interp(sq_xs, sq_fs; bc=Spl.CubicFit(), extrap=:extension) + sq = cubic_interp(sq_xs, sq_fs; bc=CubicFit(), extrap=:extension) else error("Only 'ldp' grid_type is implemented for now.") end @@ -206,7 +206,7 @@ function equilibrium_solver(input::InverseRunInput) # (Numerical operations may have broken exact periodicity) spl_fs[end, :] .= spl_fs[1, :] - spl = cubic_interp(spl_xs, spl_fs; bc=Spl.PeriodicBC()) + spl = cubic_interp(spl_xs, spl_fs; bc=PeriodicBC()) spl_fsi = FastInterpolations.cumulative_integrate(spl) spl_xs = spl_fsi[:, 5] ./ spl_fsi[mtheta+1, 5] @@ -226,7 +226,7 @@ function equilibrium_solver(input::InverseRunInput) sq_fs[ipsi+1, 4] = spl_fsi[mtheta+1, 4] * sq_fs[ipsi+1, 1] / (2 * twopi * psio) # q-profile end - sq = cubic_interp(sq_xs, sq_fs; bc=Spl.CubicFit(), extrap=:extension) + sq = cubic_interp(sq_xs, sq_fs; bc=CubicFit(), extrap=:extension) # Evaluate sq and its derivative at all grid points f_sq = zeros(Float64, mpsi+1, 4) @@ -251,18 +251,18 @@ function equilibrium_solver(input::InverseRunInput) sq_fs[ipsi+1, 4] *= ffac rzphi_fs[ipsi+1, :, 3] *= ffac end - sq = cubic_interp(sq_xs, sq_fs; bc=Spl.CubicFit(), extrap=:extension) + sq = cubic_interp(sq_xs, sq_fs; bc=CubicFit(), extrap=:extension) end qa = f_sq[mpsi+1, 4] + f1_sq[mpsi+1, 4] * (1 - sq_xs[mpsi+1]) # Create 2D interpolants for geometric quantities (rzphi) rzphi_rsquared = cubic_interp((rzphi_xs, rzphi_ys), rzphi_fs[:, :, 1]; - bc=(CubicFit(), Spl.PeriodicBC()), extrap=(:extension, :wrap)) + bc=(CubicFit(), PeriodicBC()), extrap=(:extension, :wrap)) rzphi_offset = cubic_interp((rzphi_xs, rzphi_ys), rzphi_fs[:, :, 2]; - bc=(CubicFit(), Spl.PeriodicBC()), extrap=(:extension, :wrap)) + bc=(CubicFit(), PeriodicBC()), extrap=(:extension, :wrap)) rzphi_nu = cubic_interp((rzphi_xs, rzphi_ys), rzphi_fs[:, :, 3]; - bc=(CubicFit(), Spl.PeriodicBC()), extrap=(:extension, :wrap)) + bc=(CubicFit(), PeriodicBC()), extrap=(:extension, :wrap)) rzphi_jac = cubic_interp((rzphi_xs, rzphi_ys), rzphi_fs[:, :, 4]; - bc=(CubicFit(), Spl.PeriodicBC()), extrap=(:extension, :wrap)) + bc=(CubicFit(), PeriodicBC()), extrap=(:extension, :wrap)) for ipsi in 0:mpsi f_sq = sq(sq_xs[ipsi+1]) @@ -312,11 +312,11 @@ function equilibrium_solver(input::InverseRunInput) end # Create 2D interpolants for physics quantities (eqfun) eqfun_B = cubic_interp((eqfun_xs, eqfun_ys), eqfun_fs[:, :, 1]; - bc=(CubicFit(), Spl.PeriodicBC()), extrap=(:extension, :wrap)) + bc=(CubicFit(), PeriodicBC()), extrap=(:extension, :wrap)) eqfun_metric1 = cubic_interp((eqfun_xs, eqfun_ys), eqfun_fs[:, :, 2]; - bc=(CubicFit(), Spl.PeriodicBC()), extrap=(:extension, :wrap)) + bc=(CubicFit(), PeriodicBC()), extrap=(:extension, :wrap)) eqfun_metric2 = cubic_interp((eqfun_xs, eqfun_ys), eqfun_fs[:, :, 3]; - bc=(CubicFit(), Spl.PeriodicBC()), extrap=(:extension, :wrap)) + bc=(CubicFit(), PeriodicBC()), extrap=(:extension, :wrap)) # Create ProfileSplines from sq interpolant # sq_fs columns: [F*2π, P, dV/dψ, q] diff --git a/src/Equilibrium/ReadEquilibrium.jl b/src/Equilibrium/ReadEquilibrium.jl index f676a381..f3945d78 100644 --- a/src/Equilibrium/ReadEquilibrium.jl +++ b/src/Equilibrium/ReadEquilibrium.jl @@ -211,9 +211,9 @@ function read_chease_binary(config::EquilibriumConfig) rz_in_xs = xs rz_in_ys = range(0, 2π; length=mtau) |> collect rz_in_R = cubic_interp((rz_in_xs, rz_in_ys), fs_2d[:, :, 1]; - bc=(CubicFit(), Spl.PeriodicBC()), extrap=(:extension, :wrap)) + bc=(CubicFit(), PeriodicBC()), extrap=(:extension, :wrap)) rz_in_Z = cubic_interp((rz_in_xs, rz_in_ys), fs_2d[:, :, 2]; - bc=(CubicFit(), Spl.PeriodicBC()), extrap=(:extension, :wrap)) + bc=(CubicFit(), PeriodicBC()), extrap=(:extension, :wrap)) println("--> 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) @@ -352,8 +352,9 @@ function read_chease_ascii(config::EquilibriumConfig) fs[:, 2] .= zcppr # normalized Pressure fs[:, 3] .= zq # q profile # Fit spline with extrapolation boundary condition (bctype = 3) - # Compute cumulative integral of pressure column for normalization - fsi_pressure = Spl.cumulative_integral(xs, fs[:, 2]; bc=CubicFit()) + # Compute cumulative integral of pressure column for normalization using FastInterpolations + itp_pressure = cubic_interp(xs, fs[:, 2]; bc=CubicFit()) + fsi_pressure = FastInterpolations.cumulative_integrate(itp_pressure) # Make a writable copy and normalize pressure integral fs_copy = copy(fs) fs_copy[:, 2] .= (fsi_pressure .- fsi_pressure[ma]) .* psio @@ -374,9 +375,9 @@ function read_chease_ascii(config::EquilibriumConfig) # Create separate interpolants for R and Z coordinates rz_in_xs = xs rz_in_R = cubic_interp((rz_in_xs, rz_in_ys), R_data; - bc=(CubicFit(), Spl.PeriodicBC()), extrap=(:extension, :wrap)) + bc=(CubicFit(), PeriodicBC()), extrap=(:extension, :wrap)) rz_in_Z = cubic_interp((rz_in_xs, rz_in_ys), Z_data; - bc=(CubicFit(), Spl.PeriodicBC()), extrap=(:extension, :wrap)) + bc=(CubicFit(), PeriodicBC()), extrap=(:extension, :wrap)) println("--> Finished reading CHEASE equilibrium.") println(" Magnetic axis at (ro=$ro, zo=$zo), psio=$psio") return InverseRunInput(config, sq_in, rz_in_xs, rz_in_ys, rz_in_R, rz_in_Z, ro, zo, psio) diff --git a/src/ForceFreeStates/ForceFreeStates.jl b/src/ForceFreeStates/ForceFreeStates.jl index 5d5adfb6..34601e70 100644 --- a/src/ForceFreeStates/ForceFreeStates.jl +++ b/src/ForceFreeStates/ForceFreeStates.jl @@ -11,7 +11,6 @@ using JLD2 using FastInterpolations import ..Equilibrium -import ..Spl import ..Utilities import ..Vacuum using Printf diff --git a/src/ForceFreeStates/Mercier.jl b/src/ForceFreeStates/Mercier.jl index cefd0a81..b23be520 100644 --- a/src/ForceFreeStates/Mercier.jl +++ b/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, ff_fs; bc=Spl.PeriodicBC()) + itp = cubic_interp(plasma_eq.rzphi_ys, ff_fs; bc=PeriodicBC()) avg = FastInterpolations.integrate(itp) # Evaluate Mercier criterion and related quantities diff --git a/src/JPEC.jl b/src/JPEC.jl index 26a94f0f..309f8779 100755 --- a/src/JPEC.jl +++ b/src/JPEC.jl @@ -1,10 +1,6 @@ # JPEC.jl module JPEC -include("Splines/Splines.jl") -import .SplinesMod as Spl -export SplinesMod, Spl - include("Utilities/Utilities.jl") import .Utilities as Utilities export Utilities @@ -95,7 +91,7 @@ function main(args::Vector{String}=String[]) else intr.debug_settings = DebugSettings() end - + # Set up variables # TODO: parallel threads logic ctrl.delta_mhigh *= 2 # for consistency with Fortran DCON TODO: why is this present in the Fortran? @@ -185,10 +181,14 @@ function main(args::Vector{String}=String[]) if ctrl.mat_flag || ctrl.ode_flag if ctrl.verbose println("Run parameters:") - println(" q0 = $(@sprintf("%.3f", equil.params.q0)), qmin = $(@sprintf("%.3f", equil.params.qmin)), qmax = $(@sprintf("%.3f", equil.params.qmax)), q95 = $(@sprintf("%.3f", equil.params.q95))") + println( + " q0 = $(@sprintf("%.3f", equil.params.q0)), qmin = $(@sprintf("%.3f", equil.params.qmin)), qmax = $(@sprintf("%.3f", equil.params.qmax)), q95 = $(@sprintf("%.3f", equil.params.q95))" + ) println(" qlim = $(@sprintf("%.5f", intr.qlim)), psilim = $(@sprintf("%.9f", intr.psilim))") println(" betat = $(@sprintf("%.3f", equil.params.betat)), betan = $(@sprintf("%.3f", equil.params.betan)), betap1 = $(@sprintf("%.3f", equil.params.betap1))") - println(" mlow = $(@sprintf("%4i", intr.mlow)), mhigh = $(@sprintf("%4i", intr.mhigh)), mpert = $(@sprintf("%4i", intr.mpert)), mband = $(@sprintf("%4i", intr.mband))") + println( + " mlow = $(@sprintf("%4i", intr.mlow)), mhigh = $(@sprintf("%4i", intr.mhigh)), mpert = $(@sprintf("%4i", intr.mpert)), mband = $(@sprintf("%4i", intr.mband))" + ) println(" nlow = $(@sprintf("%4i", intr.nlow)), nhigh = $(@sprintf("%4i", intr.nhigh)), npert = $(@sprintf("%4i", intr.npert))") end @@ -314,9 +314,8 @@ vacuum data if `vac_flag` is true. ### TODOs Combine spline unpacking if possible, too many extra lines - """ -function write_outputs_to_HDF5(ctrl::ForceFreeStatesControl, equil::Equilibrium.PlasmaEquilibrium, intr::ForceFreeStatesInternal, odet::OdeState, vac::Union{VacuumData, Nothing}) +function write_outputs_to_HDF5(ctrl::ForceFreeStatesControl, equil::Equilibrium.PlasmaEquilibrium, intr::ForceFreeStatesInternal, odet::OdeState, vac::Union{VacuumData,Nothing}) h5open(joinpath(intr.dir_path, ctrl.HDF5_filename), "w") do out_h5 diff --git a/src/PerturbedEquilibrium/PerturbedEquilibrium.jl b/src/PerturbedEquilibrium/PerturbedEquilibrium.jl index e38f1cfa..90d08b6b 100644 --- a/src/PerturbedEquilibrium/PerturbedEquilibrium.jl +++ b/src/PerturbedEquilibrium/PerturbedEquilibrium.jl @@ -7,7 +7,6 @@ using LinearAlgebra using Statistics # Import parent modules -import ..Spl import ..Equilibrium import ..ForceFreeStates import ..ForceFreeStates: OdeState, VacuumData, ForceFreeStatesInternal @@ -53,6 +52,7 @@ Computes plasma response to external forcing and calculates singular layer coupling metrics. ## Arguments + - `equil`: Equilibrium solution from Equilibrium module - `ForceFreeStates_results`: Stability calculation results from ForceFreeStates module - `vac_data`: Vacuum response data from ForceFreeStates free boundary calculation @@ -62,18 +62,20 @@ coupling metrics. - `intr`: Internal state variables ## Returns + - `PerturbedEquilibriumState`: Calculation results ## Workflow -1. Load forcing data from file -2. Compute plasma response (if enabled) -3. Calculate singular coupling metrics (if enabled) -4. Output eigenmode fields (if enabled) + + 1. Load forcing data from file + 2. Compute plasma response (if enabled) + 3. Calculate singular coupling metrics (if enabled) + 4. Output eigenmode fields (if enabled) """ function compute_perturbed_equilibrium( equil::Equilibrium.PlasmaEquilibrium, ForceFreeStates_results::OdeState, - vac_data::Union{VacuumData, Nothing}, + vac_data::Union{VacuumData,Nothing}, ffs_intr::ForceFreeStates.ForceFreeStatesInternal, ft_ctrl::ForcingTerms.ForcingTermsControl, ctrl::PerturbedEquilibriumControl, diff --git a/src/Splines/FastInterpolationsAdaptor.jl b/src/Splines/FastInterpolationsAdaptor.jl deleted file mode 100644 index 82bec8fb..00000000 --- a/src/Splines/FastInterpolationsAdaptor.jl +++ /dev/null @@ -1,223 +0,0 @@ -""" -FastInterpolationsAdaptor - Helpers and types for FastInterpolations.jl in JPEC - -This module provides: -- cumulative_integral: Exact spline integration (uses cubic spline coefficients) -- total_integral: Optimized final-value-only integration -- Re-exports of FastInterpolations boundary condition types - -## Boundary Conditions - -Use FastInterpolations' native BC types: -- `NaturalBC()`: f''(x) = 0 at endpoints -- `PeriodicBC()`: Periodic (requires fs[end] == fs[1]) -- `CubicFit()`: Automatic endpoint derivative estimation (4-point fit) - -## Performance Notes - -For monotonic access patterns (typical in ODE integration), use `hint=Ref(1)` with -`search=LinearBinary()` for optimal performance (~4 ns/pt). -""" - -using FastInterpolations -using LinearAlgebra - -# Re-export FastInterpolations BC types for convenience -const NaturalBC = FastInterpolations.NaturalBC -const PeriodicBC = FastInterpolations.PeriodicBC -const CubicFit = FastInterpolations.CubicFit -const BCPair = FastInterpolations.BCPair -const Deriv1 = FastInterpolations.Deriv1 -const Deriv2 = FastInterpolations.Deriv2 - -# ============================================================================= -# Exact Spline Integration -# ============================================================================= - -""" - cumulative_integral(xs, fs; bc=NaturalBC()) -> Vector/Matrix - -Compute exact cumulative integral of cubic spline fitted to (xs, fs). -Returns array where result[1] = 0 and result[i+1] = result[i] + ∫_{xs[i]}^{xs[i+1]} S(x) dx, -where S(x) is the cubic spline interpolant. - -For FastInterpolations' moment formulation, the exact integral over interval [x_i, x_{i+1}] is: -∫S(x)dx = h/2 * (y_i + y_{i+1}) - h³/24 * (z_i + z_{i+1}) -where h is the interval width, y are function values, and z are second derivatives (moments). - -This is more accurate than trapezoidal rule, matching the Fortran GPEC behavior. - -# Arguments - - - `xs`: Grid points (sorted ascending) - - `fs`: Function values (vector or matrix with columns as quantities) - - `bc`: Boundary condition for spline fitting (default: NaturalBC()) -""" -function cumulative_integral(xs::AbstractVector{Float64}, fs::AbstractVector{T}; bc=NaturalBC()) where {T} - npts = length(xs) - fsi = zeros(T, npts) - - # Create spline to get second derivative coefficients - itp = cubic_interp(xs, Vector{Float64}(fs); bc=bc) - y = itp.y - z = itp.z - h = itp.cache.spacing.h # interval widths: h[i] = xs[i+1] - xs[i] - - @inbounds for i in 1:(npts-1) - hi = h[i] - # Exact integral: h/2*(y_i + y_{i+1}) - h³/24*(z_i + z_{i+1}) - fsi[i+1] = fsi[i] + hi/2 * (y[i] + y[i+1]) - hi^3/24 * (z[i] + z[i+1]) - end - return fsi -end - -function cumulative_integral(xs::AbstractVector{Float64}, fs::AbstractMatrix{T}; bc=NaturalBC()) where {T} - npts, nqty = size(fs) - fsi = zeros(T, npts, nqty) - - @inbounds for k in 1:nqty - # Create spline for this column - itp = cubic_interp(xs, Vector{Float64}(fs[:, k]); bc=bc) - y = itp.y - z = itp.z - h = itp.cache.spacing.h - - for i in 1:(npts-1) - hi = h[i] - fsi[i+1, k] = fsi[i, k] + hi/2 * (y[i] + y[i+1]) - hi^3/24 * (z[i] + z[i+1]) - end - end - return fsi -end - -""" - integrate_spline(itp::CubicInterpolant) -> Float64 - -Compute exact integral of cubic spline over its entire domain. -""" -function integrate_spline(itp::FastInterpolations.CubicInterpolant{T}) where {T} - y = itp.y - z = itp.z - h = itp.cache.spacing.h - n = length(y) - - total = zero(T) - @inbounds for i in 1:(n-1) - hi = h[i] - total += hi/2 * (y[i] + y[i+1]) - hi^3/24 * (z[i] + z[i+1]) - end - return total -end - -# ============================================================================= -# Total Integral (Optimized for final value only) -# ============================================================================= - -""" - total_integral(xs, fs; bc=NaturalBC()) -> Float64 or Vector{Float64} - -Compute the total integral of a cubic spline over the entire domain [xs[1], xs[end]]. -This is an optimized version of `cumulative_integral(xs, fs)[end, :]` that avoids -allocating the full cumulative array. - -For the exact spline integration formula over interval [x_i, x_{i+1}]: -integral = h/2 * (y_i + y_{i+1}) - h^3/24 * (z_i + z_{i+1}) -where h is the interval width, y are function values, and z are second derivatives (moments). - -# Arguments - - - `xs::AbstractVector{Float64}`: Grid points (sorted ascending) - - `fs`: Function values - vector (returns scalar) or matrix (returns vector, one value per column) - - `bc`: Boundary condition for spline fitting (default: NaturalBC()) - -# Returns - - - For vector `fs`: Returns a scalar Float64 - - For matrix `fs`: Returns a Vector{Float64} of length `size(fs, 2)` - -# Performance - -Avoids allocating the full (npts,) or (npts, nqty) cumulative integral array. -For matrix input, uses CubicSeriesInterpolant internally for efficiency. - -# Example - -```julia -xs = collect(0.0:0.01:1.0) -fs = xs .^ 2 -total = total_integral(xs, fs) # Returns approx 1/3 - -# Matrix version -Y = hcat(xs .^ 2, sin.(xs)) -totals = total_integral(xs, Y) # Returns [1/3, 1-cos(1)] approximately # Create spline to get second derivative coefficients -``` -""" -function total_integral(xs::AbstractVector{Float64}, fs::AbstractVector{Float64}; bc=NaturalBC()) - # Create spline to get second derivative coefficients - itp = cubic_interp(xs, fs; bc=bc) - return integrate_spline(itp) -end - -""" - total_integral(xs, fs::AbstractMatrix; bc=NaturalBC()) -> Vector{Float64} - -Matrix version: compute total integral for each column of fs. -Returns a vector of length `size(fs, 2)`. -""" -function total_integral(xs::AbstractVector{Float64}, fs::AbstractMatrix{Float64}; bc=NaturalBC()) - npts, nqty = size(fs) - - # Use CubicSeriesInterpolant for efficient multi-quantity handling - itp = cubic_interp(xs, fs; bc=bc) - y = itp.y # (npts, nqty) - z = itp.z # (npts, nqty) - h = itp.cache.spacing.h # (npts-1,) - - # Allocate only the output vector, not the full (npts, nqty) cumulative array - result = zeros(Float64, nqty) - - @inbounds for i in 1:(npts-1) - hi = h[i] - hi_half = hi / 2 - hi3_24 = hi^3 / 24 - for k in 1:nqty - result[k] += hi_half * (y[i, k] + y[i+1, k]) - hi3_24 * (z[i, k] + z[i+1, k]) - end - end - - return result -end - -""" - total_integral!(result::AbstractVector{Float64}, xs, fs::AbstractMatrix; bc=NaturalBC()) -> result - -In-place matrix version: compute total integral for each column of fs into pre-allocated result. -Zero allocations beyond the spline creation. -""" -function total_integral!(result::AbstractVector{Float64}, xs::AbstractVector{Float64}, - fs::AbstractMatrix{Float64}; bc=NaturalBC()) - npts, nqty = size(fs) - @assert length(result) >= nqty "result vector must have at least $nqty elements" - - # Use CubicSeriesInterpolant for efficient multi-quantity handling - itp = cubic_interp(xs, fs; bc=bc) - y = itp.y # (npts, nqty) - z = itp.z # (npts, nqty) - h = itp.cache.spacing.h # (npts-1,) - - # Zero the result - @inbounds for k in 1:nqty - result[k] = 0.0 - end - - @inbounds for i in 1:(npts-1) - hi = h[i] - hi_half = hi / 2 - hi3_24 = hi^3 / 24 - for k in 1:nqty - result[k] += hi_half * (y[i, k] + y[i+1, k]) - hi3_24 * (z[i, k] + z[i+1, k]) - end - end - - return result -end diff --git a/src/Splines/Splines.jl b/src/Splines/Splines.jl deleted file mode 100644 index 5a96f8f0..00000000 --- a/src/Splines/Splines.jl +++ /dev/null @@ -1,18 +0,0 @@ -module SplinesMod - -# Pure Julia spline implementations -# FastInterpolationsAdaptor.jl provides helpers for FastInterpolations CubicInterpolant/CubicSeriesInterpolant -# The codebase uses native interpolants from FastInterpolations directly -include("FastInterpolationsAdaptor.jl") - -# Re-export FastInterpolations search strategies for use with CubicInterpolant hint/search kwargs -using FastInterpolations: LinearBinary, Binary, HintedBinary - -# Exports -export evaluate!, deriv1!, integrate!, cumulative_integral, integrate_spline, total_integral, total_integral! -export LinearBinary, Binary, HintedBinary # Search strategies for hint-based evaluation - -# Export boundary condition types -export NaturalBC, PeriodicBC, CubicFit, BCPair, Deriv1, Deriv2 - -end diff --git a/src/Vacuum/Vacuum.jl b/src/Vacuum/Vacuum.jl index a2696d88..9094a7a7 100644 --- a/src/Vacuum/Vacuum.jl +++ b/src/Vacuum/Vacuum.jl @@ -4,7 +4,6 @@ using TOML, SpecialFunctions, LinearAlgebra, Printf using FastInterpolations: cubic_interp, deriv1, PeriodicBC, NaturalBC # Import parent modules -import ..Spl import ..Equilibrium # Import FourierTransforms utility for coefficient calculation and transforms diff --git a/test/runtests_spline.jl b/test/runtests_spline.jl index b2de24d3..d35774f9 100644 --- a/test/runtests_spline.jl +++ b/test/runtests_spline.jl @@ -1,68 +1,17 @@ -@testset "Test Spline Module" begin - @info "Testing pure Julia spline implementations" - - @testset "MultiQuantityProfile - Polynomial" begin - # Make x^3 spline - cubic splines can exactly represent cubics - xs = collect(range(1.0; stop=2.0, length=21)) - fx3 = xs .^ 3 - fs = reshape(fx3, :, 1) # Make it a matrix - spline = JPEC.Spl.MultiQuantityProfile(xs, fs) - - # Test value interpolation accuracy (primary use case) - xs_fine = collect(range(1.1; stop=1.9, length=20)) - for x in xs_fine - f = JPEC.Spl.evaluate!(spline, x) - @test abs(f[1] - x^3) < 1e-4 # Value interpolation - end - - # Test first derivative (commonly used) - f1 = JPEC.Spl.deriv1!(spline, 1.5) - f = JPEC.Spl.evaluate!(spline, 1.5) - @test abs(f[1] - 1.5^3) < 1e-4 - @test abs(f1[1] - 3 * 1.5^2) < 0.1 - end - - @testset "MultiQuantityProfile - Sine" begin - # Make sine spline - xs = collect(range(Float64(pi) / 8; stop=3 * Float64(pi) / 8, length=50)) - fs = reshape(sin.(xs), :, 1) - spline = JPEC.Spl.MultiQuantityProfile(xs, fs) - - # Test value interpolation - xs_fine = collect(range(Float64(pi) / 6; stop=Float64(pi) / 3, length=20)) - for x in xs_fine - f = JPEC.Spl.evaluate!(spline, x) - @test abs(f[1] - sin(x)) < 1e-6 - end - - # Test first derivative - f1 = JPEC.Spl.deriv1!(spline, Float64(pi) / 4) - f = JPEC.Spl.evaluate!(spline, Float64(pi) / 4) - @test abs(f[1] - sin(Float64(pi) / 4)) < 1e-6 - @test abs(f1[1] - cos(Float64(pi) / 4)) < 1e-4 - end - - # NOTE: BicubicSpline has been removed - use native FastInterpolations.cubic_interp instead - # @testset "BicubicSpline" begin - # ... (tests removed) - # end -end - -# NOTE: BicubicSpline has been removed - use native FastInterpolations.cubic_interp instead -# @testset "BicubicSpline - PeriodicBC" begin -# ... (tests removed) -# end +# NOTE: Splines module has been removed - all spline functionality now uses FastInterpolations directly +# Legacy tests for MultiQuantityProfile and BicubicSpline have been removed @testset "Exact Spline Integration" begin - # Test that cumulative_integral uses exact spline integration formula - # For periodic functions with PeriodicBC, the spline can match well + # Test native FastInterpolations integration functions + using FastInterpolations: cubic_interp, cumulative_integrate, integrate, PeriodicBC - # Test with periodic function (sin) and PeriodicBC - spline fits well + # Test with periodic function (sin) and PeriodicBC xs_periodic = collect(range(0.0, 2π; length=65)) # Closed grid fs_periodic = sin.(xs_periodic) fs_periodic[end] = fs_periodic[1] # Ensure closed - fsi_periodic = JPEC.Spl.cumulative_integral(xs_periodic, fs_periodic; bc=JPEC.Spl.PeriodicBC()) + itp_periodic = cubic_interp(xs_periodic, fs_periodic; bc=PeriodicBC()) + fsi_periodic = cumulative_integrate(itp_periodic) # ∫sin(x)dx = -cos(x), so cumulative from 0 should be -cos(x) + cos(0) = 1 - cos(x) exact_periodic = 1.0 .- cos.(xs_periodic) @@ -82,26 +31,11 @@ end # Exact spline integration should be more accurate than trapezoidal @test spline_error < trap_error - # Test integrate_spline function for total integral - using FastInterpolations: cubic_interp - itp = cubic_interp(xs_periodic, fs_periodic; bc=JPEC.Spl.PeriodicBC()) - total_integral = JPEC.Spl.integrate_spline(itp) + # Test integrate function for total integral + total_integral = integrate(itp_periodic) # ∫₀^{2π} sin(x)dx = 0 @test abs(total_integral) < 1e-10 end -@testset "Empty Spline Constructors" begin - @info "Testing empty spline constructors for type stability" - - # Test empty MultiQuantityProfile - empty_mqp = JPEC.Spl.empty_MultiQuantityProfile() - @test length(empty_mqp.xs) >= 4 - @test typeof(empty_mqp) <: JPEC.Spl.MultiQuantityProfile - - # Test empty BicubicSpline - empty_bcs = JPEC.Spl.empty_BicubicSpline() - @test length(empty_bcs.xs) >= 2 - @test length(empty_bcs.ys) >= 2 - @test typeof(empty_bcs) <: JPEC.Spl.BicubicSpline -end +# NOTE: Empty spline constructor tests removed with Splines module From b51aeb91e6b0d284821d8043c4d4489cee4da3cd Mon Sep 17 00:00:00 2001 From: logan-nc <6198372+logan-nc@users.noreply.github.com> Date: Fri, 13 Feb 2026 14:30:40 -0500 Subject: [PATCH 5/8] SPLINES - IMPROVEMENT - Add LinearBinary search to 2D spline creation Add search=LinearBinary() parameter to all 2D cubic_interp calls for optimal interval search performance in hot loops: - Equilibrium.jl: flux interpolants - DirectEquilibrium.jl: psi_in, rzphi, and eqfun interpolants - InverseEquilibrium.jl: rz, rzphi, and eqfun interpolants - ReadEquilibrium.jl: psi_in and rz interpolants (CHEASE) - AnalyticEquilibrium.jl: rz and psi_in interpolants Setting the search strategy at spline creation time (instead of at each evaluation) follows the same pattern as 1D splines and ensures optimal O(log n) interval search for all 2D evaluations in loops. Co-Authored-By: Claude Sonnet 4.5 --- src/Equilibrium/AnalyticEquilibrium.jl | 6 +++--- src/Equilibrium/DirectEquilibrium.jl | 16 ++++++++-------- src/Equilibrium/Equilibrium.jl | 4 ++-- src/Equilibrium/InverseEquilibrium.jl | 18 +++++++++--------- src/Equilibrium/ReadEquilibrium.jl | 10 +++++----- 5 files changed, 27 insertions(+), 27 deletions(-) diff --git a/src/Equilibrium/AnalyticEquilibrium.jl b/src/Equilibrium/AnalyticEquilibrium.jl index daf49d7f..ef46c8bc 100644 --- a/src/Equilibrium/AnalyticEquilibrium.jl +++ b/src/Equilibrium/AnalyticEquilibrium.jl @@ -208,9 +208,9 @@ 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]; + 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]; + 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 @@ -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 diff --git a/src/Equilibrium/DirectEquilibrium.jl b/src/Equilibrium/DirectEquilibrium.jl index c678d59d..08a0b947 100644 --- a/src/Equilibrium/DirectEquilibrium.jl +++ b/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 @@ -524,13 +524,13 @@ 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]; + 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]; + 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]; + 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]; + 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` @@ -596,11 +596,11 @@ 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]; + 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]; + 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]; + 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, diff --git a/src/Equilibrium/Equilibrium.jl b/src/Equilibrium/Equilibrium.jl index da017d39..790103fd 100644 --- a/src/Equilibrium/Equilibrium.jl +++ b/src/Equilibrium/Equilibrium.jl @@ -442,9 +442,9 @@ 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]; + 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]; + 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 diff --git a/src/Equilibrium/InverseEquilibrium.jl b/src/Equilibrium/InverseEquilibrium.jl index 046be5c4..2ab55b7e 100644 --- a/src/Equilibrium/InverseEquilibrium.jl +++ b/src/Equilibrium/InverseEquilibrium.jl @@ -118,9 +118,9 @@ function equilibrium_solver(input::InverseRunInput) deta[:, end] .= deta[:, 1] # Create 2D interpolants for r² and dη - rz_rsq = cubic_interp((rz_in_xs, rz_in_ys), r2; + rz_rsq = cubic_interp((rz_in_xs, rz_in_ys), r2; search=LinearBinary(), bc=(CubicFit(), PeriodicBC()), extrap=(:extension, :wrap)) - rz_deta = cubic_interp((rz_in_xs, rz_in_ys), deta; + rz_deta = cubic_interp((rz_in_xs, rz_in_ys), deta; search=LinearBinary(), bc=(CubicFit(), PeriodicBC()), extrap=(:extension, :wrap)) # c----------------------------------------------------------------------- @@ -255,13 +255,13 @@ function equilibrium_solver(input::InverseRunInput) end qa = f_sq[mpsi+1, 4] + f1_sq[mpsi+1, 4] * (1 - sq_xs[mpsi+1]) # Create 2D interpolants for geometric quantities (rzphi) - rzphi_rsquared = cubic_interp((rzphi_xs, rzphi_ys), rzphi_fs[:, :, 1]; + rzphi_rsquared = cubic_interp((rzphi_xs, rzphi_ys), rzphi_fs[:, :, 1]; search=LinearBinary(), bc=(CubicFit(), PeriodicBC()), extrap=(:extension, :wrap)) - rzphi_offset = cubic_interp((rzphi_xs, rzphi_ys), rzphi_fs[:, :, 2]; + rzphi_offset = cubic_interp((rzphi_xs, rzphi_ys), rzphi_fs[:, :, 2]; search=LinearBinary(), bc=(CubicFit(), PeriodicBC()), extrap=(:extension, :wrap)) - rzphi_nu = cubic_interp((rzphi_xs, rzphi_ys), rzphi_fs[:, :, 3]; + rzphi_nu = cubic_interp((rzphi_xs, rzphi_ys), rzphi_fs[:, :, 3]; search=LinearBinary(), bc=(CubicFit(), PeriodicBC()), extrap=(:extension, :wrap)) - rzphi_jac = cubic_interp((rzphi_xs, rzphi_ys), rzphi_fs[:, :, 4]; + rzphi_jac = cubic_interp((rzphi_xs, rzphi_ys), rzphi_fs[:, :, 4]; search=LinearBinary(), bc=(CubicFit(), PeriodicBC()), extrap=(:extension, :wrap)) for ipsi in 0:mpsi @@ -311,11 +311,11 @@ function equilibrium_solver(input::InverseRunInput) end end # Create 2D interpolants for physics quantities (eqfun) - eqfun_B = cubic_interp((eqfun_xs, eqfun_ys), eqfun_fs[:, :, 1]; + eqfun_B = cubic_interp((eqfun_xs, eqfun_ys), eqfun_fs[:, :, 1]; search=LinearBinary(), bc=(CubicFit(), PeriodicBC()), extrap=(:extension, :wrap)) - eqfun_metric1 = cubic_interp((eqfun_xs, eqfun_ys), eqfun_fs[:, :, 2]; + eqfun_metric1 = cubic_interp((eqfun_xs, eqfun_ys), eqfun_fs[:, :, 2]; search=LinearBinary(), bc=(CubicFit(), PeriodicBC()), extrap=(:extension, :wrap)) - eqfun_metric2 = cubic_interp((eqfun_xs, eqfun_ys), eqfun_fs[:, :, 3]; + eqfun_metric2 = cubic_interp((eqfun_xs, eqfun_ys), eqfun_fs[:, :, 3]; search=LinearBinary(), bc=(CubicFit(), PeriodicBC()), extrap=(:extension, :wrap)) # Create ProfileSplines from sq interpolant diff --git a/src/Equilibrium/ReadEquilibrium.jl b/src/Equilibrium/ReadEquilibrium.jl index f3945d78..db45255f 100644 --- a/src/Equilibrium/ReadEquilibrium.jl +++ b/src/Equilibrium/ReadEquilibrium.jl @@ -115,7 +115,7 @@ function read_efit(config::EquilibriumConfig) psi_in_xs = collect(r_grid) psi_in_ys = collect(z_grid) - psi_in = cubic_interp((psi_in_xs, psi_in_ys), psi_proc; + psi_in = cubic_interp((psi_in_xs, psi_in_ys), psi_proc; search=LinearBinary(), bc=(CubicFit(), CubicFit()), extrap=(:extension, :extension)) # --- Bundle everything for the solver --- @@ -210,9 +210,9 @@ function read_chease_binary(config::EquilibriumConfig) # Create separate interpolants for R and Z coordinates rz_in_xs = xs rz_in_ys = range(0, 2π; length=mtau) |> collect - rz_in_R = cubic_interp((rz_in_xs, rz_in_ys), fs_2d[:, :, 1]; + rz_in_R = cubic_interp((rz_in_xs, rz_in_ys), fs_2d[:, :, 1]; search=LinearBinary(), bc=(CubicFit(), PeriodicBC()), extrap=(:extension, :wrap)) - rz_in_Z = cubic_interp((rz_in_xs, rz_in_ys), fs_2d[:, :, 2]; + rz_in_Z = cubic_interp((rz_in_xs, rz_in_ys), fs_2d[:, :, 2]; search=LinearBinary(), bc=(CubicFit(), PeriodicBC()), extrap=(:extension, :wrap)) println("--> Finished reading CHEASE equilibrium (Binary).") @@ -374,9 +374,9 @@ function read_chease_ascii(config::EquilibriumConfig) # Create separate interpolants for R and Z coordinates rz_in_xs = xs - rz_in_R = cubic_interp((rz_in_xs, rz_in_ys), R_data; + rz_in_R = cubic_interp((rz_in_xs, rz_in_ys), R_data; search=LinearBinary(), bc=(CubicFit(), PeriodicBC()), extrap=(:extension, :wrap)) - rz_in_Z = cubic_interp((rz_in_xs, rz_in_ys), Z_data; + rz_in_Z = cubic_interp((rz_in_xs, rz_in_ys), Z_data; search=LinearBinary(), bc=(CubicFit(), PeriodicBC()), extrap=(:extension, :wrap)) println("--> Finished reading CHEASE equilibrium.") println(" Magnetic axis at (ro=$ro, zo=$zo), psio=$psio") From 69be1ef80bb57ab8a0d7e9a4b328946a7ef7c8b6 Mon Sep 17 00:00:00 2001 From: logan-nc <6198372+logan-nc@users.noreply.github.com> Date: Fri, 13 Feb 2026 14:48:27 -0500 Subject: [PATCH 6/8] SPLINES - BUGFIX - Add missing PeriodicBC imports to modules Add PeriodicBC to FastInterpolations imports in modules that use it: - Equilibrium.jl: Added PeriodicBC to import list - ForceFreeStates.jl: Added explicit imports for cubic_interp, deriv1, PeriodicBC, LinearBinary These imports were missing after removing the Spl module wrapper, causing "PeriodicBC not defined" errors when running examples. Both DIII-D and Solovev examples now run successfully. Co-Authored-By: Claude Sonnet 4.5 --- src/Equilibrium/Equilibrium.jl | 2 +- src/ForceFreeStates/ForceFreeStates.jl | 1 + 2 files changed, 2 insertions(+), 1 deletion(-) diff --git a/src/Equilibrium/Equilibrium.jl b/src/Equilibrium/Equilibrium.jl index 790103fd..370f3858 100644 --- a/src/Equilibrium/Equilibrium.jl +++ b/src/Equilibrium/Equilibrium.jl @@ -5,7 +5,7 @@ module Equilibrium 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 --- diff --git a/src/ForceFreeStates/ForceFreeStates.jl b/src/ForceFreeStates/ForceFreeStates.jl index 34601e70..06a48fe4 100644 --- a/src/ForceFreeStates/ForceFreeStates.jl +++ b/src/ForceFreeStates/ForceFreeStates.jl @@ -9,6 +9,7 @@ using OrdinaryDiffEq using HDF5 using JLD2 using FastInterpolations +using FastInterpolations: cubic_interp, deriv1, PeriodicBC, LinearBinary import ..Equilibrium import ..Utilities From 55d13e7f461115cc9d6be43cdd7db5314fb99d6e Mon Sep 17 00:00:00 2001 From: logan-nc <6198372+logan-nc@users.noreply.github.com> Date: Fri, 13 Feb 2026 15:21:58 -0500 Subject: [PATCH 7/8] DOCS - CLEANUP - Remove splines.md documentation for removed Splines module The SplinesMod wrapper has been removed in favor of direct FastInterpolations usage. Remove the now-obsolete splines.md documentation and its reference from the documentation build configuration. Co-Authored-By: Claude Sonnet 4.5 --- docs/make.jl | 3 +- docs/src/splines.md | 128 -------------------------------------------- 2 files changed, 1 insertion(+), 130 deletions(-) delete mode 100644 docs/src/splines.md diff --git a/docs/make.jl b/docs/make.jl index c5f66fde..5e2b9ad3 100644 --- a/docs/make.jl +++ b/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 ) diff --git a/docs/src/splines.md b/docs/src/splines.md deleted file mode 100644 index aa83867b..00000000 --- a/docs/src/splines.md +++ /dev/null @@ -1,128 +0,0 @@ -# Splines Module - -The Splines module provides helpers and utilities for working with FastInterpolations.jl cubic splines, used for smooth representation of MHD equilibria. - -## Overview - -JPEC uses [FastInterpolations.jl](https://github.com/ThummeTo/FastInterpolations.jl) for high-performance cubic spline interpolation. The SplinesMod provides: -- Helper functions for exact spline integration (`cumulative_integral`, `total_integral`) -- Re-exports of FastInterpolations boundary condition types -- Re-exports of FastInterpolations search strategies for optimized evaluation - -## API Reference - -```@autodocs -Modules = [JPEC.SplinesMod] -``` - -## Boundary Condition Types - -The module re-exports FastInterpolations boundary condition types: -- `NaturalBC()`: Natural boundary condition (f''(x) = 0 at endpoints) -- `PeriodicBC()`: Periodic boundary condition (requires data[end] == data[1]) -- `CubicFit()`: Automatic endpoint derivative estimation using 4-point fit -- `BCPair(bc_left, bc_right)`: Different BCs at left and right boundaries -- `Deriv1(value)`: First derivative boundary condition -- `Deriv2(value)`: Second derivative boundary condition - -## Search Strategies - -For optimized evaluation (especially with monotonic access patterns): -- `LinearBinary()`: Linear search + binary fallback (~4 ns/point for monotonic access) -- `Binary()`: Binary search -- `HintedBinary()`: Binary search with hint - -## Integration Functions - -### cumulative_integral -```@docs -JPEC.SplinesMod.cumulative_integral -``` - -### total_integral -```@docs -JPEC.SplinesMod.total_integral -``` - -### integrate_spline -```@docs -JPEC.SplinesMod.integrate_spline -``` - -## Example Usage - -### 1D Cubic Spline -```julia -using JPEC -using FastInterpolations: cubic_interp, deriv1 - -# Create data points -xs = collect(range(0.0, 2π, length=21)) -fs = sin.(xs) - -# Create spline with natural boundary conditions -spline = cubic_interp(xs, fs; bc=JPEC.SplinesMod.NaturalBC()) - -# Evaluate at a single point -x = 1.0 -f = spline(x) - -# Evaluate derivative at a point -d1_spline = deriv1(spline) -f_deriv = d1_spline(x) - -# Evaluate at vector of new points -xs_fine = collect(range(0.0, 2π, length=100)) -fs_fine = spline.(xs_fine) - -# Compute cumulative integral -fs_integral = JPEC.SplinesMod.cumulative_integral(xs, fs; bc=JPEC.SplinesMod.NaturalBC()) -``` - -### 2D Cubic Spline (CubicInterpolantND) -```julia -using FastInterpolations: cubic_interp - -# Create 2D grid -xs = collect(range(0.0, 0.99, length=20)) -ys = collect(range(0.0, 2π, length=20)) - -# Create 2D function data -fs = zeros(20, 20) -for i in 1:20, j in 1:20 - fs[i, j] = sqrt(xs[i]) * cos(ys[j]) -end - -# Set up 2D cubic interpolant -# Use PeriodicBC on second dimension (θ direction) -spline_2d = cubic_interp((xs, ys), fs; - bc=(JPEC.SplinesMod.CubicFit(), JPEC.SplinesMod.PeriodicBC()), - extrap=(:extension, :wrap)) - -# Evaluate spline at a point -x_eval, y_eval = 0.5, π/4 -f = spline_2d((x_eval, y_eval)) - -# Evaluate with derivative (specify which dimension) -fx = spline_2d((x_eval, y_eval); deriv=(1, 0)) # ∂f/∂x -fy = spline_2d((x_eval, y_eval); deriv=(0, 1)) # ∂f/∂y -``` - -### Multi-series 1D Spline (CubicSeriesInterpolant) -```julia -using FastInterpolations: cubic_interp - -# Create multiple data series -xs = collect(range(0.0, 1.0, length=50)) -fs = hcat(sin.(2π*xs), cos.(2π*xs), xs.^2) # 3 series - -# Create series interpolant (more efficient than separate splines) -spline_series = cubic_interp(xs, fs; bc=JPEC.SplinesMod.CubicFit()) - -# Evaluate all series at once -x = 0.5 -f_all = spline_series(x) # Returns vector of length 3 - -# Compute total integral for all series -integrals = JPEC.SplinesMod.total_integral(xs, fs; bc=JPEC.SplinesMod.CubicFit()) -``` From 0c91e2d582971367a0c79afb44746d04ecaec38b Mon Sep 17 00:00:00 2001 From: Jake Halpern Date: Fri, 13 Feb 2026 17:11:42 -0500 Subject: [PATCH 8/8] TEST - IMPROVEMENT - removing unused spline unit tests --- test/runtests_spline.jl | 41 ----------------------------------------- 1 file changed, 41 deletions(-) delete mode 100644 test/runtests_spline.jl diff --git a/test/runtests_spline.jl b/test/runtests_spline.jl deleted file mode 100644 index d35774f9..00000000 --- a/test/runtests_spline.jl +++ /dev/null @@ -1,41 +0,0 @@ -# NOTE: Splines module has been removed - all spline functionality now uses FastInterpolations directly -# Legacy tests for MultiQuantityProfile and BicubicSpline have been removed - -@testset "Exact Spline Integration" begin - # Test native FastInterpolations integration functions - using FastInterpolations: cubic_interp, cumulative_integrate, integrate, PeriodicBC - - # Test with periodic function (sin) and PeriodicBC - xs_periodic = collect(range(0.0, 2π; length=65)) # Closed grid - fs_periodic = sin.(xs_periodic) - fs_periodic[end] = fs_periodic[1] # Ensure closed - - itp_periodic = cubic_interp(xs_periodic, fs_periodic; bc=PeriodicBC()) - fsi_periodic = cumulative_integrate(itp_periodic) - - # ∫sin(x)dx = -cos(x), so cumulative from 0 should be -cos(x) + cos(0) = 1 - cos(x) - exact_periodic = 1.0 .- cos.(xs_periodic) - - # Spline integration should be very accurate for smooth periodic functions - @test maximum(abs.(fsi_periodic .- exact_periodic)) < 1e-6 - - # Compare accuracy: spline vs trapezoidal for periodic case - fsi_trap = zeros(length(xs_periodic)) - for i in 1:(length(xs_periodic)-1) - h = xs_periodic[i+1] - xs_periodic[i] - fsi_trap[i+1] = fsi_trap[i] + h * (fs_periodic[i] + fs_periodic[i+1]) / 2 - end - trap_error = maximum(abs.(fsi_trap .- exact_periodic)) - spline_error = maximum(abs.(fsi_periodic .- exact_periodic)) - - # Exact spline integration should be more accurate than trapezoidal - @test spline_error < trap_error - - # Test integrate function for total integral - total_integral = integrate(itp_periodic) - - # ∫₀^{2π} sin(x)dx = 0 - @test abs(total_integral) < 1e-10 -end - -# NOTE: Empty spline constructor tests removed with Splines module