Repo - API! - Replace FourFitVars with immutable MatrixSplines split into ideal and kinetic (#383)

Author: jhalpern30

benchmarks/benchmark_delta_prime_methods.jl

@@ -41,16 +41,16 @@ function setup_and_run_solovev()
     intr.mpert = intr.mhigh - intr.mlow + 1
     intr.numpert_total = intr.mpert * intr.npert
     metric = FFS.make_metric(equil, intr.mpert)
-    ffit = FFS.make_matrix(equil, intr, metric)
-    odet, _, _, _ = FFS.riccati_eulerlagrange_integration(ctrl, equil, ffit, intr)
-    return ctrl, equil, ffit, intr, odet
+    mats = FFS.build_matrix_splines(equil, intr, metric)
+    odet, _, _, _ = FFS.riccati_eulerlagrange_integration(ctrl, equil, mats, intr)
+    return ctrl, equil, mats, intr, odet
 end
 
 println("\n=== compute_delta_prime_from_ca! consistency check ===")
 println("Verifies the standalone Δ' formula matches the inline Riccati crossing computation.")
 println("Expected error: exactly zero (same formula, same data).\n")
 
-ctrl, equil, ffit, intr, odet = setup_and_run_solovev()
+ctrl, equil, mats, intr, odet = setup_and_run_solovev()
 msing = intr.msing
 
 # Capture Δ' values set inline by riccati_cross_ideal_singular_surf! during integration

benchmarks/benchmark_riccati_der.jl

@@ -42,19 +42,19 @@ function setup_solovev()
     intr.mpert = intr.mhigh - intr.mlow + 1
     intr.numpert_total = intr.mpert * intr.npert
     metric = FFS.make_metric(equil, intr.mpert)
-    ffit = FFS.make_matrix(equil, intr, metric)
-    return ctrl, equil, ffit, intr
+    mats = FFS.build_matrix_splines(equil, intr, metric)
+    return ctrl, equil, mats, intr
 end
 
 # Evaluate the Riccati RHS explicitly from splines: dS = w†·F̄⁻¹·w - S·Ḡ·S
-function riccati_rhs_manual(S, psi, equil, ffit, intr)
+function riccati_rhs_manual(S, psi, equil, mats, intr)
     N = intr.numpert_total
     L = zeros(ComplexF64, N, N)
     Kmat = zeros(ComplexF64, N, N)
     Gmat = zeros(ComplexF64, N, N)
-    ffit.fmats_lower(vec(L), psi; hint=ffit._hint)
-    ffit.kmats(vec(Kmat), psi; hint=ffit._hint)
-    ffit.gmats(vec(Gmat), psi; hint=ffit._hint)
+    mats.ideal.F_spline_lower(vec(L), psi; hint=mats._hint)
+    mats.ideal.K_spline(vec(Kmat), psi; hint=mats._hint)
+    mats.ideal.G_spline(vec(Gmat), psi; hint=mats._hint)
 
     q = equil.profiles.q_spline(psi)
     singfac = vec(1.0 ./ ((intr.mlow:intr.mhigh) .- q .* (intr.nlow:intr.nhigh)'))
@@ -77,7 +77,7 @@ println("\n=== riccati_der! formula verification ===")
 println("Verifies riccati_der! output matches manual evaluation of Glasser 2018 Eq. 19.")
 println("Test state: Hermitian S (physical constraint). Expected error: ~machine epsilon.\n")
 
-ctrl, equil, ffit, intr = setup_solovev()
+ctrl, equil, mats, intr = setup_solovev()
 N = intr.numpert_total
 
 odet = FFS.OdeState(N, ctrl.numsteps_init, ctrl.numunorms_init, intr.msing)
@@ -101,15 +101,15 @@ max_err = let max_err = 0.0
         S = (A + A') / 2   # Hermitian by construction
 
         # Manual RHS
-        dS_manual = riccati_rhs_manual(S, psi, equil, ffit, intr)
+        dS_manual = riccati_rhs_manual(S, psi, equil, mats, intr)
 
         # riccati_der! RHS
         u_ric = zeros(ComplexF64, N, N, 2)
         du_ric = zeros(ComplexF64, N, N, 2)
         u_ric[:, :, 1] .= S
         u_ric[:, :, 2] .= Matrix{ComplexF64}(I, N, N)
         dummy_chunk = FFS.IntegrationChunk(psi, psi, false, 0, 1)
-        params = (ctrl, equil, ffit, intr, odet, dummy_chunk)
+        params = (ctrl, equil, mats, intr, odet, dummy_chunk)
         FFS.riccati_der!(du_ric, u_ric, params, psi)
         dS_ric = du_ric[:, :, 1]
 

benchmarks/benchmark_threads.jl

@@ -30,9 +30,9 @@ function run_ffs(ex; integrator)
     intr.mpert = intr.mhigh - intr.mlow + 1
     intr.numpert_total = intr.mpert * intr.npert
     metric = GeneralizedPerturbedEquilibrium.ForceFreeStates.make_metric(equil, intr.mpert)
-    ffit = GeneralizedPerturbedEquilibrium.ForceFreeStates.make_matrix(equil, intr, metric)
-    odet, _, _, _ = GeneralizedPerturbedEquilibrium.ForceFreeStates.eulerlagrange_integration(ctrl, equil, ffit, intr)
-    vac = GeneralizedPerturbedEquilibrium.ForceFreeStates.free_run(odet, ctrl, equil, ffit, intr)
+    mats = GeneralizedPerturbedEquilibrium.ForceFreeStates.build_matrix_splines(equil, intr, metric)
+    odet, _, _, _ = GeneralizedPerturbedEquilibrium.ForceFreeStates.eulerlagrange_integration(ctrl, equil, mats, intr)
+    vac = GeneralizedPerturbedEquilibrium.ForceFreeStates.free_run(odet, ctrl, equil, mats, intr)
     return real(vac.et[1]), intr.numpert_total
 end
 

docs/development/architecture.md

@@ -64,7 +64,7 @@ Splines are provided by the external `FastInterpolations` package rather than by
      - `Riccati/` - Chunked fundamental-matrix (STRIDE) driver and Δ' boundary-value problem
      - `Galerkin/` - RDCON outer-region singular Galerkin Δ' solver
      - `Matching/` - Outer↔inner resistive matching (`DeltaPrimeData`, `resonant_match_rpec`)
-     - `Fourfit.jl` - Fourier fitting routines (`FourFitVars`)
+     - `Fourfit.jl` - Fourier fitting routines (`MatrixSplines`)
      - `FixedBoundaryStability.jl` - Fixed boundary analysis
      - `Free.jl` - Free boundary stability
    - Status: Stable, core DCON functionality implemented

docs/src/stability.md

@@ -324,15 +324,15 @@ intr.mpert = intr.mhigh - intr.mlow + 1
 intr.numpert_total = intr.mpert * intr.npert
 
 metric = FFS.make_metric(equil, intr.mpert)
-ffit   = FFS.make_matrix(equil, intr, metric)
+mats   = FFS.build_matrix_splines(equil, intr, metric)
 
 # Choose integration driver.  The top-level `eulerlagrange_integration` dispatches
 # on ctrl.integrator and always returns a 4-tuple
 # (odet, propagators, chunks, S_at_surface_left).  The trailing three are `nothing`
 # on the forward path.
-odet, _, _, _ = FFS.eulerlagrange_integration(ctrl, equil, ffit, intr)
+odet, _, _, _ = FFS.eulerlagrange_integration(ctrl, equil, mats, intr)
 
-vac = FFS.free_run(odet, ctrl, equil, ffit, intr)
+vac = FFS.free_run(odet, ctrl, equil, mats, intr)
 println("Energy eigenvalue et[1] = ", real(vac.et[1]))
 ```