DOCS/CONFIG - REFACTOR - Rename TOML config variables to descriptive spellings (closes #287)

.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|power_bp|power_b|power_r|power_rc)\s*='
+        entry: '^(mer_flag|force_wv_symmetry|ode_flag|cyl_flag|mat_flag|power_bp|power_b|power_r|power_rc|nstep|diagnose_ca|use_parallel|use_riccati|parallel_threads|psiedge|newq0|use_galgrid|reg_spot|dc_type|dr_val|dgeo_val|xnom|ynom|znom)\s*='
         files: ^(examples/.*\.toml|test/test_data/.*\.toml)$

CLAUDE.md

@@ -146,6 +146,10 @@ Additional file hygiene (enforced by pre-commit hooks):
 - Files must end with exactly one newline
 - LF line endings only (no CRLF)
 
+### HDF5 Output Conventions
+
+The `gpec.h5` schema follows one physics-first convention (CamelCase groups at all levels, snake_case datasets, data-driven tokens verbatim, inputs only under `Input/`, five named top-level physics-topic exceptions). **Do not invent new group names or echo inputs into output groups** — read **[`docs/development/hdf5-conventions.md`](docs/development/hdf5-conventions.md)** before adding or moving any HDF5 output, and keep the regression-harness extractor's legacy map in sync on renames.
+
 ### TOML Annotation Conventions
 
 Config-style TOML files (`examples/*/gpec.toml`, `examples/*/sol.toml`, `test/test_data/*` fixtures, `regression-harness/cases/*.toml`) follow one shared annotation style (header comment block, inline `# description` on every variable line sourced from the matching config struct's docstring, no Fortran references, no deprecated variables). **Do not invent a new convention** — read **[`docs/development/toml-conventions.md`](docs/development/toml-conventions.md)** in full before adding or editing one of these files.

Project.toml

@@ -7,6 +7,7 @@ version = "0.1.0"
 [deps]
 AdaptiveArrayPools = "4f381ef7-9af0-4cbe-99d4-cf36d7b0f233"
 Contour = "d38c429a-6771-53c6-b99e-75d170b6e991"
+Dates = "ade2ca70-3891-5945-98fb-dc099432e06a"
 DelaunayTriangulation = "927a84f5-c5f4-47a5-9785-b46e178433df"
 DelimitedFiles = "8bb1440f-4735-579b-a4ab-409b98df4dab"
 DiffEqCallbacks = "459566f4-90b8-5000-8ac3-15dfb0a30def"

benchmarks/benchmark_against_fortran_run.jl

@@ -149,6 +149,7 @@ function parse_fortran_run(dir::String)
     eq_file = _find_string(eq_text, "eq_filename"; default="")
     jac_type = _find_string(eq_text, "jac_type"; default="hamada")
     grid_type = _find_string(eq_text, "grid_type"; default="ldp")
+    grid_type == "ldp" && (grid_type = "rational_packed")  # map the legacy Fortran value to the Julia spelling
     psilow = _find_scalar(eq_text, "psilow"; default=1e-4)
     psihigh = _find_scalar(eq_text, "psihigh"; default=0.993)
     mpsi = _find_int(eq_text, "mpsi"; default=128)
@@ -312,17 +313,17 @@ function load_julia_outputs(h5_path::String)
 
     julia = Dict{String,Any}()
     h5open(h5_path, "r") do f
-        julia["psilim"] = read(f, "info/psilim")
-        julia["qlim"] = read(f, "info/qlim")
-        julia["et"] = read(f, "vacuum/et")
-        julia["psi_q"] = read(f, "splines/profiles/xs")
-        julia["q"] = read(f, "splines/profiles/q")
-        julia["di"] = haskey(f, "locstab/di") ? read(f, "locstab/di") : Float64[]
-        julia["dr"] = haskey(f, "locstab/dr") ? read(f, "locstab/dr") : Float64[]
+        julia["psilim"] = read(f, "Info/psilim")
+        julia["qlim"] = read(f, "Info/qlim")
+        julia["et"] = read(f, "ForceFreeStates/FreeBoundaryStability/eigenmode_energies")
+        julia["psi_q"] = read(f, "Equilibrium/Profiles/xs")
+        julia["q"] = read(f, "Equilibrium/Profiles/q")
+        julia["di"] = haskey(f, "LocalStability/di") ? read(f, "LocalStability/di") : Float64[]
+        julia["dr"] = haskey(f, "LocalStability/dr") ? read(f, "LocalStability/dr") : Float64[]
 
-        julia["psio"] = haskey(f, "equil/psio") ? read(f, "equil/psio") : NaN
+        julia["psio"] = haskey(f, "Equilibrium/psio") ? read(f, "Equilibrium/psio") : NaN
 
-        sc = "perturbed_equilibrium/singular_coupling"
+        sc = "PerturbedEquilibrium/SingularCoupling"
         julia["rational_psi"] = haskey(f, "$sc/rational_psi") ? read(f, "$sc/rational_psi") : Float64[]
         julia["rational_q"] = haskey(f, "$sc/rational_q") ? read(f, "$sc/rational_q") : Float64[]
         julia["rational_n"] = haskey(f, "$sc/rational_n") ? read(f, "$sc/rational_n") : Int[]
@@ -333,30 +334,31 @@ function load_julia_outputs(h5_path::String)
         julia["chirikov_parameter"] = haskey(f, "$sc/chirikov_parameter") ? read(f, "$sc/chirikov_parameter") : Float64[]
         julia["delta_prime"] = haskey(f, "$sc/delta_prime") ? read(f, "$sc/delta_prime") : ComplexF64[]
 
-        pe = "perturbed_equilibrium"
+        pe = "PerturbedEquilibrium"
         # Fortran Phi_x/Phi_tot are the area-weighted field b̄ (tesla), matching forcing/response_b_area directly.
         julia["forcing_vec"] = haskey(f, "$pe/forcing_b_area") ? read(f, "$pe/forcing_b_area") : ComplexF64[]   # Phi_x   = b̄_x
         julia["response_vec"] = haskey(f, "$pe/response_b_area") ? read(f, "$pe/response_b_area") : ComplexF64[]   # Phi_tot = b̄_tot
-        julia["b_n"] = haskey(f, "$pe/response/b_n") ? read(f, "$pe/response/b_n") : Matrix{ComplexF64}(undef, 0, 0)
-        julia["Jbgradpsi"] = haskey(f, "$pe/response/b_psi_area_weighted") ? read(f, "$pe/response/b_psi_area_weighted") : Matrix{ComplexF64}(undef, 0, 0)
-        julia["xi_psi"] = haskey(f, "$pe/response/xi_psi") ? read(f, "$pe/response/xi_psi") : Matrix{ComplexF64}(undef, 0, 0)
-        julia["xi_n"] = haskey(f, "$pe/response/xi_n") ? read(f, "$pe/response/xi_n") : Matrix{ComplexF64}(undef, 0, 0)
-        julia["clebsch_psi1"] = haskey(f, "$pe/response/clebsch_psi1") ? read(f, "$pe/response/clebsch_psi1") : Matrix{ComplexF64}(undef, 0, 0)
-        julia["clebsch_alpha"] = haskey(f, "$pe/response/clebsch_alpha") ? read(f, "$pe/response/clebsch_alpha") : Matrix{ComplexF64}(undef, 0, 0)
-        julia["psi_grid"] = haskey(f, "integration/psi") ? read(f, "integration/psi") : Float64[]
+        julia["b_n"] = haskey(f, "$pe/Response/b_n") ? read(f, "$pe/Response/b_n") : Matrix{ComplexF64}(undef, 0, 0)
+        julia["Jbgradpsi"] = haskey(f, "$pe/Response/b_psi_area_weighted") ? read(f, "$pe/Response/b_psi_area_weighted") : Matrix{ComplexF64}(undef, 0, 0)
+        julia["xi_psi"] = haskey(f, "$pe/Response/xi_psi") ? read(f, "$pe/Response/xi_psi") : Matrix{ComplexF64}(undef, 0, 0)
+        julia["xi_n"] = haskey(f, "$pe/Response/xi_n") ? read(f, "$pe/Response/xi_n") : Matrix{ComplexF64}(undef, 0, 0)
+        julia["clebsch_psi1"] = haskey(f, "$pe/Response/clebsch_psi1") ? read(f, "$pe/Response/clebsch_psi1") : Matrix{ComplexF64}(undef, 0, 0)
+        julia["clebsch_alpha"] = haskey(f, "$pe/Response/clebsch_alpha") ? read(f, "$pe/Response/clebsch_alpha") : Matrix{ComplexF64}(undef, 0, 0)
+        julia["psi_grid"] = haskey(f, "ForceFreeStates/Solutions/ForwardIntegration/psi") ? read(f, "ForceFreeStates/Solutions/ForwardIntegration/psi") : Float64[]
 
         # R,Z,φ: loaded via modes_to_theta helper below (not raw modes)
         julia["h5_path"] = h5_path  # stash for modes_to_theta
         # mn_index[:, 1] = m values, mn_index[:, 2] = n values for each mode index
-        julia["m_modes"] = haskey(f, "info/mn_index") ? Int.(read(f, "info/mn_index")[:, 1]) : Int[]
+        julia["m_modes"] = haskey(f, "Info/mn_index") ? Int.(read(f, "Info/mn_index")[:, 1]) : Int[]
 
         # Control surface matrices
-        rm = "$pe/response_matrices"
+        rm = "$pe/ResponseMatrices"
         julia["permeability"] = haskey(f, "$rm/permeability") ? read(f, "$rm/permeability") : Matrix{ComplexF64}(undef, 0, 0)
         julia["plasma_inductance"] = haskey(f, "$rm/plasma_inductance") ? read(f, "$rm/plasma_inductance") : Matrix{ComplexF64}(undef, 0, 0)
         julia["surface_inductance"] = haskey(f, "$rm/surface_inductance") ? read(f, "$rm/surface_inductance") : Matrix{ComplexF64}(undef, 0, 0)
         julia["reluctance"] = haskey(f, "$rm/reluctance") ? read(f, "$rm/reluctance") : Matrix{ComplexF64}(undef, 0, 0)
-        julia["wt0"] = haskey(f, "vacuum/wt0") ? read(f, "vacuum/wt0") : Matrix{ComplexF64}(undef, 0, 0)
+        fbs = "ForceFreeStates/FreeBoundaryStability"
+        julia["wt0"] = haskey(f, "$fbs/W_freeboundary") ? read(f, "$fbs/W_freeboundary") : Matrix{ComplexF64}(undef, 0, 0)
     end
     return julia
 end
@@ -415,7 +417,7 @@ function write_gpec_toml_coil(
         # Match the Fortran sas_flag truncation: integration stops at q = qhigh
         # (= outermost rational q + dmlim). psiedge from dcon.in (=1.0 → no edge dW scan).
         @printf(io, "qhigh    = %.4f\n", qhigh)
-        @printf(io, "psiedge  = %.4f\n", p.psiedge)
+        @printf(io, "dW_edge_scan_start = %.4f\n", p.psiedge)
         println(io, "nn_low   = $(p.nn)")
         println(io, "nn_high  = $(p.nn)")
         println(io, "delta_mlow  = $(p.delta_mlow)")
@@ -475,7 +477,7 @@ function write_gpec_toml_file(
         println(io, "local_stability_flag = true")
         println(io, "vac_flag = true")
         @printf(io, "qlow     = %.4f\n", p.qlow)
-        @printf(io, "psiedge  = %.4f\n", p.psiedge)
+        @printf(io, "dW_edge_scan_start = %.4f\n", p.psiedge)
         println(io, "nn_low   = $(p.nn)")
         println(io, "nn_high  = $(p.nn)")
         println(io, "delta_mlow  = $(p.delta_mlow)")
@@ -1194,12 +1196,12 @@ function generate_plots(fort, julia, bench_dir, nn)
     h5_path = julia["h5_path"]
     rzphi_panels = []
     for (comp_label, f_key, h5_var) in [
-        ("ξ_R", "xi_r_fun", "perturbed_equilibrium/response/xi_R"),
-        ("ξ_Z", "xi_z_fun", "perturbed_equilibrium/response/xi_Z"),
-        ("ξ_φ", "xi_phi_fun", "perturbed_equilibrium/response/xi_phi"),
-        ("b_R", "b_r_fun", "perturbed_equilibrium/response/b_R"),
-        ("b_Z", "b_z_fun", "perturbed_equilibrium/response/b_Z"),
-        ("b_φ", "b_phi_fun", "perturbed_equilibrium/response/b_phi")]
+        ("ξ_R", "xi_r_fun", "PerturbedEquilibrium/Response/xi_R"),
+        ("ξ_Z", "xi_z_fun", "PerturbedEquilibrium/Response/xi_Z"),
+        ("ξ_φ", "xi_phi_fun", "PerturbedEquilibrium/Response/xi_phi"),
+        ("b_R", "b_r_fun", "PerturbedEquilibrium/Response/b_R"),
+        ("b_Z", "b_z_fun", "PerturbedEquilibrium/Response/b_Z"),
+        ("b_φ", "b_phi_fun", "PerturbedEquilibrium/Response/b_phi")]
         f_fun = get(fort, f_key, Matrix{ComplexF64}(undef, 0, 0))
 
         # Reconstruct theta-space from mode-space via modes_to_theta (applies ν phase + helicity)

benchmarks/benchmark_coil_ForcingTerms_against_fortran.jl

@@ -108,6 +108,7 @@ function parse_fortran_run(dir::String)::FortranRunParams
     psihigh   = _find_scalar(eq_text, "psihigh";   default=0.99)
     mtheta_eq = _find_int(eq_text,    "mtheta";    default=256)
     grid_type = _find_string(eq_text, "grid_type"; default="ldp")
+    grid_type == "ldp" && (grid_type = "rational_packed")  # map the legacy Fortran value to the Julia spelling
 
     # Toroidal mode number and m-range expansion
     nn          = _find_int(dcon_text, "nn";          default=1)

benchmarks/benchmark_delta_prime_methods.jl

@@ -23,7 +23,7 @@ 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"]["use_riccati"] = true
+    inputs["ForceFreeStates"]["integrator"] = "riccati"
     intr = FFS.ForceFreeStatesInternal(; dir_path=ex)
     ctrl = FFS.ForceFreeStatesControl(;
         (Symbol(k) => v for (k, v) in inputs["ForceFreeStates"])...)

benchmarks/benchmark_diiid_ideal_ntv_torque.jl

@@ -118,12 +118,12 @@ function build_benchmark_tomldir(eq_file::String)
 eq_filename = "$eq_name"
 eq_type = "efit"
 jac_type = "hamada"
-grid_type = "ldp"
+grid_type = "rational_packed"
 psilow = 1e-4
 psihigh = 0.993
 mpsi = 128
 mtheta = 256
-newq0 = 0
+q0_override = 0.0
 etol = 1e-7
 
 [Wall]
@@ -134,7 +134,7 @@ local_stability_flag = true
 vac_flag = true
 force_termination = false
 
-psiedge = 1.00                # No edge-scan truncation (dmlim mechanism removed in develop)
+dW_edge_scan_start = 1.00                # No edge-scan truncation (dmlim mechanism removed in develop)
 qlow = 1.02
 qhigh = 1e3
 sing_start = 0

benchmarks/benchmark_diiid_kinetic_stability.jl

@@ -8,7 +8,7 @@ Fortran GPEC's kinetic DCON reference.
 Runs `GPE.main()` with `kinetic_source="calculated"` and `kinetic_factor=1.0`
 against the EFIT g-file and `.kin` profile taken from a Fortran GPEC
 kinetic example directory, then compares the least-stable total-energy
-eigenvalue `vacuum/et[1]` against `W_t_eigenvalue[:, 0]` in the Fortran
+eigenvalue `ForceFreeStates/FreeBoundaryStability/eigenmode_energies[1]` against `W_t_eigenvalue[:, 0]` in the Fortran
 `dcon_output_n1.nc`.
 
 No inputs are duplicated into this repo — everything is read from the
@@ -129,7 +129,7 @@ psilow = 0.01
 psihigh = 0.993
 mpsi = 128
 mtheta = 256
-newq0 = 0
+q0_override = 0.0
 etol = 1e-7
 
 [Wall]
@@ -140,7 +140,7 @@ local_stability_flag = true
 vac_flag = true
 force_termination = true       # Skip PE+KF post-processing — we only need FFS eigenvalues
 
-psiedge = 1.0                 # No edge-scan truncation (dmlim mechanism removed in develop)
+dW_edge_scan_start = 1.0                 # No edge-scan truncation (dmlim mechanism removed in develop)
 qlow = 1.02
 qhigh = 1e3
 sing_start = 0
@@ -202,9 +202,9 @@ function run_benchmark(fortran_dir::String=default_fortran_dir())
     isfile(h5path) || error("Expected Julia output not found: $h5path")
 
     et = h5open(h5path, "r") do h5
-        read(h5["vacuum/et"])
+        read(h5["ForceFreeStates/FreeBoundaryStability/eigenmode_energies"])
     end
-    isempty(et) && error("vacuum/et is empty in $h5path")
+    isempty(et) && error("ForceFreeStates/FreeBoundaryStability/eigenmode_energies is empty in $h5path")
     # et is stored as a length-2*N real array (re,im interleaved) by HDF5.jl
     # when the underlying Julia array is ComplexF64. NCDatasets and HDF5 give
     # us a ComplexF64 array directly here.

benchmarks/benchmark_git_branches.jl

@@ -194,8 +194,8 @@ function run_example_benchmark(example_path, num_runs)
         end
 
         h5 = h5open(gpec_path, "r")
-        et = read(h5["FreeBoundaryStability/eigenmode_energies"])
-        nsteps = read(h5["integration/nstep"])
+        et = read(h5["ForceFreeStates/FreeBoundaryStability/eigenmode_energies"])
+        nsteps = read(h5["ForceFreeStates/Solutions/ForwardIntegration/nstep"])
         close(h5)
 
         avg_runtime = sum(runtimes) / length(runtimes)

benchmarks/benchmark_q_vs_iota_edge.jl

@@ -21,7 +21,7 @@ const EXAMPLE_DIR = joinpath(@__DIR__, "..", "examples", "DIIID-like_ideal_examp
 # Dense ldp reference equilibrium: treat its q(ψ) as ground truth
 function reference_q()
     _, eq_config, additional_input = GPE.build_inputs_from_toml(EXAMPLE_DIR)
-    eq_config.grid_type = "ldp"
+    eq_config.grid_type = "rational_packed"
     eq_config.mpsi = 1024
     equil = GPE.Equilibrium.setup_equilibrium(eq_config, additional_input)
     return equil, eq_config

benchmarks/benchmark_solovev_kinetic_stability.jl

@@ -155,7 +155,7 @@ function run_julia_reference()
     GPE.main([rundir])
     wall = time() - t0
     et = h5open(joinpath(rundir, "gpec.h5"), "r") do h5
-        read(h5["FreeBoundaryStability/eigenmode_energies"])
+        read(h5["ForceFreeStates/FreeBoundaryStability/eigenmode_energies"])
     end
     return real(et[1]), imag(et[1]), wall
 end

benchmarks/benchmark_threads.jl

@@ -1,17 +1,16 @@
 # Thread-scaling benchmark for the bidirectional parallel FM integration.
-# Runs the Solovev (N=8) and DIIID-like (N=26) examples with use_parallel=true
+# Runs the Solovev (N=8) and DIIID-like (N=26) examples with integrator="stride"
 # across 1, 2, 4, 8 threads and compares against the serial Riccati path.
 #
 # Usage (from JPEC_main root):
 #   for t in 1 2 4 8; do julia -t $t --project=. benchmarks/benchmark_threads.jl; done
 
 using GeneralizedPerturbedEquilibrium, TOML, Printf, Statistics
 
-function run_ffs(ex; use_parallel, use_riccati=false)
+function run_ffs(ex; integrator)
     inputs = TOML.parsefile(joinpath(ex, "gpec.toml"))
     inputs["ForceFreeStates"]["verbose"] = false
-    inputs["ForceFreeStates"]["use_parallel"] = use_parallel
-    inputs["ForceFreeStates"]["use_riccati"] = use_riccati
+    inputs["ForceFreeStates"]["integrator"] = integrator
     inputs["ForceFreeStates"]["write_outputs_to_HDF5"] = false
     intr = GeneralizedPerturbedEquilibrium.ForceFreeStates.ForceFreeStatesInternal(; dir_path=ex)
     ctrl = GeneralizedPerturbedEquilibrium.ForceFreeStates.ForceFreeStatesControl(;
@@ -36,17 +35,17 @@ function run_ffs(ex; use_parallel, use_riccati=false)
     return real(vac.et[1]), intr.numpert_total
 end
 
-function timed_run(ex; use_parallel, use_riccati=false, nwarm=1, nrep=2)
+function timed_run(ex; integrator, nwarm=1, nrep=2)
     # Warmup
     for _ in 1:nwarm
-        run_ffs(ex; use_parallel, use_riccati)
+        run_ffs(ex; integrator)
     end
     # Timed runs
     times = Float64[]
     local et1, N
     for _ in 1:nrep
         t0 = time()
-        et1, N = run_ffs(ex; use_parallel, use_riccati)
+        et1, N = run_ffs(ex; integrator)
         push!(times, time() - t0)
     end
     return mean(times), et1, N
@@ -60,9 +59,9 @@ 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_std, et_std, N = timed_run(ex; use_parallel=false, use_riccati=false)
-    t_ric, et_ric, _ = timed_run(ex; use_parallel=false, use_riccati=true)
-    t_par, et_par, _ = timed_run(ex; use_parallel=true, use_riccati=false)
+    t_std, et_std, N = timed_run(ex; integrator="serial")
+    t_ric, et_ric, _ = timed_run(ex; integrator="riccati")
+    t_par, et_par, _ = timed_run(ex; integrator="stride")
 
     err_ric = abs(et_ric - et_std) / abs(et_std) * 100
     err_par = abs(et_par - et_std) / abs(et_std) * 100

benchmarks/compare_gal_vs_el.jl

@@ -1,7 +1,7 @@
 # Overlay the IDEAL gal matched ξ(ψ) against the EL total-energy eigenmode ξ(ψ), for the same edge
 # eigenvector w. The ideal gal solution should reproduce the EL (DCON) ideal solution.
 #
-# EL :  ξ_EL(ψ)  = U_EL(ψ) · (U_EL_edge \ w)          (fundamental matrix, integration/xi_psi)
+# EL :  ξ_EL(ψ)  = U_EL(ψ) · (U_EL_edge \ w)          (fundamental matrix, ForceFreeStates/Solutions/ForwardIntegration/xi_psi)
 # gal:  ξ_gal(ψ) = U_gal(ψ) · w                       (identity-at-edge ⇒ coefficient is w itself)
 # w = eigenvector of the total energy operator W = W_plasma + W_vacuum (FreeBoundaryStability).
 #
@@ -19,12 +19,12 @@ 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
-    (to_c(read(f["FreeBoundaryStability/eigenmode_energies"])),
-        to_c(read(f["FreeBoundaryStability/W_freeboundary_eigenmodes"])),
-        to_c(read(f["integration/xi_psi"])), read(f["integration/psi"]),
-        to_c(read(f["galerkin/match/xi"])), read(f["galerkin/solution/psi"]),
-        Bool.(read(f["galerkin/solution/issing"])), read(f["info/mlow"]),
-        read(f["galerkin/sing_psi"]))
+    (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/issing"])), read(f["Info/mlow"]),
+        read(f["SingularSurfaces/GalerkinDeltaPrime/sing_psi"]))
 end
 
 mpert = size(u1, 1)