Repo - API! - Consolidate the v2 deck sections, driver and HDF5 writers

.pre-commit-config.yaml

@@ -77,5 +77,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|gal_flag|nstep|diagnose_ca|diagnose|ksing|power_bp|power_b|power_r|power_rc|use_galgrid|moment|eq_out|theta_out|xlmda_out|eqpsi_out|save_records|gpec_file|fnml_flag|ellip_flag|diag_flag|force_xialpha|fixed_boundary|output_eigenmodes|filter_modes|singular_point_method|compute_response|local_stability_flag|profile_file|profile_group)\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|nstep|diagnose_ca|diagnose|ksing|power_bp|power_b|power_r|power_rc|use_galgrid|moment|eq_out|theta_out|xlmda_out|eqpsi_out|save_records|gpec_file|fnml_flag|ellip_flag|diag_flag|force_xialpha|fixed_boundary|output_eigenmodes|filter_modes|singular_point_method|compute_response|local_stability_flag|gal_solver|gal_sing_order|gap_kHz_threshold|force_termination|HDF5_filename|write_outputs_to_HDF5|output_filename|kinetic_source|nn_low|nn_high|gal_[A-Za-z0-9_]+|output_benchmark_data|ion_charge|ion_mass|chi_perp_e|chi_perp_i|chi_tor|profile_file|profile_group)\s*=|^\[(DEBUG|SOL_INPUT|LAR_INPUT|TJ_ANALYTIC_INPUT|SLAYER\.scan_grid|SLAYER\.CriticalResonantField)\]|^eq_type\s*=\s*"(sol|lar|tj_analytic|tj_analytic_direct)"'
         files: ^(examples/.*\.toml|test/test_data/.*\.toml)$

benchmarks/benchmark_against_fortran_run.jl

@@ -406,9 +406,6 @@ function write_gpec_toml_coil(
         println(io, "mtheta      = $(p.mtheta_equil)")
         println(io, "flux_surface_rtol = 1e-7")
         println(io)
-        println(io, "[Wall]")
-        println(io, "shape = \"nowall\"")
-        println(io)
         println(io, "[LocalStability]")
         println(io)
         println(io, "[ForceFreeStates]")
@@ -418,13 +415,14 @@ function write_gpec_toml_coil(
         # (= outermost rational q + dmlim). psiedge from dcon.in (=1.0 → no edge dW scan).
         @printf(io, "qhigh    = %.4f\n", qhigh)
         @printf(io, "dW_edge_scan_start = %.4f\n", p.psiedge)
-        println(io, "nn_low   = $(p.nn)")
-        println(io, "nn_high  = $(p.nn)")
+        println(io, "nn       = $(p.nn)")
         println(io, "delta_mlow  = $(p.delta_mlow)")
         println(io, "delta_mhigh = $(p.delta_mhigh)")
         println(io, "mtheta_vacuum = 512")
-        println(io, "eulerlagrange_tolerance = 1e-7")
         println(io, "singfac_min = 1e-4")
+        println(io)
+        println(io, "[Riccati]")
+        println(io, "eulerlagrange_tolerance = 1e-7")
         println(io, "save_interval = $(p.euler_stride)")
         println(io)
         println(io, "[ForcingTerms]")
@@ -441,8 +439,6 @@ function write_gpec_toml_coil(
         println(io)
         println(io, "[PerturbedEquilibrium]")
         println(io, "compute_singular_coupling = true")
-        println(io, "verbose                   = true")
-        println(io, "write_outputs_to_HDF5     = true")
     end
 end
 
@@ -468,22 +464,20 @@ function write_gpec_toml_file(
         println(io, "mtheta      = $(p.mtheta_equil)")
         println(io, "flux_surface_rtol = 1e-7")
         println(io)
-        println(io, "[Wall]")
-        println(io, "shape = \"nowall\"")
-        println(io)
         println(io, "[LocalStability]")
         println(io)
         println(io, "[ForceFreeStates]")
         println(io, "vac_flag = true")
         @printf(io, "qlow     = %.4f\n", p.qlow)
         @printf(io, "dW_edge_scan_start = %.4f\n", p.psiedge)
-        println(io, "nn_low   = $(p.nn)")
-        println(io, "nn_high  = $(p.nn)")
+        println(io, "nn       = $(p.nn)")
         println(io, "delta_mlow  = $(p.delta_mlow)")
         println(io, "delta_mhigh = $(p.delta_mhigh)")
         println(io, "mtheta_vacuum = 512")
-        println(io, "eulerlagrange_tolerance = 1e-7")
         println(io, "singfac_min = 1e-4")
+        println(io)
+        println(io, "[Riccati]")
+        println(io, "eulerlagrange_tolerance = 1e-7")
         println(io, "save_interval = $(p.euler_stride)")
         println(io)
         println(io, "[ForcingTerms]")
@@ -492,8 +486,6 @@ function write_gpec_toml_file(
         println(io)
         println(io, "[PerturbedEquilibrium]")
         println(io, "compute_singular_coupling = true")
-        println(io, "verbose                   = true")
-        println(io, "write_outputs_to_HDF5     = true")
     end
 end
 
@@ -511,7 +503,7 @@ function setup_equil(p::FortranRunParams, fortran_dir::String)
         "grid_type" => p.grid_type,
         "flux_surface_rtol" => 1e-7
     )
-    eq_config = Equilibrium.EquilibriumConfig(eq_dict, fortran_dir)
+    eq_config = GeneralizedPerturbedEquilibrium.equilibrium_config(Dict{String,Any}("Equilibrium" => eq_dict), fortran_dir)
     return Equilibrium.PlasmaEquilibrium(eq_config)
 end
 

benchmarks/benchmark_coil_ForcingTerms_against_fortran.jl

@@ -206,7 +206,7 @@ t_equil = @elapsed begin
         "psi_accuracy" => 0.001,
         "flux_surface_rtol" => 1e-7,
     )
-    eq_config = Equilibrium.EquilibriumConfig(eq_dict, run_dir)
+    eq_config = GeneralizedPerturbedEquilibrium.equilibrium_config(Dict{String,Any}("Equilibrium" => eq_dict), run_dir)
     equil = Equilibrium.PlasmaEquilibrium(eq_config)
 end
 @printf "    Done in %.1f s  (bt_sign=%d)\n" t_equil equil.params.bt_sign

benchmarks/benchmark_delta_prime_methods.jl

@@ -22,15 +22,13 @@ const FFS = GeneralizedPerturbedEquilibrium.ForceFreeStates
 function setup_and_run_solovev()
     ex = joinpath(@__DIR__, "..", "test", "test_data", "regression_solovev_ideal_example")
     inputs = TOML.parsefile(joinpath(ex, "gpec.toml"))
-    inputs["ForceFreeStates"]["verbose"] = false
-    inputs["ForceFreeStates"]["integrator"] = "riccati"
-    intr = FFS.ForceFreeStatesInternal(; dir_path=ex)
-    ctrl = FFS.ForceFreeStatesControl(;
-        (Symbol(k) => v for (k, v) in inputs["ForceFreeStates"])...)
-    eq_config = GeneralizedPerturbedEquilibrium.Equilibrium.EquilibriumConfig(inputs["Equilibrium"], ex)
+    inputs["ForceFreeStates"]["integrator"] = "Riccati"
+    inputs["Riccati"] = pop!(inputs, "Forward")
+    intr = FFS.ForceFreeStatesInternal()
+    ctrl = GeneralizedPerturbedEquilibrium.force_balance_control(inputs, ex)
+    eq_config = GeneralizedPerturbedEquilibrium.equilibrium_config(inputs, ex)
     equil = GeneralizedPerturbedEquilibrium.Equilibrium.PlasmaEquilibrium(eq_config)
-    intr.wall_settings = GeneralizedPerturbedEquilibrium.Vacuum.WallShapeSettings(;
-        (Symbol(k) => v for (k, v) in inputs["Wall"])...)
+    intr.wall_settings = GeneralizedPerturbedEquilibrium.Vacuum.parse_wall(get(inputs, "Wall", nothing))
     FFS.sing_lim!(intr, ctrl, equil)
     intr.nlow = ctrl.nn_low;
     intr.nhigh = ctrl.nn_high;

benchmarks/benchmark_diiid_ideal_ntv_torque.jl

@@ -125,31 +125,25 @@ mtheta = 256
 newq0 = 0
 flux_surface_rtol = 1e-7
 
-[Wall]
-shape = "nowall"
-
 [LocalStability]
 
 [ForceFreeStates]
 vac_flag = true
-force_termination = false
 
 dW_edge_scan_start = 1.0      # No edge-scan truncation (rational_offset_fraction truncation off)
 qlow = 1.02
 qhigh = 1e3
 sing_start = 0
 
-nn_low = 1
-nn_high = 1
+nn = 1
 delta_mlow = 8
 delta_mhigh = 8
 mtheta_vacuum = 512
+singfac_min = 1e-4
 
-kinetic_source = "fixed"
-kinetic_factor = 0.0
+[Riccati]
 eulerlagrange_tolerance = 1e-7
 save_interval = 3
-singfac_min = 1e-4
 renorm_threshold = 1e4
 """
         )

benchmarks/benchmark_diiid_kinetic_stability.jl

@@ -5,7 +5,7 @@
 Benchmark the self-consistent kinetic-DCON path (KF → FFS → PE) against
 Fortran GPEC's kinetic DCON reference.
 
-Runs `GPE.main()` with `kinetic_source="calculated"` and `kinetic_factor=1.0`
+Runs `GPE.main()` with `kinetic = "KineticForces"`
 against the EFIT g-file and `.kin` profile taken from a Fortran GPEC
 kinetic example directory, then compares the least-stable total-energy
 eigenvalue `ForceFreeStates/FreeBoundaryStability/eigenmode_energies[1]` against `W_t_eigenvalue[:, 0]` in the Fortran
@@ -123,8 +123,6 @@ function build_benchmark_tomldir(eq_file::String, kin_file::String)
 [Equilibrium]
 eq_filename = "$eq_name"
 kinetic_file = "$kin_file"
-eq_type = "efit"
-jac_type = "hamada"
 grid_type = "pow1"             # Matches Fortran DIIID_kinetic_example (equil/grid.f90 powspace)
 psilow = 0.01
 psihigh = 0.993
@@ -133,39 +131,31 @@ mtheta = 256
 newq0 = 0
 flux_surface_rtol = 1e-7
 
-[Wall]
-shape = "nowall"
-
 [LocalStability]
 
 [ForceFreeStates]
 vac_flag = true
-force_termination = true       # Skip PE+KF post-processing — we only need FFS eigenvalues
 
 dW_edge_scan_start = 1.0      # No edge-scan truncation (rational_offset_fraction truncation off)
 qlow = 1.02
 qhigh = 1e3
 sing_start = 0
 
-nn_low = 1
-nn_high = 1
+nn = 1
 delta_mlow = 8
 delta_mhigh = 8
 mtheta_vacuum = 512
-thmax0 = 1
 
-kinetic_source = "calculated"
-kinetic_factor = 1.0
-eulerlagrange_tolerance = 1e-7
+integrator = "Forward"
+kinetic = "KineticForces"
 singfac_min = 1e-4
+
+[Forward]
+eulerlagrange_tolerance = 1e-7
 renorm_threshold = 1e4
-write_outputs_to_HDF5 = true
 
 [KineticForces]
-nn = 1
-nl = 4
-zi = 1
-mi = 2
+bounce_harmonic_max = 4
 nutype = "harmonic"
 f0type = "maxwellian"
 """
@@ -192,7 +182,7 @@ function run_benchmark(fortran_dir::String=default_fortran_dir())
 
     tomldir = build_benchmark_tomldir(inputs.eq_file, inputs.kin_file)
 
-    _p("\n--- Running GPE.main with kinetic_source=\"calculated\" ---")
+    _p("\n--- Running GPE.main with kinetic = \"KineticForces\" ---")
     t0 = time()
     GPE.main([tomldir])
     wall = time() - t0

benchmarks/benchmark_q_vs_iota_edge.jl

@@ -22,7 +22,7 @@ const EXAMPLE_DIR = joinpath(@__DIR__, "..", "examples", "DIIID-like_ideal_examp
 function reference_q()
     inputs, _, additional_input = GPE.build_inputs_from_toml(EXAMPLE_DIR)
     equil_dict = merge(inputs["Equilibrium"], Dict{String,Any}("grid_type" => "core_edge_packed", "mpsi" => 1024))
-    eq_config = GPE.Equilibrium.EquilibriumConfig(equil_dict, EXAMPLE_DIR)
+    eq_config = GPE.equilibrium_config(Dict{String,Any}("Equilibrium" => equil_dict), EXAMPLE_DIR)
     equil = GPE.Equilibrium.PlasmaEquilibrium(eq_config, additional_input)
     # The dict rides along so the mpsi scan can rebuild configs; EquilibriumConfig is immutable.
     return equil, equil_dict
@@ -65,7 +65,7 @@ function main()
 
     for N in (16, 32, 64, 128, 256)
         # Same log_asymptotic knot layout the auto grid uses for fixed mpsi
-        cfg = GPE.Equilibrium.EquilibriumConfig(merge(equil_dict, Dict{String,Any}("mpsi" => N)), EXAMPLE_DIR)
+        cfg = GPE.equilibrium_config(Dict{String,Any}("Equilibrium" => merge(equil_dict, Dict{String,Any}("mpsi" => N))), EXAMPLE_DIR)
         knots = GPE.Equilibrium._build_psi_grid(cfg, psilow, psihigh)
         q_nodes = [q_ref_itp(p) for p in knots]
 

benchmarks/benchmark_riccati_der.jl

@@ -24,14 +24,11 @@ const FFS = GeneralizedPerturbedEquilibrium.ForceFreeStates
 function setup_solovev()
     ex = joinpath(@__DIR__, "..", "test", "test_data", "regression_solovev_ideal_example")
     inputs = TOML.parsefile(joinpath(ex, "gpec.toml"))
-    inputs["ForceFreeStates"]["verbose"] = false
-    intr = FFS.ForceFreeStatesInternal(; dir_path=ex)
-    ctrl = FFS.ForceFreeStatesControl(;
-        (Symbol(k) => v for (k, v) in inputs["ForceFreeStates"])...)
-    eq_config = GeneralizedPerturbedEquilibrium.Equilibrium.EquilibriumConfig(inputs["Equilibrium"], ex)
+    intr = FFS.ForceFreeStatesInternal()
+    ctrl = GeneralizedPerturbedEquilibrium.force_balance_control(inputs, ex)
+    eq_config = GeneralizedPerturbedEquilibrium.equilibrium_config(inputs, ex)
     equil = GeneralizedPerturbedEquilibrium.Equilibrium.PlasmaEquilibrium(eq_config)
-    intr.wall_settings = GeneralizedPerturbedEquilibrium.Vacuum.WallShapeSettings(;
-        (Symbol(k) => v for (k, v) in inputs["Wall"])...)
+    intr.wall_settings = GeneralizedPerturbedEquilibrium.Vacuum.parse_wall(get(inputs, "Wall", nothing))
     FFS.sing_lim!(intr, ctrl, equil)
     intr.nlow = ctrl.nn_low;
     intr.nhigh = ctrl.nn_high;
@@ -80,7 +77,7 @@ println("Test state: Hermitian S (physical constraint). Expected error: ~machine
 ctrl, equil, mats, intr = setup_solovev()
 N = intr.numpert_total
 
-odet = FFS.OdeState(N, ctrl.numsteps_init, ctrl.numunorms_init, intr.msing)
+odet = FFS.OdeState(N, ctrl.algorithm.numsteps_init, ctrl.algorithm.numunorms_init, intr.msing)
 FFS.initialize_el_at_axis!(odet, ctrl, equil.profiles, intr)
 chunks = FFS.chunk_el_integration_bounds(odet, ctrl, intr)
 

benchmarks/benchmark_solovev_kinetic_stability.jl

@@ -2,7 +2,7 @@
 """
     benchmark_solovev_kinetic_stability.jl
 
-Benchmark the self-consistent kinetic-DCON path (`kinetic_source="calculated"`)
+Benchmark the self-consistent kinetic-DCON path (`kinetic = "KineticForces"`)
 for the analytic **Solovev** equilibrium against Fortran GPEC's kinetic DCON.
 
 Unlike the DIIID benchmark (which reads a pre-existing `dcon_output_n1.nc`),
@@ -103,7 +103,7 @@ function build_matched_fortran_deck(workdir::String)
 
     # dcon.in — kinetic flags already match; align mode band, edge truncation, n.
     dconin = joinpath(workdir, "dcon.in")
-    replace_namelist_value!(dconin, "nn",          string(ffs["nn_low"]))
+    replace_namelist_value!(dconin, "nn",          string(first(ffs["nn"])))
     replace_namelist_value!(dconin, "delta_mlow",  string(get(ffs, "delta_mlow", 0)))
     replace_namelist_value!(dconin, "delta_mhigh", string(get(ffs, "delta_mhigh", 0)))
     replace_namelist_value!(dconin, "qlow",        string(ffs["qlow"]))
@@ -171,7 +171,7 @@ function run_benchmark()
     fort_re, fort_im = run_fortran_reference(workdir)
     _p(@sprintf("  Fortran W_t[1]: Re = %.6f  Im = %.6f", fort_re, fort_im))
 
-    _p("\n--- Running GPE.main with kinetic_source=\"calculated\" ---")
+    _p("\n--- Running GPE.main with kinetic = \"KineticForces\" ---")
     jul_re, jul_im, wall = run_julia_reference()
     _p(@sprintf("  GPE.main completed in %.1f s", wall))
 

benchmarks/benchmark_threads.jl

@@ -1,5 +1,5 @@
 # Thread-scaling benchmark for the Riccati chunked propagator integration.
-# Runs the Solovev (N=8) and DIIID-like (N=26) examples with integrator="riccati"
+# Runs the Solovev (N=8) and DIIID-like (N=26) examples with the Riccati integrator
 # across 1, 2, 4, 8 threads and compares against the forward path. The chunk list is
 # thread-independent, so et[1] must not move with the thread count — only wall-clock.
 #
@@ -8,32 +8,13 @@
 
 using GeneralizedPerturbedEquilibrium, TOML, Printf, Statistics
 
+# One solve of the example deck with `integrator` (Forward or Riccati), sharing the deck's ODE numerics.
 function run_ffs(ex; integrator)
-    inputs = TOML.parsefile(joinpath(ex, "gpec.toml"))
-    inputs["ForceFreeStates"]["verbose"] = false
-    inputs["ForceFreeStates"]["integrator"] = integrator
-    inputs["ForceFreeStates"]["write_outputs_to_HDF5"] = false
-    intr = GeneralizedPerturbedEquilibrium.ForceFreeStates.ForceFreeStatesInternal(; dir_path=ex)
-    ctrl = GeneralizedPerturbedEquilibrium.ForceFreeStates.ForceFreeStatesControl(;
-        (Symbol(k) => v for (k, v) in inputs["ForceFreeStates"])...)
-    eq_config = GeneralizedPerturbedEquilibrium.Equilibrium.EquilibriumConfig(inputs["Equilibrium"], ex)
-    equil = GeneralizedPerturbedEquilibrium.Equilibrium.PlasmaEquilibrium(eq_config)
-    intr.wall_settings = GeneralizedPerturbedEquilibrium.Vacuum.WallShapeSettings(;
-        (Symbol(k) => v for (k, v) in inputs["Wall"])...)
-    GeneralizedPerturbedEquilibrium.ForceFreeStates.sing_lim!(intr, ctrl, equil)
-    intr.nlow = ctrl.nn_low;
-    intr.nhigh = ctrl.nn_high;
-    intr.npert = 1
-    GeneralizedPerturbedEquilibrium.ForceFreeStates.sing_find!(intr, equil)
-    intr.mlow = min(intr.nlow * equil.params.qmin, 0) - 4 - ctrl.delta_mlow
-    intr.mhigh = trunc(Int, intr.nhigh * equil.params.qmax) + ctrl.delta_mhigh
-    intr.mpert = intr.mhigh - intr.mlow + 1
-    intr.numpert_total = intr.mpert * intr.npert
-    metric = GeneralizedPerturbedEquilibrium.ForceFreeStates.make_metric(equil, intr.mpert)
-    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
+    GPE = GeneralizedPerturbedEquilibrium
+    inputs, config, analytic = GPE.build_inputs_from_toml(ex)
+    alg = integrator(; GPE.section_keywords(inputs["Forward"], "Forward", GPE._struct_fields(Forward))...)
+    ffs = solve(PlasmaEquilibrium(config, analytic), alg; GPE.force_balance_keywords(inputs, ex)...)
+    return real(ffs.free_boundary.et[1]), ffs.numpert_total
 end
 
 function timed_run(ex; integrator, nwarm=1, nrep=2)
@@ -60,8 +41,8 @@ diiid_ex = joinpath(root, "examples", "DIIID-like_ideal_example")
 println("\n=== Thread-scaling benchmark ($(nthreads) thread(s)) ===\n")
 
 for (label, ex) in [("Solovev", sol_ex), ("DIIID-like", diiid_ex)]
-    t_fwd, et_fwd, N = timed_run(ex; integrator="forward")
-    t_ric, et_ric, _ = timed_run(ex; integrator="riccati")
+    t_fwd, et_fwd, N = timed_run(ex; integrator=Forward)
+    t_ric, et_ric, _ = timed_run(ex; integrator=Riccati)
 
     err_ric = abs(et_ric - et_fwd) / abs(et_fwd) * 100
 

benchmarks/equil_method_comparison.jl

@@ -38,7 +38,7 @@ function make_config(path::String, eq_type::String; psihigh=nothing)
     raw["Equilibrium"]["eq_type"] = eq_type
     psihigh !== nothing && (raw["Equilibrium"]["psihigh"] = psihigh)
     base = dirname(path)
-    return Equilibrium.EquilibriumConfig(raw["Equilibrium"], base)
+    return GeneralizedPerturbedEquilibrium.equilibrium_config(raw, base)
 end
 
 # ─── Run each method (3 warm runs, report average of last 2) ─────────────────