ALL - IMPROVEMENT - Self-describing HDF5 metadata (issue #226, PR B)

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

@@ -315,23 +315,23 @@ function load_julia_outputs(h5_path::String)
         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["psi_q"] = read(f, "Equilibrium/Profiles/psi")
         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["di"] = haskey(f, "LocalStability/D_I") ? read(f, "LocalStability/D_I") : Float64[]
+        julia["dr"] = haskey(f, "LocalStability/D_R") ? read(f, "LocalStability/D_R") : Float64[]
 
-        julia["psio"] = haskey(f, "Equilibrium/psio") ? read(f, "Equilibrium/psio") : NaN
+        julia["psio"] = haskey(f, "Equilibrium/psi_total") ? read(f, "Equilibrium/psi_total") : NaN
 
         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[]
-        julia["rational_m_res"] = haskey(f, "$sc/rational_m_res") ? read(f, "$sc/rational_m_res") : Int[]
+        julia["rational_m_res"] = haskey(f, "$sc/rational_m") ? read(f, "$sc/rational_m") : Int[]
         julia["resonant_area_weighted_field"] = haskey(f, "$sc/resonant_area_weighted_field") ? read(f, "$sc/resonant_area_weighted_field") : ComplexF64[]
         julia["resonant_current"] = haskey(f, "$sc/resonant_current") ? read(f, "$sc/resonant_current") : ComplexF64[]
         julia["island_half_width"] = haskey(f, "$sc/island_half_width") ? read(f, "$sc/island_half_width") : Float64[]
         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[]
+        julia["delta_prime"] = haskey(f, "$sc/Delta_prime") ? read(f, "$sc/Delta_prime") : ComplexF64[]
 
         pe = "PerturbedEquilibrium"
         # Fortran Phi_x/Phi_tot are the area-weighted field b̄ (tesla), matching forcing/response_b_area directly.
@@ -341,8 +341,8 @@ function load_julia_outputs(h5_path::String)
         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["clebsch_psi1"] = haskey(f, "$pe/Response/dxi_clebsch_psidpsi") ? read(f, "$pe/Response/dxi_clebsch_psidpsi") : Matrix{ComplexF64}(undef, 0, 0)
+        julia["clebsch_alpha"] = haskey(f, "$pe/Response/xi_clebsch_alpha") ? read(f, "$pe/Response/xi_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)

benchmarks/compare_gal_vs_el.jl

@@ -23,8 +23,8 @@ et, wt, u1, psiE, gxi, psiG, issing, mlow, sing_psi = h5open(h5path) do f
         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"]))
+        Bool.(read(f["ForceFreeStates/Solutions/GalerkinIntegration/Solution/is_rational"])), read(f["Info/mlow"]),
+        read(f["SingularSurfaces/GalerkinDeltaPrime/rational_psi"]))
 end
 
 mpert = size(u1, 1)

benchmarks/compare_jbgradpsi_m2.jl

@@ -19,9 +19,9 @@ 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/issing"]))
+    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/sing_psi"]), read(f["SingularSurfaces/GalerkinDeltaPrime/sing_m"]))
+        read(f["SingularSurfaces/GalerkinDeltaPrime/rational_psi"]), read(f["SingularSurfaces/GalerkinDeltaPrime/rational_m"]))
 end
 # shooting run: PE grid = ForceFreeStates/Solutions/ForwardIntegration/psi
 pa_s, psi_s = h5open(sh_h5) do f

benchmarks/plot_xi_eigenmode.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/sing_psi") ? read(f["SingularSurfaces/GalerkinDeltaPrime/sing_psi"]) : Float64[])
+        haskey(f, "SingularSurfaces/GalerkinDeltaPrime/rational_psi") ? read(f["SingularSurfaces/GalerkinDeltaPrime/rational_psi"]) : Float64[])
 end
 
 mpert, _, nstep = size(u1)

benchmarks/scan_resistivity_m2.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/issing"]))] :
+        psi = gal ? read(f["ForceFreeStates/Solutions/GalerkinIntegration/Solution/psi"])[.!Bool.(read(f["ForceFreeStates/Solutions/GalerkinIntegration/Solution/is_rational"]))] :
               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/sing_psi"]), read(f["SingularSurfaces/GalerkinDeltaPrime/sing_m"]))
+    (read(f["SingularSurfaces/GalerkinDeltaPrime/rational_psi"]), read(f["SingularSurfaces/GalerkinDeltaPrime/rational_m"]))
 end
 psi_res = mtarget in sing_m ? sing_psi[findfirst(==(mtarget), sing_m)] : NaN
 

benchmarks/scan_rotation_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/issing"]))] :
+        psi = gal ? read(f["ForceFreeStates/Solutions/GalerkinIntegration/Solution/psi"])[.!Bool.(read(f["ForceFreeStates/Solutions/GalerkinIntegration/Solution/is_rational"]))] :
               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/sing_psi"]), read(f["SingularSurfaces/GalerkinDeltaPrime/sing_m"]))
+    (read(f["SingularSurfaces/GalerkinDeltaPrime/rational_psi"]), read(f["SingularSurfaces/GalerkinDeltaPrime/rational_m"]))
 end
 psi_res = mtarget in sing_m ? sing_psi[findfirst(==(mtarget), sing_m)] : NaN
 

benchmarks/verify_gal_ideal.jl

@@ -8,10 +8,10 @@ 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
-    (to_c(read(f["ForceFreeStates/Solutions/GalerkinIntegration/Match/cout"])), to_c(read(f["ForceFreeStates/Solutions/GalerkinIntegration/Match/deltar"])),
-        to_c(read(f["ForceFreeStates/Solutions/GalerkinIntegration/Match/xi"])), to_c(read(f["ForceFreeStates/Solutions/GalerkinIntegration/Match/xi_deriv"])),
-        to_c(read(f["ForceFreeStates/Solutions/GalerkinIntegration/Solution/xi"])), to_c(read(f["ForceFreeStates/Solutions/GalerkinIntegration/Solution/xi_deriv"])),
-        read(f["SingularSurfaces/GalerkinDeltaPrime/sing_psi"]))
+    (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"]))
 end
 msing = length(sing_psi)
 mpert, ngrid, mcoil = size(mxi)

benchmarks/verify_gal_match.jl

@@ -11,10 +11,10 @@ h5path = length(ARGS) >= 1 ? ARGS[1] : "examples/DIIID-like_gal_resistive_exampl
 @info "Reading $h5path"
 
 xi, dxi, cout, cin, deltar, eig, resid, sing_psi = h5open(h5path) do f
-    (read(f["ForceFreeStates/Solutions/GalerkinIntegration/Match/xi"]), read(f["ForceFreeStates/Solutions/GalerkinIntegration/Match/xi_deriv"]),
+    (read(f["ForceFreeStates/Solutions/GalerkinIntegration/Match/xi"]), read(f["ForceFreeStates/Solutions/GalerkinIntegration/Match/dxidpsi"]),
         read(f["ForceFreeStates/Solutions/GalerkinIntegration/Match/cout"]), read(f["ForceFreeStates/Solutions/GalerkinIntegration/Match/cin"]),
-        read(f["ForceFreeStates/Solutions/GalerkinIntegration/Match/deltar"]), read(f["ForceFreeStates/Solutions/GalerkinIntegration/Match/rpec_eig"]),
-        read(f["ForceFreeStates/Solutions/GalerkinIntegration/Match/residual"]), read(f["SingularSurfaces/GalerkinDeltaPrime/sing_psi"]))
+        read(f["ForceFreeStates/Solutions/GalerkinIntegration/Match/Delta_r"]), read(f["ForceFreeStates/Solutions/GalerkinIntegration/Match/rpec_eig"]),
+        read(f["ForceFreeStates/Solutions/GalerkinIntegration/Match/residual"]), read(f["SingularSurfaces/GalerkinDeltaPrime/rational_psi"]))
 end
 # HDF5 stores ComplexF64 as a compound (re,im); convert if needed
 to_c(a) = eltype(a) <: Complex ? a : map(x -> ComplexF64(x.re, x.im), a)

benchmarks/verify_gal_solution.jl

@@ -8,10 +8,10 @@ using HDF5, Printf, Statistics
 h5path = length(ARGS) >= 1 ? ARGS[1] : "examples/DIIID-like_gal_resistive_example/gpec.h5"
 @info "Reading $h5path"
 
-psi, q, issing, xi, dxi, sing_psi = h5open(h5path) do f
-    (read(f["ForceFreeStates/Solutions/GalerkinIntegration/Solution/psi"]), read(f["ForceFreeStates/Solutions/GalerkinIntegration/Solution/q"]),
-        read(f["ForceFreeStates/Solutions/GalerkinIntegration/Solution/issing"]), read(f["ForceFreeStates/Solutions/GalerkinIntegration/Solution/xi"]),
-        read(f["ForceFreeStates/Solutions/GalerkinIntegration/Solution/xi_deriv"]), read(f["SingularSurfaces/GalerkinDeltaPrime/sing_psi"]))
+psi, issing, xi, dxi, sing_psi = h5open(h5path) do f
+    (read(f["ForceFreeStates/Solutions/GalerkinIntegration/Solution/psi"]),
+        read(f["ForceFreeStates/Solutions/GalerkinIntegration/Solution/is_rational"]), read(f["ForceFreeStates/Solutions/GalerkinIntegration/Solution/xi_psi"]),
+        read(f["ForceFreeStates/Solutions/GalerkinIntegration/Solution/dxi_psidpsi"]), read(f["SingularSurfaces/GalerkinDeltaPrime/rational_psi"]))
 end
 issing = Bool.(issing)
 mpert, ngrid, nsol = size(xi)