OpenFUSIONToolkit · github-actions · Jul 30, 2026 · Jul 30, 2026 · Jul 30, 2026 · Jul 30, 2026
diff --git a/benchmarks/benchmark_against_fortran_run.jl b/benchmarks/benchmark_against_fortran_run.jl
diff --git a/benchmarks/benchmark_diiid_ideal_ntv_torque.jl b/benchmarks/benchmark_diiid_ideal_ntv_torque.jl
@@ -129,8 +129,6 @@ shape = "nowall"
 
 [ForceFreeStates]
 local_stability_flag = true
-mat_flag = true
-ode_flag = true
 vac_flag = true
 force_termination = false
 

diff --git a/benchmarks/benchmark_diiid_kinetic_stability.jl b/benchmarks/benchmark_diiid_kinetic_stability.jl
@@ -135,8 +135,6 @@ shape = "nowall"
 
 [ForceFreeStates]
 local_stability_flag = true
-mat_flag = true
-ode_flag = true
 vac_flag = true
 force_termination = true       # Skip PE+KF post-processing — we only need FFS eigenvalues
 

diff --git a/benchmarks/gal_validation/tokamaker_beta/gal_tok_worker.jl b/benchmarks/gal_validation/tokamaker_beta/gal_tok_worker.jl
@@ -24,7 +24,7 @@ function base_inputs(dir)
             "etol" => 1e-10, "force_termination" => false),
         "Wall" => Dict("shape" => "nowall", "a" => 0.2415),
         "ForceFreeStates" => Dict(
-            "local_stability_flag" => true, "mat_flag" => true, "ode_flag" => true,
+            "local_stability_flag" => true,
             "vac_flag" => true,
             "qlow" => 1.02, "qhigh" => 1e3, "sing_start" => 0,
             "nn_low" => 1, "nn_high" => 1, "delta_mlow" => 8, "delta_mhigh" => 8,
@@ -49,9 +49,13 @@ function run_pipeline(psihigh)
     eq_config = EQ.EquilibriumConfig(inputs["Equilibrium"], dir)
     equil = EQ.setup_equilibrium(eq_config, nothing)
     FFS.sing_lim!(intr, ctrl, equil)
-    xs = equil.profiles.xs; ls = zeros(length(xs), 5); FFS.compute_ballooning_stability!(ctrl, ls, equil)
+    xs = equil.profiles.xs;
+    ls = zeros(length(xs), 5);
+    FFS.compute_ballooning_stability!(ctrl, ls, equil)
     intr.locstab = cubic_interp(xs, Series(ls); extrap=ExtendExtrap())
-    intr.nlow = ctrl.nn_low; intr.nhigh = ctrl.nn_high; intr.npert = 1
+    intr.nlow = ctrl.nn_low;
+    intr.nhigh = ctrl.nn_high;
+    intr.npert = 1
     FFS.sing_find!(intr, equil)
     # Replicate main()'s post-find surface filter: drop surfaces outside [qlow, qlim].
     # Without this the worker keeps the q=1 internal-kink surface (and mis-aligns the STRIDE
@@ -81,13 +85,19 @@ end
 P = nothing
 for ph in (0.999, 0.99, 0.98, 0.97, 0.95)
     try
-        global P = run_pipeline(ph); break
+        global P = run_pipeline(ph);
+        break
     catch e
         @warn "pipeline failed at psihigh=$ph for $(basename(GEQDSK)); retrying lower" exception=e
     end
 end
 P === nothing && error("all psihigh retries failed for $(basename(GEQDSK))")
-ctrl = P.ctrl; equil = P.equil; ffit = P.ffit; intr = P.intr; vac_data = P.vac_data; dir = P.dir
+ctrl = P.ctrl;
+equil = P.equil;
+ffit = P.ffit;
+intr = P.intr;
+vac_data = P.vac_data;
+dir = P.dir
 @printf("  using psihigh=%.3f for %s\n", P.psihigh, basename(GEQDSK))
 dpm = intr.delta_prime_matrix
 m_s = [s.m[1] for s in intr.sing if s.psifac < intr.psilim]
@@ -100,7 +110,8 @@ stride_diag = [i <= ndpm ? real(dpm[i, i]) : NaN for i in 1:length(m_s)]
 open(OUTCSV, "a") do io
     for pfac in PFACS
         ctrl.gal_pfac = pfac
-        gdiag = fill(NaN, length(m_s)); gm = Int[]
+        gdiag = fill(NaN, length(m_s));
+        gm = Int[]
         try
             res = FFS.galerkin_solve(ctrl, equil, ffit, intr; vac_data=vac_data)
             gm = res.sing_m

diff --git a/examples/DIIID-like_gal_resistive_example/gpec.toml b/examples/DIIID-like_gal_resistive_example/gpec.toml
@@ -24,8 +24,6 @@ equal_arc_wall = true          # Equal arc length distribution of nodes on wall
 
 [ForceFreeStates]
 local_stability_flag = true    # Perform local stability analysis (Mercier and ballooning) across the ψ profile
-mat_flag = true                # Construct coefficient matrices for diagnostic purposes (required by gal_flag)
-ode_flag = true                # Integrate ODE's for determining stability of internal long-wavelength mode (must be true for GPEC)
 vac_flag = true                # Compute plasma, vacuum, and total energies for free-boundary modes
 
 psiedge = 0.99                 # Edge dW scan band: dW(ψ) computed for ψ ∈ [psiedge, psilim], integration truncated at peak
@@ -55,7 +53,7 @@ dmlim               = 0.2      # Truncate integration at (last_rational_q + dmli
 
 # Outer-region singular Galerkin Δ′ solver with rpec coil columns.
 # This case exercises the gal solve AND the rpec coil/edge path (delta_coil) for regression tracking.
-gal_flag               = true  # Enable the outer-region singular Galerkin Δ′ solve (requires mat_flag = true)
+gal_flag               = true  # Enable the outer-region singular Galerkin Δ′ solve
 gal_solver             = "LU"  # Banded solver: "LU" (zgbtrf/zgbtrs) or "cholesky" (zpbtrf/zpbtrs). rpec requires "LU".
 gal_nx                 = 256   # Hermite-cubic elements per interval between singular surfaces
 gal_nq                 = 6     # Gauss-Lobatto quadrature order per element

diff --git a/examples/DIIID-like_gal_resistive_pe_example/gpec.toml b/examples/DIIID-like_gal_resistive_pe_example/gpec.toml
@@ -24,8 +24,6 @@ equal_arc_wall = true          # Equal arc length distribution of nodes on wall
 
 [ForceFreeStates]
 local_stability_flag = true    # Perform local stability analysis (Mercier and ballooning) across the ψ profile
-mat_flag = true                # Construct coefficient matrices for diagnostic purposes (required by gal_flag)
-ode_flag = true                # Integrate ODE's for determining stability of internal long-wavelength mode (must be true for GPEC)
 vac_flag = true                # Compute plasma, vacuum, and total energies for free-boundary modes
 
 psiedge = 0.99                 # Edge dW scan band: dW(ψ) computed for ψ ∈ [psiedge, psilim], integration truncated at peak
@@ -55,7 +53,7 @@ dmlim               = 0.2      # Truncate integration at (last_rational_q + dmli
 
 # Outer-region singular Galerkin Δ′ solver with rpec coil columns.
 # This case exercises the gal solve AND the rpec coil/edge path (delta_coil) for regression tracking.
-gal_flag               = true  # Enable the outer-region singular Galerkin Δ′ solve (requires mat_flag = true)
+gal_flag               = true  # Enable the outer-region singular Galerkin Δ′ solve
 gal_solver             = "LU"  # Banded solver: "LU" (zgbtrf/zgbtrs) or "cholesky" (zpbtrf/zpbtrs). rpec requires "LU".
 gal_nx                 = 256   # Hermite-cubic elements per interval between singular surfaces
 gal_nq                 = 6     # Gauss-Lobatto quadrature order per element

diff --git a/examples/DIIID-like_ideal_example/gpec.toml b/examples/DIIID-like_ideal_example/gpec.toml
@@ -29,8 +29,6 @@ equal_arc_wall = true          # Equal arc length distribution of nodes on wall
 
 [ForceFreeStates]
 local_stability_flag = true    # Perform local stability analysis (Mercier and ballooning) across the ψ profile
-mat_flag = true                # Construct coefficient matrices for diagnostic purposes
-ode_flag = true                # Integrate ODEs for stability of the internal long-wavelength mode (must be true for GPEC)
 vac_flag = true                # Compute plasma, vacuum, and total energies for free-boundary modes
 
 psiedge = 0.99                 # Edge dW(ψ) diagnostic scan band [psiedge, psilim]; set ≥ psilim to disable

diff --git a/examples/DIIID-like_ideal_example_IMAS/gpec.toml b/examples/DIIID-like_ideal_example_IMAS/gpec.toml
@@ -24,8 +24,6 @@ shape = "nowall"               # Wall shape (nowall, conformal, elliptical, dee,
 nn_low               = 1       # Smallest toroidal mode number to include
 nn_high              = 1       # Largest toroidal mode number to include
 vac_flag             = true    # Compute plasma, vacuum, and total energies for free-boundary modes
-ode_flag             = true    # Integrate ODEs for stability of the internal long-wavelength mode (must be true for GPEC)
-mat_flag             = true    # Construct coefficient matrices for diagnostic purposes
 local_stability_flag = false   # Perform local stability analysis (Mercier and ballooning) across the ψ profile
 write_outputs_to_HDF5 = false  # Write outputs to the HDF5 file
 verbose              = false   # Enable verbose logging

diff --git a/examples/DIIID-like_ideal_example_IMAS/run_imas_example.jl b/examples/DIIID-like_ideal_example_IMAS/run_imas_example.jl
@@ -26,41 +26,42 @@ const GEQDSK = joinpath(@__DIR__, "../DIIID-like_ideal_example/TkMkr_D3Dlike_Hmo
 # Only the [Equilibrium] section differs between IMAS and g-file paths.
 const SHARED_CONFIG = Dict(
     "ForceFreeStates" => Dict(
-        "nn_low"               => 1,
-        "nn_high"              => 1,
-        "vac_flag"             => true,
-        "ode_flag"             => true,
-        "mat_flag"             => true,
+        "nn_low" => 1,
+        "nn_high" => 1,
+        "vac_flag" => true,
         "local_stability_flag" => false,
-        "write_outputs_to_HDF5"=> false,
-        "verbose"              => false,
-        "delta_mhigh"          => 16,
-        "qlow"                 => 1.02,
+        "write_outputs_to_HDF5" => false,
+        "verbose" => false,
+        "delta_mhigh" => 16,
+        "qlow" => 1.02
     ),
-    "Wall" => Dict("shape" => "nowall"),
+    "Wall" => Dict("shape" => "nowall")
 )
 
 const EQ_COMMON = Dict(
-    "jac_type"   => "boozer",
-    "psilow"     => 0.01,
-    "psihigh"    => 0.9995,
-    "mpsi"       => 128,
-    "mtheta"     => 256,
+    "jac_type" => "boozer",
+    "psilow" => 0.01,
+    "psihigh" => 0.9995,
+    "mpsi" => 128,
+    "mtheta" => 256
 )
 
 @info "Loading g-file and converting to IMAS (COCOS 11)..."
-g  = EFIT.readg(GEQDSK)
+g = EFIT.readg(GEQDSK)
 dd = IMASdd.dd()
 EFIT.geqdsk2imas!([g], dd; dd_cocos=11)
 
 @info "Running GPEC from IMAS dd..."
-tmpdir_imas = mktempdir(prefix="gpec_imas_")
+tmpdir_imas = mktempdir(; prefix="gpec_imas_")
 try
     config_imas = merge(SHARED_CONFIG, Dict("Equilibrium" =>
         merge(EQ_COMMON, Dict("eq_type" => "imas", "eq_filename" => "from_dd", "imas_cocos" => 11))))
-    open(joinpath(tmpdir_imas, "gpec.toml"), "w") do io; TOML.print(io, config_imas); end
+    open(joinpath(tmpdir_imas, "gpec.toml"), "w") do io
+        ;
+        TOML.print(io, config_imas);
+    end
     result_imas = GPEC.main([tmpdir_imas]; dd=dd)
-    global et_imas  = real(result_imas.vac_data.et[1])
+    global et_imas = real(result_imas.vac_data.et[1])
     global mpert_imas = result_imas.intr.mpert
     GPEC.write_imas(dd, result_imas)
     @assert dd.mhd_linear.time_slice[1].toroidal_mode[1].energy_perturbed ≈ et_imas
@@ -70,13 +71,16 @@ end
 @info @sprintf "IMAS path:   δW(n=1) = %.6f  [mpert = %d]" et_imas mpert_imas
 
 @info "Running GPEC from g-file (reference)..."
-tmpdir_gfile = mktempdir(prefix="gpec_gfile_")
+tmpdir_gfile = mktempdir(; prefix="gpec_gfile_")
 try
     config_gfile = merge(SHARED_CONFIG, Dict("Equilibrium" =>
         merge(EQ_COMMON, Dict("eq_type" => "efit", "eq_filename" => GEQDSK))))
-    open(joinpath(tmpdir_gfile, "gpec.toml"), "w") do io; TOML.print(io, config_gfile); end
+    open(joinpath(tmpdir_gfile, "gpec.toml"), "w") do io
+        ;
+        TOML.print(io, config_gfile);
+    end
     result_gfile = GPEC.main([tmpdir_gfile])
-    global et_gfile  = real(result_gfile.vac_data.et[1])
+    global et_gfile = real(result_gfile.vac_data.et[1])
     global mpert_gfile = result_gfile.intr.mpert
 finally
     rm(tmpdir_gfile; recursive=true)
@@ -89,7 +93,7 @@ println("\n" * "="^60)
 println("  IMAS vs g-file comparison (n=1, nowall, boozer)")
 println("="^60)
 @printf "  G-file δW(n=1) : %+.6f  [mpert = %d]\n" et_gfile mpert_gfile
-@printf "  IMAS   δW(n=1) : %+.6f  [mpert = %d]\n" et_imas  mpert_imas
+@printf "  IMAS   δW(n=1) : %+.6f  [mpert = %d]\n" et_imas mpert_imas
 @printf "  Relative diff  : %.2f%%\n" rel_diff * 100
 println("="^60)
 

diff --git a/examples/LAR_beta_scan/gpec.toml b/examples/LAR_beta_scan/gpec.toml
@@ -35,8 +35,6 @@ a     = 20                     # Distance from plasma (conformal) or shape param
 
 [ForceFreeStates]
 local_stability_flag = true    # Perform local stability analysis (Mercier and ballooning) across the ψ profile
-mat_flag = true                # Construct coefficient matrices for diagnostic purposes
-ode_flag = true                # Integrate ODEs for stability of the internal long-wavelength mode (must be true for GPEC)
 vac_flag = true                # Compute plasma, vacuum, and total energies for free-boundary modes
 
 qlow         = 1.02            # Integration initiated at q determined by min(q0, qlow)

diff --git a/examples/LAR_epsilon_scan/gpec.toml b/examples/LAR_epsilon_scan/gpec.toml
@@ -36,8 +36,6 @@ a     = 20                     # Distance from plasma (conformal) or shape param
 
 [ForceFreeStates]
 local_stability_flag = true    # Perform local stability analysis (Mercier and ballooning) across the ψ profile
-mat_flag = true                # Construct coefficient matrices for diagnostic purposes
-ode_flag = true                # Integrate ODEs for stability of the internal long-wavelength mode (must be true for GPEC)
 vac_flag = true                # Compute plasma, vacuum, and total energies for free-boundary modes
 
 qlow         = 1.02            # Integration initiated at q determined by min(q0, qlow)

diff --git a/examples/Solovev_ideal_example/gpec.toml b/examples/Solovev_ideal_example/gpec.toml
@@ -46,8 +46,6 @@ reg_spot = 0.05                         # Regularization width for singular surf
 
 [ForceFreeStates]
 local_stability_flag = true   # Perform local stability analysis (Mercier and ballooning) across the ψ profile
-mat_flag = true               # Construct coefficient matrices for diagnostic purposes
-ode_flag = true               # Integrate ODEs for stability of the internal long-wavelength mode (must be true for GPEC)
 vac_flag = true               # Compute plasma, vacuum, and total energies for free-boundary modes
 
 psiedge = 0.99                # Edge dW(ψ) diagnostic scan band [psiedge, psilim]; set ≥ psilim to disable

diff --git a/examples/Solovev_ideal_example_3D/gpec.toml b/examples/Solovev_ideal_example_3D/gpec.toml
@@ -16,8 +16,6 @@ force_termination = false               # Terminate after equilibrium setup (ski
 
 [ForceFreeStates]
 local_stability_flag = true   # Perform local stability analysis (Mercier and ballooning) across the ψ profile
-mat_flag = true               # Construct coefficient matrices for diagnostic purposes
-ode_flag = true               # Integrate ODEs for stability of the internal long-wavelength mode (must be true for GPEC)
 vac_flag = true               # Compute plasma, vacuum, and total energies for free-boundary modes
 
 psiedge = 1.0                # Edge dW(ψ) diagnostic scan band [psiedge, psilim]; set ≥ psilim to disable

diff --git a/examples/Solovev_ideal_example_multi_n/gpec.toml b/examples/Solovev_ideal_example_multi_n/gpec.toml
@@ -28,8 +28,6 @@ equal_arc_wall = true                   # Equal arc length distribution of nodes
 
 [ForceFreeStates]
 local_stability_flag = true   # Perform local stability analysis (Mercier and ballooning) across the ψ profile
-mat_flag = true               # Construct coefficient matrices for diagnostic purposes
-ode_flag = true               # Integrate ODEs for stability of the internal long-wavelength mode (must be true for GPEC)
 vac_flag = true               # Compute plasma, vacuum, and total energies for free-boundary modes
 
 psiedge = 0.99                # Edge dW(ψ) diagnostic scan band [psiedge, psilim]; set ≥ psilim to disable

diff --git a/examples/Solovev_ideal_example_multi_n/single_n_1/gpec.toml b/examples/Solovev_ideal_example_multi_n/single_n_1/gpec.toml
@@ -27,8 +27,6 @@ equal_arc_wall = true                   # Equal arc length distribution of nodes
 
 [ForceFreeStates]
 local_stability_flag = true   # Perform local stability analysis (Mercier and ballooning) across the ψ profile
-mat_flag = true               # Construct coefficient matrices for diagnostic purposes
-ode_flag = true               # Integrate ODEs for stability of the internal long-wavelength mode (must be true for GPEC)
 vac_flag = true               # Compute plasma, vacuum, and total energies for free-boundary modes
 
 psiedge = 0.99                # Edge dW(ψ) diagnostic scan band [psiedge, psilim]; set ≥ psilim to disable

diff --git a/examples/Solovev_ideal_example_multi_n/single_n_2/gpec.toml b/examples/Solovev_ideal_example_multi_n/single_n_2/gpec.toml
@@ -27,8 +27,6 @@ equal_arc_wall = true                   # Equal arc length distribution of nodes
 
 [ForceFreeStates]
 local_stability_flag = true   # Perform local stability analysis (Mercier and ballooning) across the ψ profile
-mat_flag = true               # Construct coefficient matrices for diagnostic purposes
-ode_flag = true               # Integrate ODEs for stability of the internal long-wavelength mode (must be true for GPEC)
 vac_flag = true               # Compute plasma, vacuum, and total energies for free-boundary modes
 
 psiedge = 0.99                # Edge dW(ψ) diagnostic scan band [psiedge, psilim]; set ≥ psilim to disable

diff --git a/examples/Solovev_kinetic_NTV_example/gpec.toml b/examples/Solovev_kinetic_NTV_example/gpec.toml
@@ -48,8 +48,6 @@ reg_spot = 0.05                         # Regularization width for singular surf
 
 [ForceFreeStates]
 local_stability_flag = true   # Perform local stability analysis (Mercier and ballooning) across the ψ profile
-mat_flag = true               # Construct coefficient matrices for diagnostic purposes
-ode_flag = true               # Integrate ODEs for stability of the internal long-wavelength mode (must be true for GPEC)
 vac_flag = true               # Compute plasma, vacuum, and total energies for free-boundary modes
 
 psiedge = 0.99                # Edge dW(ψ) diagnostic scan band [psiedge, psilim]; set ≥ psilim to disable

diff --git a/examples/Solovev_kinetic_calculated_example/gpec.toml b/examples/Solovev_kinetic_calculated_example/gpec.toml
@@ -27,8 +27,6 @@ equal_arc_wall = true                   # Equal arc length distribution of nodes
 
 [ForceFreeStates]
 local_stability_flag = true   # Perform local stability analysis (Mercier and ballooning) across the ψ profile
-mat_flag = false              # Construct coefficient matrices for diagnostic purposes
-ode_flag = true               # Integrate ODEs for stability of the internal long-wavelength mode (must be true for GPEC)
 vac_flag = true               # Compute plasma, vacuum, and total energies for free-boundary modes
 
 psiedge = 0.99                # Edge dW scan band: dW(ψ) computed for ψ ∈ [psiedge, psilim], integration truncated at peak

diff --git a/examples/a10_kinetic_example/gpec.toml b/examples/a10_kinetic_example/gpec.toml
@@ -23,8 +23,6 @@ shape = "nowall"               # Wall shape (nowall, conformal, elliptical, dee,
 
 [ForceFreeStates]
 local_stability_flag = true    # Perform local stability analysis (Mercier and ballooning) across the ψ profile
-mat_flag = true                # Construct coefficient matrices for diagnostic purposes
-ode_flag = true                # Integrate ODEs for stability of the internal long-wavelength mode (must be true for GPEC)
 vac_flag = true                # Compute plasma, vacuum, and total energies for free-boundary modes
 
 set_psilim_via_dmlim = false   # FALSE for limited/analytical equilibria — rationals sparse, dmlim would chop too much edge

diff --git a/src/Analysis/ForceFreeStates.jl b/src/Analysis/ForceFreeStates.jl
@@ -41,7 +41,7 @@ function plot_mode_displacement(h5path; modes=1:5, save_path=nothing)
     mhigh = mlow + mpert - 1
     dW = isempty(et) ? nothing : et[1]
     title_str = isnothing(dW) ? "Least stable mode" :
-        "Least stable mode, δW = $(round(real(dW), sigdigits=4))"
+                "Least stable mode, δW = $(round(real(dW), sigdigits=4))"
 
     p = plot(;
         xlims=(0, 1),
@@ -425,7 +425,7 @@ end
 
 Plot `cond(F̄)` vs ψ from the kinetic-singular-surface scan stored in
 `singular/kinetic/` (populated when ForceFreeStates runs with
-`kinetic_factor > 0`, `ode_flag = true`, `singfac_min > 0`).
+`kinetic_factor > 0`, `singfac_min > 0`).
 
 `F̄` is the kinetic Euler-Lagrange matrix formed by Schur-reducing the six
 kinetic matrices against the ideal A/B/C/D/E/H blocks (Logan 2015 Appendix
@@ -471,7 +471,7 @@ function plot_cond_fbar(h5path; save_path=nothing, zoom=false)
     end
 
     if isempty(scan_psi)
-        @warn "No kinetic-singular-surface scan in $h5path — rerun with kinetic_factor>0, ode_flag=true, singfac_min>0"
+        @warn "No kinetic-singular-surface scan in $h5path — rerun with kinetic_factor>0, singfac_min>0"
         return nothing
     end