FFS - REFACTOR - Publish solves as ForceFreeStatesResult; staged main; solve API (refactor plan interface PR) (#393)

Author: matt-pharr

.pre-commit-config.yaml

@@ -69,5 +69,5 @@ repos:
       - id: toml-no-deprecated-keys
         name: 'TOML conventions: no deprecated config keys'
         language: pygrep
-        entry: '^(mer_flag|force_wv_symmetry|ode_flag|cyl_flag|mat_flag|reform_eq_with_psilim|use_riccati|use_parallel|parallel_threads|populate_dense_xi|power_bp|power_b|power_r|power_rc)\s*='
+        entry: '^(mer_flag|force_wv_symmetry|ode_flag|cyl_flag|mat_flag|reform_eq_with_psilim|use_riccati|use_parallel|parallel_threads|populate_dense_xi|gal_flag|power_bp|power_b|power_r|power_rc)\s*='
         files: ^(examples/.*\.toml|test/test_data/.*\.toml)$

Project.toml

@@ -6,6 +6,7 @@ version = "0.1.0"
 
 [deps]
 AdaptiveArrayPools = "4f381ef7-9af0-4cbe-99d4-cf36d7b0f233"
+CommonSolve = "38540f10-b2f7-11e9-35d8-d573e4eb0ff2"
 Contour = "d38c429a-6771-53c6-b99e-75d170b6e991"
 Dates = "ade2ca70-3891-5945-98fb-dc099432e06a"
 DelaunayTriangulation = "927a84f5-c5f4-47a5-9785-b46e178433df"
@@ -38,6 +39,7 @@ Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
 
 [compat]
 AdaptiveArrayPools = "0.3.5"
+CommonSolve = "0.2"
 Contour = "0.6.3"
 DelaunayTriangulation = "1.6.6"
 DelimitedFiles = "1.9.1"

REFACTOR_PLAN.md

@@ -31,7 +31,6 @@ using Plots
 # Load GPEC
 using GeneralizedPerturbedEquilibrium
 const GPE = GeneralizedPerturbedEquilibrium
-const FFS = GPE.ForceFreeStates
 const KF = GPE.KineticForces
 const Eq = GPE.Equilibrium
 const PE = GPE.PerturbedEquilibrium
@@ -237,15 +236,12 @@ function run_benchmark(fortran_dir::String=default_fortran_dir())
     _p("\n--- Equilibrium + ForceFreeStates (via main()) ---")
     t0 = time()
     result = GPE.main([tomldir])
-    equil = result.equil
-    intr = result.intr
-    ctrl = result.ctrl
+    ffs = result.ffs
+    equil = ffs.equil
+    metric = ffs.metric
 
     _pf("  FFS completed in %.1f s\n", time() - t0)
-    _pf("  mpert=%d, mlow=%d, mhigh=%d\n", intr.mpert, intr.mlow, intr.mhigh)
-
-    # Build metric (needed for JBB deweighting)
-    metric = FFS.make_metric(equil, intr.mpert)
+    _pf("  mpert=%d, mlow=%d, mhigh=%d\n", ffs.mpert, ffs.mlow, ffs.mhigh)
 
     # Load Fortran xclebsch data
     _p("\n--- Load Fortran xclebsch ---")
@@ -257,8 +253,8 @@ function run_benchmark(fortran_dir::String=default_fortran_dir())
         npsi_f, mpert_f, mlow_f)
     _pf("  ψ range: [%.6f, %.6f]\n", psi_grid_f[1], psi_grid_f[end])
 
-    if mpert_f != intr.mpert || mlow_f != intr.mlow
-        @warn "Mode ranges differ: Fortran mpert=$mpert_f,mlow=$mlow_f vs Julia mpert=$(intr.mpert),mlow=$(intr.mlow)"
+    if mpert_f != ffs.mpert || mlow_f != ffs.mlow
+        @warn "Mode ranges differ: Fortran mpert=$mpert_f,mlow=$mlow_f vs Julia mpert=$(ffs.mpert),mlow=$(ffs.mlow)"
     end
 
     # Build PE state from Fortran data and run JBB deweighting
@@ -291,7 +287,7 @@ function run_benchmark(fortran_dir::String=default_fortran_dir())
     kf_intr = KF.KineticForcesInternal(equil; verbose=false)
 
     # Run set_perturbation_data! — builds dbob_m, divx_m, xs_m via JBB deweighting
-    KF.set_perturbation_data!(kf_intr, pe_state, intr, equil, metric)
+    KF.set_perturbation_data!(kf_intr, pe_state, ffs, equil, metric)
 
     _pf("  JBB deweighting completed in %.1f s\n", time() - t1)
 

benchmarks/benchmark_diiid_ideal_ntv_torque.jl

@@ -18,13 +18,13 @@ ksel = length(ARGS) >= 3 ? ARGS[3] : "highest"
 
 to_c(a) = eltype(a) <: Complex ? ComplexF64.(a) : map(x -> ComplexF64(x.re, x.im), a)
 
-et, wt, u1, psiE, gxi, psiG, issing, mlow, sing_psi = h5open(h5path) do f
+et, wt, u1, psiE, gxi, psiG, mlow, sing_psi = h5open(h5path) do f
     (to_c(read(f["ForceFreeStates/FreeBoundaryStability/eigenmode_energies"])),
         to_c(read(f["ForceFreeStates/FreeBoundaryStability/W_freeboundary_eigenmodes"])),
         to_c(read(f["ForceFreeStates/Solutions/ForwardIntegration/xi_psi"])), read(f["ForceFreeStates/Solutions/ForwardIntegration/psi"]),
-        to_c(read(f["ForceFreeStates/Solutions/GalerkinIntegration/Match/xi"])), read(f["ForceFreeStates/Solutions/GalerkinIntegration/Solution/psi"]),
-        Bool.(read(f["ForceFreeStates/Solutions/GalerkinIntegration/Solution/is_rational"])), read(f["Info/mlow"]),
-        read(f["SingularSurfaces/GalerkinDeltaPrime/rational_psi"]))
+        to_c(read(f["ForceFreeStates/Solutions/GalerkinIntegration/xi_psi"])), read(f["ForceFreeStates/Solutions/GalerkinIntegration/psi"]),
+        read(f["Info/mlow"]),
+        read(f["ForceFreeStates/Solutions/GalerkinIntegration/rational_psi"]))
 end
 
 mpert = size(u1, 1)
@@ -38,10 +38,10 @@ w = wt[:, k]
 cEL = u1[:, :, end] \ w
 xiE = reduce(hcat, (u1[:, :, ip] * cEL for ip in 1:size(u1, 3)))   # (mpert, nE)
 
-# gal ideal profile (identity-at-edge ⇒ coefficient = w); drop on-surface points
-keep = .!issing
-psiGk = psiG[keep]
-xiG = reduce(hcat, (gxi[:, ip, :] * w for ip in findall(keep)))    # (mpert, nGk)
+# gal ideal profile (identity-at-edge ⇒ coefficient = w); the closed-profile grid
+# already excludes on-surface points
+psiGk = psiG
+xiG = reduce(hcat, (gxi[:, :, ip] * w for ip in eachindex(psiG)))    # (mpert, nGk)
 
 ms = mlow .+ (0:mpert-1)
 peak = [maximum(abs, @view xiE[i, :]) for i in 1:mpert]

benchmarks/compare_gal_vs_el.jl

@@ -3,7 +3,7 @@
 #   (1) IDEAL galerkin matched ξ   (gal_match_flag=true, gal_ideal_flag=true)
 #   (2) FORWARD ξ                  (gal_match_flag=false)
 #
-# PE writes no ψ grid, so it's reconstructed: gal-ideal → ForceFreeStates/Solutions/GalerkinIntegration/Solution/psi minus issing points;
+# PE writes no ψ grid, so it's reconstructed: gal-ideal → ForceFreeStates/Solutions/GalerkinIntegration/psi (already excludes on-surface points);
 # forward → ForceFreeStates/Solutions/ForwardIntegration/psi.
 # Usage: julia --project=. benchmarks/compare_jbgradpsi_m2.jl [gal_h5] [shoot_h5] [out.png] [m]
 
@@ -19,9 +19,8 @@ to_c(a) = eltype(a) <: Complex ? ComplexF64.(a) : map(x -> ComplexF64(x.re, x.im
 # gal-ideal run: PE grid = gal solution grid with the on-surface (issing) points dropped
 pa_g, psi_g, mlow, sing_psi, sing_m = h5open(gal_h5) do f
     pa = to_c(read(f["PerturbedEquilibrium/Response/psi_area"]))   # [npsi, mpert]
-    iss = Bool.(read(f["ForceFreeStates/Solutions/GalerkinIntegration/Solution/is_rational"]))
-    (pa, read(f["ForceFreeStates/Solutions/GalerkinIntegration/Solution/psi"])[.!iss], read(f["Info/mlow"]),
-        read(f["SingularSurfaces/GalerkinDeltaPrime/rational_psi"]), read(f["SingularSurfaces/GalerkinDeltaPrime/rational_m"]))
+    (pa, read(f["ForceFreeStates/Solutions/GalerkinIntegration/psi"]), read(f["Info/mlow"]),
+        read(f["ForceFreeStates/Solutions/GalerkinIntegration/rational_psi"]), read(f["ForceFreeStates/Solutions/GalerkinIntegration/rational_m"]))
 end
 # forward run: PE grid = ForceFreeStates/Solutions/ForwardIntegration/psi
 pa_s, psi_s = h5open(sh_h5) do f

benchmarks/compare_jbgradpsi_m2.jl

@@ -22,7 +22,7 @@ et, wt, u1, psi, mlow, sing_psi = h5open(h5path) do f
         to_c(read(f["ForceFreeStates/Solutions/ForwardIntegration/xi_psi"])),
         read(f["ForceFreeStates/Solutions/ForwardIntegration/psi"]),
         read(f["Info/mlow"]),
-        haskey(f, "SingularSurfaces/GalerkinDeltaPrime/rational_psi") ? read(f["SingularSurfaces/GalerkinDeltaPrime/rational_psi"]) : Float64[])
+        haskey(f, "ForceFreeStates/Solutions/GalerkinIntegration/rational_psi") ? read(f["ForceFreeStates/Solutions/GalerkinIntegration/rational_psi"]) : Float64[])
 end
 
 mpert, _, nstep = size(u1)

benchmarks/plot_xi_eigenmode.jl

@@ -16,7 +16,7 @@ function read_m2(h5; gal::Bool)
     h5open(h5) do f
         pa = to_c(read(f["PerturbedEquilibrium/Response/psi_area"]))   # [npsi, mpert]
         col = mtarget - read(f["Info/mlow"]) + 1
-        psi = gal ? read(f["ForceFreeStates/Solutions/GalerkinIntegration/Solution/psi"])[.!Bool.(read(f["ForceFreeStates/Solutions/GalerkinIntegration/Solution/is_rational"]))] :
+        psi = gal ? read(f["ForceFreeStates/Solutions/GalerkinIntegration/psi"]) :
               read(f["ForceFreeStates/Solutions/ForwardIntegration/psi"])
         (psi, pa[:, col])
     end
@@ -34,7 +34,7 @@ eta_ref = 8e-8
 
 # rational surface for m=target
 sing_psi, sing_m = h5open(joinpath(scandirs[1], "gpec.h5")) do f
-    (read(f["SingularSurfaces/GalerkinDeltaPrime/rational_psi"]), read(f["SingularSurfaces/GalerkinDeltaPrime/rational_m"]))
+    (read(f["ForceFreeStates/Solutions/GalerkinIntegration/rational_psi"]), read(f["ForceFreeStates/Solutions/GalerkinIntegration/rational_m"]))
 end
 psi_res = mtarget in sing_m ? sing_psi[findfirst(==(mtarget), sing_m)] : NaN
 

benchmarks/scan_resistivity_m2.jl

@@ -15,7 +15,7 @@ function read_m2(h5; gal::Bool)
     h5open(h5) do f
         pa = to_c(read(f["PerturbedEquilibrium/Response/psi_area"]))
         col = mtarget - read(f["Info/mlow"]) + 1
-        psi = gal ? read(f["ForceFreeStates/Solutions/GalerkinIntegration/Solution/psi"])[.!Bool.(read(f["ForceFreeStates/Solutions/GalerkinIntegration/Solution/is_rational"]))] :
+        psi = gal ? read(f["ForceFreeStates/Solutions/GalerkinIntegration/psi"]) :
               read(f["ForceFreeStates/Solutions/ForwardIntegration/psi"])
         (psi, pa[:, col])
     end
@@ -29,7 +29,7 @@ scandirs, rots = scandirs[ord], rots[ord]
 @printf("%d scan runs: rotation f = %s Hz  (η fixed = 8e-8)\n", length(rots), join((@sprintf("%g", r) for r in rots), ", "))
 
 sing_psi, sing_m = h5open(joinpath(scandirs[1], "gpec.h5")) do f
-    (read(f["SingularSurfaces/GalerkinDeltaPrime/rational_psi"]), read(f["SingularSurfaces/GalerkinDeltaPrime/rational_m"]))
+    (read(f["ForceFreeStates/Solutions/GalerkinIntegration/rational_psi"]), read(f["ForceFreeStates/Solutions/GalerkinIntegration/rational_m"]))
 end
 psi_res = mtarget in sing_m ? sing_psi[findfirst(==(mtarget), sing_m)] : NaN
 

benchmarks/scan_rotation_m2.jl

@@ -2,27 +2,30 @@
 #   1. cout / cin / deltar all zero (no resistive plasma combination, no inner layer)
 #   2. matched ξ_j == the gal coil column sols(:,:,2·msing+j)   (and ξ′ likewise)
 # Usage: julia --project=. benchmarks/verify_gal_ideal.jl [path/to/gpec.h5]
+# Requires [DEBUG] gal_basis_output = true in the deck (compares against the raw-basis dump).
 using HDF5, Printf, LinearAlgebra
 
 h5path = length(ARGS) >= 1 ? ARGS[1] : "/tmp/gal_ideal_test/gpec.h5"
 to_c(a) = eltype(a) <: Complex ? ComplexF64.(a) : map(x -> ComplexF64(x.re, x.im), a)
 
-cout, deltar, mxi, mdxi, sols, sols_d, sing_psi = h5open(h5path) do f
+cout, deltar, mxi, mdxi, sols, sols_d, iss, sing_psi = h5open(h5path) do f
     (to_c(read(f["ForceFreeStates/Solutions/GalerkinIntegration/Match/cout"])), to_c(read(f["ForceFreeStates/Solutions/GalerkinIntegration/Match/Delta_r"])),
-        to_c(read(f["ForceFreeStates/Solutions/GalerkinIntegration/Match/xi"])), to_c(read(f["ForceFreeStates/Solutions/GalerkinIntegration/Match/dxidpsi"])),
-        to_c(read(f["ForceFreeStates/Solutions/GalerkinIntegration/Solution/xi_psi"])), to_c(read(f["ForceFreeStates/Solutions/GalerkinIntegration/Solution/dxi_psidpsi"])),
-        read(f["SingularSurfaces/GalerkinDeltaPrime/rational_psi"]))
+        to_c(read(f["ForceFreeStates/Solutions/GalerkinIntegration/xi_psi"])), to_c(read(f["ForceFreeStates/Solutions/GalerkinIntegration/dxi_psidpsi"])),
+        to_c(read(f["ForceFreeStates/Solutions/GalerkinIntegration/Basis/xi_psi"])), to_c(read(f["ForceFreeStates/Solutions/GalerkinIntegration/Basis/dxi_psidpsi"])),
+        Bool.(read(f["ForceFreeStates/Solutions/GalerkinIntegration/Basis/is_rational"])),
+        read(f["ForceFreeStates/Solutions/GalerkinIntegration/rational_psi"]))
 end
 msing = length(sing_psi)
-mpert, ngrid, mcoil = size(mxi)
+mpert, mcoil, ngrid = size(mxi)
+keep = .!iss
 
 @printf("[1] ‖cout‖ = %.2e, ‖deltar‖ = %.2e   %s\n", norm(cout), norm(deltar),
     (norm(cout) == 0 && norm(deltar) == 0) ? "✓ no resistive combination (ideal)" : "✗")
 
 # matched ξ_j should equal the coil column sols(:,:,2msing+j)
 err = maximum(1:mcoil) do j
     csol = 2msing + j
-    max(norm(mxi[:, :, j] - sols[:, :, csol]), norm(mdxi[:, :, j] - sols_d[:, :, csol]))
+    max(norm(mxi[:, j, :] - sols[:, csol, keep]), norm(mdxi[:, j, :] - sols_d[:, csol, keep]))
 end
 @printf("[2] max‖ξ_matched − coil column‖ = %.2e   %s\n", err,
     err < 1e-12 ? "✓ matched ξ == bare ideal coil column" : "✗")