diff --git a/Project.toml b/Project.toml index 04edf6d0b..bb4d32cc5 100644 --- a/Project.toml +++ 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" diff --git a/benchmarks/benchmark_against_fortran_run.jl b/benchmarks/benchmark_against_fortran_run.jl index 74793ea59..6116a8058 100644 --- a/benchmarks/benchmark_against_fortran_run.jl +++ b/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) diff --git a/benchmarks/compare_gal_vs_el.jl b/benchmarks/compare_gal_vs_el.jl index cff8b26d4..16428c468 100644 --- a/benchmarks/compare_gal_vs_el.jl +++ b/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) diff --git a/benchmarks/compare_jbgradpsi_m2.jl b/benchmarks/compare_jbgradpsi_m2.jl index dd30afbd2..c31452c2d 100644 --- a/benchmarks/compare_jbgradpsi_m2.jl +++ b/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 diff --git a/benchmarks/plot_xi_eigenmode.jl b/benchmarks/plot_xi_eigenmode.jl index c4be18af7..bd3875c23 100644 --- a/benchmarks/plot_xi_eigenmode.jl +++ b/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) diff --git a/benchmarks/scan_resistivity_m2.jl b/benchmarks/scan_resistivity_m2.jl index 6a5695efd..dcb5204c9 100644 --- a/benchmarks/scan_resistivity_m2.jl +++ b/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 diff --git a/benchmarks/scan_rotation_m2.jl b/benchmarks/scan_rotation_m2.jl index bd01dfe4e..b8b539a18 100644 --- a/benchmarks/scan_rotation_m2.jl +++ b/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 diff --git a/benchmarks/verify_gal_ideal.jl b/benchmarks/verify_gal_ideal.jl index 45cf79127..9614807f0 100644 --- a/benchmarks/verify_gal_ideal.jl +++ b/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) diff --git a/benchmarks/verify_gal_match.jl b/benchmarks/verify_gal_match.jl index ee6a06290..084a6552c 100644 --- a/benchmarks/verify_gal_match.jl +++ b/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) diff --git a/benchmarks/verify_gal_solution.jl b/benchmarks/verify_gal_solution.jl index 66dac7931..d2890df39 100644 --- a/benchmarks/verify_gal_solution.jl +++ b/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) diff --git a/docs/development/hdf5-conventions.md b/docs/development/hdf5-conventions.md index 00c3fde09..a28dcdc39 100644 --- a/docs/development/hdf5-conventions.md +++ b/docs/development/hdf5-conventions.md @@ -16,12 +16,20 @@ Physics-topic groups elevated to top level (rather than nested under their produ These rules govern `gpec.h5` (and any future GPEC-produced HDF5 output); harness-internal synthetic fixtures (e.g. the `ggj/*` reference files written by `regression-harness/src/runner.jl`) are out of scope. - **Groups are CamelCase at every level** (`ForceFreeStates/`, `PerSurface/`, `GalerkinDeltaPrime/`). -- **Datasets (leaves) are snake_case** (`eigenmode_energies`, `delta_prime_matrix`). Established physics symbols keep their natural case (`E`, `F`, `Q_root_real`, `pest3_Delta`, `2piF`). +- **Datasets (leaves) are snake_case** (`eigenmode_energies`, `delta_prime_matrix`). Established physics symbols keep their natural case (`E`, `F`, `Q_root`, `pest3_Delta`, `2piF`). - **Data-driven tokens are stored verbatim**: coil-set names under `Input/RawInputs/Coils/`, KineticForces method tokens (`fgar`, …), scan indices (`Surface_`, `psi_`). +- **Word-valued names and boolean flags**: multi-word dataset names are snake_case English (`resonance_psi`, `trajectory_offsets`, `layer_widths`), never CamelCase — CamelCase is reserved for groups. A boolean flag is named for the state it asserts when true, with an `is_` prefix only where the bare word would read as a noun or collide with a data family: `is_rational` (bare `rational` would clash with the `rational_*` coordinate family) versus `enabled`, `truncated`, `no_root`, which already read as predicates. +- **Literature capitalization for physics symbols**: names match the standard literature — `D_I`, `D_R`, `Delta_prime`, `tau_R`, `tau_A` (not `di`, `dr`, `delta_prime`, `taur`); lowercase stays where the literature is lowercase (`alpha`, `q`, `beta*`, `delta_s`). +- **Scalar equilibrium parameters spell out the physics**: `R_axis`, `Z_axis`, `B_T_axis`, `a_mean`, `aspect_ratio`, `I_p`, `q_edge`, `beta_N`, `delta_upper`/`delta_lower`. The qualifier is a trailing subscript (`_axis`, `_edge`, `_wall`, `_min`, `_max`, `_upper`, `_lower`, `_extremum`), and a numbered literature definition keeps its number as the last subscript (`beta_p_1`, `l_i_2`). Fortran-era contractions (`bt0`, `amean`, `crnt`, `qa`, `betan`, `li1`) survive only as `Equilibrium.EquilibriumParameters` struct fields: `EQUIL_H5_NAMES` in `src/HDF5Schema.jl` maps each field to its dataset name, and `EQUIL_H5_SKIP` drops fields that duplicate another dataset or echo a control flag. +- **Derivatives are `ddpsi`** (`dTdpsi`, `dVdpsi`, `dqdpsi`, `dxidpsi`) — never Fortran `1` suffixes or `_deriv`. +- **"rational" over "singular"** in dataset names (`rational_psi`, `rational_q`, `rational_m`, `rational_n`, `rational_index`, `rational_count`) — kinetic/resistive runs are not singular at the rationals. Specifier order is standardized specifier-first (`rational_psi`, never `psi_rational`). +- **Vector components** follow `[d][_]_[dpsi]` — the variable always comes first and the coordinate is always the trailing subscript. A bare coordinate suffix is the **contravariant** component (`xi_psi` = ξ^ψ), `_cov_` marks the **covariant** one (`b_cov_theta` = b_θ), a leading `J` marks a **Jacobian-weighted** component (`Jxi_theta` = J·ξ^θ), and other representations sit in the same slot (`xi_clebsch_psi`, `dxi_clebsch_psidpsi`). Never drop the variable: `clebsch_psi` is wrong because it does not say *what* is being represented. There is no HDF5/netCDF standard for super- vs subscripts — flat `_` names are universal — so the typeset form always appears in the dataset's `long_name`. +- **Coordinates**: the radial abscissa is `psi` (normalized poloidal flux ψ_N) and the poloidal one is `theta` in every group; never `psi_n`, `xs`, or `ys`. +- **One name per physical quantity**: a quantity written in several groups carries the identical leaf name everywhere (`rational_psi` in `SingularSurfaces/`, `GalerkinDeltaPrime/`, and `SingularCoupling/`; `Delta_prime_matrix` in `SingularSurfaces/` and `Tearing/PerSurface/`; `dVdpsi` in `Profiles/`, `SingularSurfaces/`, and `KineticForces//`) — the group supplies the context, the leaf supplies the identity. ## Inputs live only under `Input/` -`Input/gpec_toml_raw` stores the full merged TOML, and `Input/RawInputs/` stores the raw equilibrium/forcing/coil data — together they make `gpec.h5` a self-contained rerun snapshot (`Rerun.jl` reconstructs every control struct from them; the writer/reader path pair is locked by shared `H5_*` consts in `GeneralizedPerturbedEquilibrium.jl`). **Never echo TOML flags or control-struct values into any other group** — every group outside `Input/` is derived output. (The former `kinetic/` and `slayer/settings/` echoes were removed under this rule.) +`Input/gpec_toml_raw` stores the full merged TOML, and `Input/RawInputs/` stores the raw equilibrium/forcing/coil data — together they make `gpec.h5` a self-contained rerun snapshot (`Rerun.jl` reconstructs every control struct from them; the writer and rerun reader carry cross-reference comments marking the mirrored path pair). **Never echo TOML flags or control-struct values into any other group** — every group outside `Input/` is derived output. (The former `kinetic/` and `slayer/settings/` echoes were removed under this rule.) ## Schema @@ -31,10 +39,10 @@ Top level (10 groups): |---|---| | `Info/` | Run metadata: `git_version`, mode-number ranges (`mpert`, `mlow`, …, `mn_index`), `psilim`, `qlim` | | `Input/` | Rerun snapshot: `gpec_toml_raw`, `RawInputs/{Equilibrium, ForcingTerms, Coils/}` | -| `Equilibrium/` | Scalars (β, q₀, q95, …) plus `Profiles/` (1-D: xs, 2piF, mu0p, dVdpsi, q) and `Geometry/` (2-D: rcoords, offset, nu, jac) | +| `Equilibrium/` | Scalars (`beta_N`, `q_axis`, `q_95`, `I_p`, …) plus `Profiles/` (1-D on `psi`: 2piF, mu0p, dVdpsi, q) and `Geometry/` (2-D on `psi`×`theta`: rcoords, offset, nu, jac) | | `ForceFreeStates/` | `Solutions/ForwardIntegration/` (u-solutions), `Solutions/GalerkinIntegration/` (`Solution/`, `Match/`, `msing`), `EulerLagrangeMatrices/{Ideal,Kinetic}`, `FreeBoundaryStability/`, `EdgeScan/` | -| `LocalStability/` | Mercier `di`, resistive interchange `dr`, `ballooning_Delta_prime`, ballooning α boundary | -| `SingularSurfaces/` | Per-rational-surface data: ψ, q, m/n, GGJ coefficients, `delta_prime_matrix`/`delta_prime_raw`/`delta_coil`, `GalerkinDeltaPrime/`, `Kinetic/` | +| `LocalStability/` | Mercier `D_I`, resistive interchange `D_R`, `ballooning_Delta_prime` on `psi`; the ballooning α boundary on `ballooning_psi` | +| `SingularSurfaces/` | Per-rational-surface data: `rational_psi`/`rational_q`/`rational_m`/`rational_n`, GGJ coefficients, `Delta_prime_matrix`/`Delta_prime_raw`/`Delta_coil`, `GalerkinDeltaPrime/`, `Kinetic/` | | `PerturbedEquilibrium/` | `ForcingModes/`, `Response/`, `ResponseMatrices/`, `SingularCoupling/`, `Energies/`, control-surface spectra | | `KineticForces/` | `/` (torque/energy profiles, `EnergyIntegrals/`, `KineticMatrices/`) | | `Tearing/` | `PerSurface/` (+ `DpMatrix/`), `Roots/`, `LayerWidths/`, `Diagnostics/{ValidRoots,Poles,FilteredRoots}`, `Scan/Surface_/` | @@ -42,9 +50,22 @@ Top level (10 groups): Reserved (documented, not yet written): `ForceFreeStates/Solutions/RiccatiIntegration/` — the third integrator backend slot alongside `ForwardIntegration` and `GalerkinIntegration`. +## Metadata contract (self-describing datasets) + +Every dataset outside `Input/` (raw snapshot) and `GalerkinIntegration/Match/` (debug-only) must answer "what is this, in what units, plotted against what" without opening the source — enforced by `test/runtests_h5_schema.jl`: + +- **`long_name`** — plain-text physics description. +- **`units`** — SI string (`"T"`, `"Wb/rad"`, `"A"`, `"m"`, `"J"`, `"N*m"`, `"Hz"`, `"Ohm*m"`); `"1"` for dimensionless (CF convention). Normalized quantities state the normalization in `long_name` (e.g. the power-normalized stability energies are per unit ⟨|ξ|²⟩, not joules). +- **`dims`** — required on rank ≥ 2 datasets: a greppable string like `"(psi, m)"` listing axis names in **Julia (column-major) order, axis 1 first**. Note h5py/HDFView users see file dimensions in the reversed (row-major) order. Square-matrix axes use distinct `_row`/`_col` names (`(mode_row, mode_col)`): netCDF permits repeated dimension names, but xarray mangles them on load. +- **HDF5 Dimension Scales** (the netCDF-4 coordinate mechanism): shared coordinate datasets (the `psi` grids, rational-surface `rational_psi`, geometry `psi`/`theta`) are marked with `h5ds_set_scale` and attached per-axis with `h5ds_attach_scale`/`h5ds_set_label`, so h5py `.dims`, xarray, and HDFView resolve axes natively. The H5DS C API indexes file (row-major) dimensions: Julia axis `k` of an `N`-d dataset is C index `N - k`. + +Root-level file attributes: `schema_version` (currently `"2.0"`; bump on breaking schema changes — readers dispatch on it), `Conventions = "GPEC-HDF5-2.0"`, `references`, `title` (run description), `date_created` (ISO 8601 UTC). The code version stays in `Info/git_version`. + +Mechanism: writers stay table-driven — each writer keeps a `path => (; long_name, units, dims, scale, attach)` table next to it (`scale` marks a coordinate dataset as a dimension scale; `attach` binds axes to scales — see the `Utilities.HDF5Annotations.annotate!` docstring) (`src/HDF5Schema.jl` for the main writer; alongside `write_galerkin!`, the PerturbedEquilibrium writer, `KineticForces/Output.jl`, and `Tearing/Runner/HDF5Output.jl` for the rest) and applies it post-write via `Utilities.HDF5Annotations.annotate!`. Entries for conditionally-written datasets are simply skipped when absent. When adding a dataset, add its table entry in the same commit — the schema test fails otherwise. + ## File-wide conventions -- Complex numbers are stored as the native HDF5.jl compound type (readable by h5py as a compound dtype). +- Complex quantities are stored as the native HDF5.jl compound type (readable by h5py as a compound dtype) — **never split into `*_real`/`*_imag` dataset pairs**. Sole sanctioned exception: `Input/RawInputs/ForcingTerms/amplitude_{real,imag}`, which mirrors the external forcing ingest-file format and keeps pre-existing snapshots replayable. - `NaN` is the not-computed sentinel in numeric datasets (e.g. auto-derived settings, rootless growth-rate entries). - Ragged (variable-length) data uses the flat-plus-`offsets` companion pattern (`offsets[k+1] - offsets[k]` = length of row `k`) rather than HDF5 VLEN types, e.g. `KineticForces//EnergyIntegrals/` and `Tearing/Diagnostics/*`. diff --git a/docs/src/kinetic_forces.md b/docs/src/kinetic_forces.md index 9787d9698..bfdf825ca 100644 --- a/docs/src/kinetic_forces.md +++ b/docs/src/kinetic_forces.md @@ -89,6 +89,36 @@ terms respectively, and do not modify the stored kinetic profile splines. when computing the `toroidal_rotation_factor` back-solve. Julia uses a clean reimplementation with consistent pre-scaling derivatives throughout. +## HDF5 outputs: complex torque convention and the EnergyIntegrals layout + +The method level of `KineticForces//` reports the two physical scalars a user +wants first: `total_torque` = ``T_\phi`` (N·m) and `total_energy` = ``\delta W_k`` (J). +Internally both are halves of one complex quantity ``T = T_\phi + 2in\,\delta W_k`` +(Logan 2013), and the ψ-profiles `dTdpsi` and `T` store that complex ``T`` directly — +so `imag(T)` carries the ``2n`` factor while `total_energy` has it divided out. The +per-record `EnergyIntegrals/torque` and `EnergyIntegrals/kinetic_energy` are separate +complex diagnostics of the two integrand halves at each ``(\psi, \lambda, \ell)`` +evaluation, which is why they are not packed into one number there. + +`EnergyIntegrals/` stores the variable-length integration trajectories in the +flat-plus-offsets ragged layout (chosen over HDF5 VLEN types for cross-language +support; `Tearing/Diagnostics/*` uses the same pattern). Record `k` spans +`offsets[k]+1 : offsets[k+1]` (Julia, 1-based) of each `*_all` array: + +```julia +h5open("gpec.h5", "r") do f + g = f["KineticForces/fgar/EnergyIntegrals"] + off = read(g["trajectory_offsets"]) + x_k = read(g["x_all"])[off[k]+1:off[k+1]] # record k's abscissae + I_k = read(g["integrand_all"])[off[k]+1:off[k+1]] # its complex integrand +end +``` + +```python +g = f["KineticForces/fgar/EnergyIntegrals"] # h5py, 0-based +x_k = g["x_all"][g["trajectory_offsets"][k]:g["trajectory_offsets"][k + 1]] +``` + ```@autodocs Modules = [GeneralizedPerturbedEquilibrium.KineticForces] ``` diff --git a/docs/src/perturbed_equilibrium.md b/docs/src/perturbed_equilibrium.md index dd4332f1b..e588d150f 100644 --- a/docs/src/perturbed_equilibrium.md +++ b/docs/src/perturbed_equilibrium.md @@ -16,3 +16,23 @@ GeneralizedPerturbedEquilibrium.PerturbedEquilibrium.PerturbedEquilibriumState GeneralizedPerturbedEquilibrium.PerturbedEquilibrium.compute_perturbed_equilibrium GeneralizedPerturbedEquilibrium.PerturbedEquilibrium.write_outputs_to_HDF5 ``` + +## Plotting per-surface results against ψ or q + +`SingularCoupling/` quantities are indexed by rational-surface **index**, not by q: with +multi-n runs a single q value can host several resonances, so the index is the only +unambiguous axis. Both `rational_psi` and `rational_q` are attached to that axis as HDF5 +dimension scales, so plotting against either is direct: + +```julia +h5open("gpec.h5", "r") do f + g = f["PerturbedEquilibrium/SingularCoupling"] + q = read(g["rational_q"]) + b_res = abs.(read(g["resonant_area_weighted_field"])) + scatter(q, b_res; xlabel="q", ylabel="|b^r| [T]") # or read(g["rational_psi"]) for ψ_N +end +``` + +In Python the same scales are visible through `h5py`'s dimension API +(`dset.dims[0]["psi_rational"]`, `dset.dims[0]["q_rational"]`), so xarray-style tooling can +label the axis automatically. diff --git a/docs/src/stability.md b/docs/src/stability.md index 3e0474482..ae96f9aed 100644 --- a/docs/src/stability.md +++ b/docs/src/stability.md @@ -145,12 +145,12 @@ For the derivation and implementation details behind these diagnostics, see Three diagnostics are produced and stored under the `LocalStability/` HDF5 group, each a profile in normalized poloidal flux ``\psi``: -- **Mercier criterion ``D_I``** (`LocalStability/di`) — the ideal interchange criterion. A surface +- **Mercier criterion ``D_I``** (`LocalStability/D_I`) — the ideal interchange criterion. A surface is Mercier-unstable where ``D_I > 0``. It is evaluated from the ``\det(\bar{d}_0)`` of the integrated local-mode matrix. -- **Resistive interchange ``D_R``** (`LocalStability/dr`) — the Glasser–Greene–Johnson resistive +- **Resistive interchange ``D_R``** (`LocalStability/D_R`) — the Glasser–Greene–Johnson resistive interchange criterion ``D_R = D_I + (H - 1/2)^2``. The ``D_I`` term is the same - ``\det(\bar{d}_0)`` value reported in `LocalStability/di`; ``H`` is computed from the legacy + ``\det(\bar{d}_0)`` value reported in `LocalStability/D_I`; ``H`` is computed from the legacy Mercier/GGJ flux-surface averages of the field and metric quantities. ``D_R > 0`` indicates resistive interchange instability. - **Ballooning ``\Delta'``** (`LocalStability/ballooning_Delta_prime`) — the high-``n`` ballooning @@ -161,7 +161,7 @@ in normalized poloidal flux ``\psi``: !!! note "Two different Δ' quantities" `LocalStability/ballooning_Delta_prime` is the **local high-``n`` ballooning** index and is distinct from the **resistive tearing** ``\Delta'`` described in the next section, which - is written under `SingularSurfaces/` and `PerturbedEquilibrium/SingularCoupling/delta_prime`. + is written under `SingularSurfaces/` and `PerturbedEquilibrium/SingularCoupling/Delta_prime`. They measure different instabilities; do not confuse them. ### s–α diagram @@ -231,7 +231,7 @@ where ``\Phi_R[j]`` is the forward FM product from ``\psi_{R,j-1}`` to the junct ``\Phi_L[j]`` is the backward crossing FM from ``\psi_{L,j}`` to the junction. The matrix is only populated by the parallel FM path and is written to the HDF5 output -under `SingularSurfaces/delta_prime_matrix`. +under `SingularSurfaces/Delta_prime_matrix`. ## Configuration reference diff --git a/docs/src/utilities.md b/docs/src/utilities.md index 44b80a931..f4ece086b 100644 --- a/docs/src/utilities.md +++ b/docs/src/utilities.md @@ -29,6 +29,15 @@ neoclassical models) used to set the Lundquist number in the tearing stack. Modules = [GeneralizedPerturbedEquilibrium.Utilities.NeoclassicalResistivity] ``` +## HDF5 Annotations + +Self-describing metadata for `gpec.h5` (long_name/units/dims attributes and HDF5 +Dimension Scales); see the metadata contract in `docs/development/hdf5-conventions.md`. + +```@autodocs +Modules = [GeneralizedPerturbedEquilibrium.Utilities.HDF5Annotations] +``` + ## IMAS Output ```@docs diff --git a/docs/src/workflow.md b/docs/src/workflow.md index 6ae9eef2a..228459ed3 100644 --- a/docs/src/workflow.md +++ b/docs/src/workflow.md @@ -176,7 +176,7 @@ All results are written to a single HDF5 file (default: `gpec.h5`). The top-leve |---|---| | `Info/` | Run metadata: git version, mode-number ranges, ψ limit | | `Input/` | Self-contained rerun snapshot: merged TOML blob, raw equilibrium/forcing/coil inputs | -| `Equilibrium/` | Equilibrium scalars (β, q₀, q95, …), 1-D profiles (`Profiles/`), 2-D geometry (`Geometry/`) | +| `Equilibrium/` | Equilibrium scalars (`beta_N`, `q_axis`, `q_95`, …), 1-D profiles (`Profiles/`), 2-D geometry (`Geometry/`) | | `ForceFreeStates/` | Stability solve: `Solutions/{ForwardIntegration,GalerkinIntegration}`, `EulerLagrangeMatrices/`, `FreeBoundaryStability/`, `EdgeScan/` | | `LocalStability/` | Mercier D_I, resistive interchange D_R, ballooning Δ' profiles | | `SingularSurfaces/` | Per-rational-surface data: ψ_s, q, m/n, GGJ coefficients, Δ' matrices (`GalerkinDeltaPrime/`), kinetic surfaces (`Kinetic/`) | diff --git a/examples/LAR_beta_scan/run_scan.jl b/examples/LAR_beta_scan/run_scan.jl index 95a924a62..6249ff30d 100644 --- a/examples/LAR_beta_scan/run_scan.jl +++ b/examples/LAR_beta_scan/run_scan.jl @@ -74,11 +74,11 @@ function extract_results(h5_path::String) h5open(h5_path, "r") do f fbs = "ForceFreeStates/FreeBoundaryStability" ep = read(f, "$fbs/eigenmode_plasma_energies"); ev = read(f, "$fbs/eigenmode_vacuum_energies"); et = read(f, "$fbs/eigenmode_energies") - msing = read(f, "SingularSurfaces/msing") - m_sing = read(f, "SingularSurfaces/m") - dp_mat = haskey(f, "SingularSurfaces/delta_prime_matrix") ? read(f, "SingularSurfaces/delta_prime_matrix") : nothing - qlim = haskey(f, "Info/qlim") ? read(f, "Info/qlim") : read(f, "Equilibrium/qmax") - q0 = read(f, "Equilibrium/q0"); qmax = read(f, "Equilibrium/qmax") + msing = read(f, "SingularSurfaces/rational_count") + m_sing = read(f, "SingularSurfaces/rational_m") + dp_mat = haskey(f, "SingularSurfaces/Delta_prime_matrix") ? read(f, "SingularSurfaces/Delta_prime_matrix") : nothing + qlim = haskey(f, "Info/qlim") ? read(f, "Info/qlim") : read(f, "Equilibrium/q_max") + q0 = read(f, "Equilibrium/q_axis"); qmax = read(f, "Equilibrium/q_max") dp_21 = NaN + NaN*im; dp_31 = NaN + NaN*im if dp_mat !== nothing && msing > 0 diff --git a/examples/LAR_epsilon_scan/run_scan.jl b/examples/LAR_epsilon_scan/run_scan.jl index 73f52b7c4..e3327ab95 100644 --- a/examples/LAR_epsilon_scan/run_scan.jl +++ b/examples/LAR_epsilon_scan/run_scan.jl @@ -81,11 +81,11 @@ function extract_results(h5_path::String) h5open(h5_path, "r") do f fbs = "ForceFreeStates/FreeBoundaryStability" ep = read(f, "$fbs/eigenmode_plasma_energies"); ev = read(f, "$fbs/eigenmode_vacuum_energies"); et = read(f, "$fbs/eigenmode_energies") - msing = read(f, "SingularSurfaces/msing") - m_sing = read(f, "SingularSurfaces/m") - dp_mat = haskey(f, "SingularSurfaces/delta_prime_matrix") ? read(f, "SingularSurfaces/delta_prime_matrix") : nothing - qlim = haskey(f, "Info/qlim") ? read(f, "Info/qlim") : read(f, "Equilibrium/qmax") - q0 = read(f, "Equilibrium/q0"); qmax = read(f, "Equilibrium/qmax") + msing = read(f, "SingularSurfaces/rational_count") + m_sing = read(f, "SingularSurfaces/rational_m") + dp_mat = haskey(f, "SingularSurfaces/Delta_prime_matrix") ? read(f, "SingularSurfaces/Delta_prime_matrix") : nothing + qlim = haskey(f, "Info/qlim") ? read(f, "Info/qlim") : read(f, "Equilibrium/q_max") + q0 = read(f, "Equilibrium/q_axis"); qmax = read(f, "Equilibrium/q_max") dp_21 = NaN + NaN*im; dp_31 = NaN + NaN*im if dp_mat !== nothing && msing > 0 diff --git a/regression-harness/cases/diiid_n1.toml b/regression-harness/cases/diiid_n1.toml index 6b1161ab6..a0e0a8396 100644 --- a/regression-harness/cases/diiid_n1.toml +++ b/regression-harness/cases/diiid_n1.toml @@ -93,7 +93,7 @@ order = 31 # Equilibrium [quantities.q0] -h5path = "Equilibrium/q0" +h5path = "Equilibrium/q_axis" type = "real_scalar" extract = "value" label = "q0" @@ -101,7 +101,7 @@ noise_threshold = 1e-12 order = 40 [quantities.q95] -h5path = "Equilibrium/q95" +h5path = "Equilibrium/q_95" type = "real_scalar" extract = "value" label = "q95" @@ -109,7 +109,7 @@ noise_threshold = 1e-12 order = 41 [quantities.betat] -h5path = "Equilibrium/betat" +h5path = "Equilibrium/beta_t" type = "real_scalar" extract = "value" label = "beta_t" @@ -117,7 +117,7 @@ noise_threshold = 1e-12 order = 42 [quantities.betan] -h5path = "Equilibrium/betan" +h5path = "Equilibrium/beta_N" type = "real_scalar" extract = "value" label = "beta_n" @@ -125,7 +125,7 @@ noise_threshold = 1e-12 order = 43 [quantities.li1] -h5path = "Equilibrium/li1" +h5path = "Equilibrium/l_i_1" type = "real_scalar" extract = "value" label = "internal inductance li1" @@ -133,7 +133,7 @@ noise_threshold = 1e-12 order = 44 [quantities.li2] -h5path = "Equilibrium/li2" +h5path = "Equilibrium/l_i_2" type = "real_scalar" extract = "value" label = "internal inductance li2" @@ -141,7 +141,7 @@ noise_threshold = 1e-12 order = 45 [quantities.li3] -h5path = "Equilibrium/li3" +h5path = "Equilibrium/l_i_3" type = "real_scalar" extract = "value" label = "internal inductance li3" @@ -149,7 +149,7 @@ noise_threshold = 1e-12 order = 46 [quantities.betap1] -h5path = "Equilibrium/betap1" +h5path = "Equilibrium/beta_p_1" type = "real_scalar" extract = "value" label = "poloidal beta betap1" @@ -157,7 +157,7 @@ noise_threshold = 1e-12 order = 47 [quantities.betap2] -h5path = "Equilibrium/betap2" +h5path = "Equilibrium/beta_p_2" type = "real_scalar" extract = "value" label = "poloidal beta betap2" @@ -165,7 +165,7 @@ noise_threshold = 1e-12 order = 48 [quantities.betap3] -h5path = "Equilibrium/betap3" +h5path = "Equilibrium/beta_p_3" type = "real_scalar" extract = "value" label = "poloidal beta betap3" @@ -173,7 +173,7 @@ noise_threshold = 1e-12 order = 49 [quantities.betaj] -h5path = "Equilibrium/betaj" +h5path = "Equilibrium/beta_j" type = "real_scalar" extract = "value" label = "current beta betaj" @@ -189,7 +189,7 @@ noise_threshold = 1e-12 order = 58 [quantities.crnt] -h5path = "Equilibrium/crnt" +h5path = "Equilibrium/I_p" type = "real_scalar" extract = "value" label = "plasma current" @@ -197,7 +197,7 @@ noise_threshold = 1e-12 order = 59 [quantities.bt0] -h5path = "Equilibrium/bt0" +h5path = "Equilibrium/B_T_axis" type = "real_scalar" extract = "value" label = "toroidal field bt0" @@ -205,7 +205,7 @@ noise_threshold = 1e-12 order = 62 [quantities.bwall] -h5path = "Equilibrium/bwall" +h5path = "Equilibrium/B_T_wall" type = "real_scalar" extract = "value" label = "wall field bwall" @@ -213,7 +213,7 @@ noise_threshold = 1e-12 order = 63 [quantities.aratio] -h5path = "Equilibrium/aratio" +h5path = "Equilibrium/aspect_ratio" type = "real_scalar" extract = "value" label = "aspect ratio" @@ -230,7 +230,7 @@ order = 65 # Singular surfaces [quantities.msing] -h5path = "SingularSurfaces/msing" +h5path = "SingularSurfaces/rational_count" type = "int_scalar" extract = "value" label = "# singular surfaces" @@ -238,7 +238,7 @@ noise_threshold = 0 order = 50 [quantities.sing_psi] -h5path = "SingularSurfaces/psi" +h5path = "SingularSurfaces/rational_psi" type = "real_vector" extract = "all_real" label = "singular psi locations" @@ -246,7 +246,7 @@ noise_threshold = 1e-8 order = 51 [quantities.sing_q] -h5path = "SingularSurfaces/q" +h5path = "SingularSurfaces/rational_q" type = "real_vector" extract = "all_real" label = "singular q values" @@ -271,11 +271,11 @@ noise_threshold = 0 order = 61 # Tearing stability Δ' — canonical STRIDE BVP matrix diagonal (replaces the -# previous `PerturbedEquilibrium/SingularCoupling/delta_prime` track, which +# previous `PerturbedEquilibrium/SingularCoupling/Delta_prime` track, which # was a per-surface stub computed by SingularCoupling from (rbwp1-lbwp1)/(2π·χ'). # Per-surface Δ' is de-emphasized; SingularCoupling reads this BVP matrix diagonal. [quantities.delta_prime] -h5path = "SingularSurfaces/delta_prime_matrix" +h5path = "SingularSurfaces/Delta_prime_matrix" type = "complex_matrix" extract = "diagonal_complex" label = "delta prime (BVP diagonal)" @@ -390,7 +390,7 @@ order = 71 # Local stability profiles (local_stability_flag = true). Checksums hash the full # profile so they are robust to NaN edge-surface entries and catch any change. [quantities.locstab_di] -h5path = "LocalStability/di" +h5path = "LocalStability/D_I" type = "real_vector" extract = "checksum" label = "Mercier D_I profile (checksum)" @@ -398,7 +398,7 @@ noise_threshold = 0 order = 72 [quantities.locstab_dr] -h5path = "LocalStability/dr" +h5path = "LocalStability/D_R" type = "real_vector" extract = "checksum" label = "resistive interchange D_R profile (checksum)" diff --git a/regression-harness/cases/diiid_slayer_n1.toml b/regression-harness/cases/diiid_slayer_n1.toml index 6682b97de..540ca7e50 100644 --- a/regression-harness/cases/diiid_slayer_n1.toml +++ b/regression-harness/cases/diiid_slayer_n1.toml @@ -10,7 +10,7 @@ example_dir = "examples/DIIID-like_SLAYER_example" # Per-surface SLAYER layer parameters (geometry + dimensionless) [quantities.slayer_ising] -h5path = "Tearing/PerSurface/ising" +h5path = "Tearing/PerSurface/rational_index" type = "real_vector" extract = "all_real" label = "SLAYER surface indices" @@ -74,7 +74,7 @@ noise_threshold = 1e-8 order = 17 [quantities.slayer_tauk] -h5path = "Tearing/PerSurface/tauk" +h5path = "Tearing/PerSurface/tau_k" type = "real_vector" extract = "all_real" label = "SLAYER tauk" @@ -98,24 +98,16 @@ order = 19 # surfaces only (2/1, 3/1, 4/1) via `first_3`: the Δ'/γ contour search is # numerically unreliable on the outermost surfaces (e.g. 5/1, 6/1, 7/1 near the # edge), so those are deliberately not golden-tracked. -[quantities.slayer_Q_re] -h5path = "Tearing/Roots/Q_root_real" -type = "real_vector" -extract = "first_3" -label = "SLAYER Re(Q_root) [2/1,3/1,4/1]" +[quantities.slayer_Q] +h5path = "Tearing/Roots/Q_root" +type = "complex_vector" +extract = "first_3_complex" +label = "SLAYER Q_root [2/1,3/1,4/1]" noise_threshold = 1e-4 order = 30 -[quantities.slayer_Q_im] -h5path = "Tearing/Roots/Q_root_imag" -type = "real_vector" -extract = "first_3" -label = "SLAYER Im(Q_root) [2/1,3/1,4/1]" -noise_threshold = 1e-4 -order = 31 - [quantities.slayer_omega_Hz] -h5path = "Tearing/Roots/omega_Hz" +h5path = "Tearing/Roots/omega" type = "real_vector" extract = "first_3" label = "SLAYER ω_Hz [2/1,3/1,4/1]" @@ -123,7 +115,7 @@ noise_threshold = 1.0 order = 32 [quantities.slayer_gamma_Hz] -h5path = "Tearing/Roots/gamma_Hz" +h5path = "Tearing/Roots/gamma" type = "real_vector" extract = "first_3" label = "SLAYER γ_Hz [2/1,3/1,4/1]" diff --git a/regression-harness/cases/efit_fixedbdy_separatrix.toml b/regression-harness/cases/efit_fixedbdy_separatrix.toml index 6f5d023bf..d64f09ba7 100644 --- a/regression-harness/cases/efit_fixedbdy_separatrix.toml +++ b/regression-harness/cases/efit_fixedbdy_separatrix.toml @@ -13,7 +13,7 @@ description = "Fixed-boundary EFIT (box hugs LCFS), Newton separatrix-find regre kind = "computed" [quantities.psio] -h5path = "Equilibrium/psio" +h5path = "Equilibrium/psi_total" type = "real_scalar" extract = "value" label = "psio (total poloidal flux)" @@ -21,7 +21,7 @@ noise_threshold = 1e-12 order = 10 [quantities.q0] -h5path = "Equilibrium/q0" +h5path = "Equilibrium/q_axis" type = "real_scalar" extract = "value" label = "q0" @@ -29,7 +29,7 @@ noise_threshold = 1e-12 order = 11 [quantities.q95] -h5path = "Equilibrium/q95" +h5path = "Equilibrium/q_95" type = "real_scalar" extract = "value" label = "q95" @@ -37,7 +37,7 @@ noise_threshold = 1e-12 order = 12 [quantities.betat] -h5path = "Equilibrium/betat" +h5path = "Equilibrium/beta_t" type = "real_scalar" extract = "value" label = "beta_t" @@ -45,7 +45,7 @@ noise_threshold = 1e-12 order = 13 [quantities.betan] -h5path = "Equilibrium/betan" +h5path = "Equilibrium/beta_N" type = "real_scalar" extract = "value" label = "beta_n" diff --git a/regression-harness/cases/gal_resistive_diiid.toml b/regression-harness/cases/gal_resistive_diiid.toml index e1fc42227..e9f1b23a6 100644 --- a/regression-harness/cases/gal_resistive_diiid.toml +++ b/regression-harness/cases/gal_resistive_diiid.toml @@ -11,7 +11,7 @@ example_dir = "examples/DIIID-like_gal_resistive_example" # Number of resonant surfaces in the gal domain [quantities.gal_msing] -h5path = "ForceFreeStates/Solutions/GalerkinIntegration/msing" +h5path = "ForceFreeStates/Solutions/GalerkinIntegration/rational_count" type = "int_scalar" extract = "value" label = "gal # singular surfaces" @@ -19,7 +19,7 @@ noise_threshold = 0 order = 10 [quantities.gal_sing_q] -h5path = "SingularSurfaces/GalerkinDeltaPrime/sing_q" +h5path = "SingularSurfaces/GalerkinDeltaPrime/rational_q" type = "real_vector" extract = "all_real" label = "gal singular q values" @@ -37,7 +37,7 @@ order = 20 # Full outer Δ′ matrix (nsol × 2·msing) — Frobenius norm catches any element drift [quantities.gal_delta_norm] -h5path = "SingularSurfaces/GalerkinDeltaPrime/delta" +h5path = "SingularSurfaces/GalerkinDeltaPrime/Delta_prime_raw" type = "complex_matrix" extract = "norm" label = "||gal Δ′ matrix||" @@ -46,7 +46,7 @@ order = 21 # rpec coil-response block (mpert × 2·msing) — Frobenius norm [quantities.gal_delta_coil_norm] -h5path = "SingularSurfaces/GalerkinDeltaPrime/delta_coil" +h5path = "SingularSurfaces/GalerkinDeltaPrime/Delta_coil" type = "complex_matrix" extract = "norm" label = "||gal Δ_coil block||" @@ -55,7 +55,7 @@ order = 22 # Mercier index per surface [quantities.gal_di] -h5path = "SingularSurfaces/GalerkinDeltaPrime/di" +h5path = "SingularSurfaces/GalerkinDeltaPrime/D_I" type = "real_vector" extract = "all_real" label = "gal D_I per surface" @@ -79,7 +79,7 @@ order = 31 # Inner-layer matching data Δ(Q) per surface (resist_eval geometry + GGJ inner solver). (msing × 2) [quantities.gal_match_deltar_norm] -h5path = "ForceFreeStates/Solutions/GalerkinIntegration/Match/deltar" +h5path = "ForceFreeStates/Solutions/GalerkinIntegration/Match/Delta_r" type = "complex_matrix" extract = "norm" label = "||gal inner-layer Δ||" diff --git a/regression-harness/cases/gal_resistive_pe.toml b/regression-harness/cases/gal_resistive_pe.toml index a9623f92f..9a04b30ac 100644 --- a/regression-harness/cases/gal_resistive_pe.toml +++ b/regression-harness/cases/gal_resistive_pe.toml @@ -68,7 +68,7 @@ order = 32 # Per-surface Δ' from the PE singular-coupling analysis [quantities.pe_delta_prime] -h5path = "PerturbedEquilibrium/SingularCoupling/delta_prime" +h5path = "PerturbedEquilibrium/SingularCoupling/Delta_prime" type = "complex_vector" extract = "all_complex" label = "PE Δ' per surface" diff --git a/regression-harness/cases/solovev_kinetic_calculated.toml b/regression-harness/cases/solovev_kinetic_calculated.toml index 3f10c8599..172416458 100644 --- a/regression-harness/cases/solovev_kinetic_calculated.toml +++ b/regression-harness/cases/solovev_kinetic_calculated.toml @@ -61,14 +61,14 @@ noise_threshold = 0 # Equilibrium [quantities.q0] -h5path = "Equilibrium/q0" +h5path = "Equilibrium/q_axis" type = "real_scalar" extract = "value" label = "q0" noise_threshold = 1e-12 [quantities.q95] -h5path = "Equilibrium/q95" +h5path = "Equilibrium/q_95" type = "real_scalar" extract = "value" label = "q95" @@ -76,21 +76,21 @@ noise_threshold = 1e-12 # Singular surfaces [quantities.msing] -h5path = "SingularSurfaces/msing" +h5path = "SingularSurfaces/rational_count" type = "int_scalar" extract = "value" label = "# singular surfaces" noise_threshold = 0 [quantities.sing_psi] -h5path = "SingularSurfaces/psi" +h5path = "SingularSurfaces/rational_psi" type = "real_vector" extract = "all_real" label = "singular psi locations" noise_threshold = 1e-8 [quantities.sing_q] -h5path = "SingularSurfaces/q" +h5path = "SingularSurfaces/rational_q" type = "real_vector" extract = "all_real" label = "singular q values" diff --git a/regression-harness/cases/solovev_kinetic_ntv.toml b/regression-harness/cases/solovev_kinetic_ntv.toml index cfdac6492..a9a4d5324 100644 --- a/regression-harness/cases/solovev_kinetic_ntv.toml +++ b/regression-harness/cases/solovev_kinetic_ntv.toml @@ -37,7 +37,7 @@ noise_threshold = 1e-10 order = 20 [quantities.msing] -h5path = "SingularSurfaces/msing" +h5path = "SingularSurfaces/rational_count" type = "int_scalar" extract = "value" label = "# singular surfaces" @@ -45,7 +45,7 @@ noise_threshold = 0 order = 21 [quantities.sing_psi] -h5path = "SingularSurfaces/psi" +h5path = "SingularSurfaces/rational_psi" type = "real_vector" extract = "all_real" label = "singular psi locations" @@ -53,7 +53,7 @@ noise_threshold = 1e-8 order = 22 [quantities.q0] -h5path = "Equilibrium/q0" +h5path = "Equilibrium/q_axis" type = "real_scalar" extract = "value" label = "q0" diff --git a/regression-harness/cases/solovev_kinetic_nuzero.toml b/regression-harness/cases/solovev_kinetic_nuzero.toml index d18ac2be6..ec66ab7c0 100644 --- a/regression-harness/cases/solovev_kinetic_nuzero.toml +++ b/regression-harness/cases/solovev_kinetic_nuzero.toml @@ -67,14 +67,14 @@ noise_threshold = 0 # Equilibrium [quantities.q0] -h5path = "Equilibrium/q0" +h5path = "Equilibrium/q_axis" type = "real_scalar" extract = "value" label = "q0" noise_threshold = 1e-12 [quantities.q95] -h5path = "Equilibrium/q95" +h5path = "Equilibrium/q_95" type = "real_scalar" extract = "value" label = "q95" @@ -82,21 +82,21 @@ noise_threshold = 1e-12 # Singular surfaces [quantities.msing] -h5path = "SingularSurfaces/msing" +h5path = "SingularSurfaces/rational_count" type = "int_scalar" extract = "value" label = "# singular surfaces" noise_threshold = 0 [quantities.sing_psi] -h5path = "SingularSurfaces/psi" +h5path = "SingularSurfaces/rational_psi" type = "real_vector" extract = "all_real" label = "singular psi locations" noise_threshold = 1e-8 [quantities.sing_q] -h5path = "SingularSurfaces/q" +h5path = "SingularSurfaces/rational_q" type = "real_vector" extract = "all_real" label = "singular q values" diff --git a/regression-harness/cases/solovev_multi_n.toml b/regression-harness/cases/solovev_multi_n.toml index 68aebae5c..2a6692540 100644 --- a/regression-harness/cases/solovev_multi_n.toml +++ b/regression-harness/cases/solovev_multi_n.toml @@ -76,7 +76,7 @@ order = 31 # Equilibrium [quantities.q0] -h5path = "Equilibrium/q0" +h5path = "Equilibrium/q_axis" type = "real_scalar" extract = "value" label = "q0" @@ -84,7 +84,7 @@ noise_threshold = 1e-12 order = 40 [quantities.q95] -h5path = "Equilibrium/q95" +h5path = "Equilibrium/q_95" type = "real_scalar" extract = "value" label = "q95" @@ -93,7 +93,7 @@ order = 41 # Singular surfaces [quantities.msing] -h5path = "SingularSurfaces/msing" +h5path = "SingularSurfaces/rational_count" type = "int_scalar" extract = "value" label = "# singular surfaces" @@ -101,7 +101,7 @@ noise_threshold = 0 order = 50 [quantities.sing_psi] -h5path = "SingularSurfaces/psi" +h5path = "SingularSurfaces/rational_psi" type = "real_vector" extract = "all_real" label = "singular psi locations" diff --git a/regression-harness/cases/solovev_n1.toml b/regression-harness/cases/solovev_n1.toml index d6fab13e3..1c625f22a 100644 --- a/regression-harness/cases/solovev_n1.toml +++ b/regression-harness/cases/solovev_n1.toml @@ -92,7 +92,7 @@ order = 31 # Equilibrium [quantities.q0] -h5path = "Equilibrium/q0" +h5path = "Equilibrium/q_axis" type = "real_scalar" extract = "value" label = "q0" @@ -100,7 +100,7 @@ noise_threshold = 1e-12 order = 40 [quantities.q95] -h5path = "Equilibrium/q95" +h5path = "Equilibrium/q_95" type = "real_scalar" extract = "value" label = "q95" @@ -108,7 +108,7 @@ noise_threshold = 1e-12 order = 41 [quantities.betat] -h5path = "Equilibrium/betat" +h5path = "Equilibrium/beta_t" type = "real_scalar" extract = "value" label = "beta_t" @@ -116,7 +116,7 @@ noise_threshold = 1e-12 order = 42 [quantities.betan] -h5path = "Equilibrium/betan" +h5path = "Equilibrium/beta_N" type = "real_scalar" extract = "value" label = "beta_n" @@ -125,7 +125,7 @@ order = 43 # Singular surfaces [quantities.msing] -h5path = "SingularSurfaces/msing" +h5path = "SingularSurfaces/rational_count" type = "int_scalar" extract = "value" label = "# singular surfaces" @@ -133,7 +133,7 @@ noise_threshold = 0 order = 50 [quantities.sing_psi] -h5path = "SingularSurfaces/psi" +h5path = "SingularSurfaces/rational_psi" type = "real_vector" extract = "all_real" label = "singular psi locations" @@ -141,7 +141,7 @@ noise_threshold = 1e-8 order = 51 [quantities.sing_q] -h5path = "SingularSurfaces/q" +h5path = "SingularSurfaces/rational_q" type = "real_vector" extract = "all_real" label = "singular q values" diff --git a/regression-harness/src/extractor.jl b/regression-harness/src/extractor.jl index 0c27e55ce..33a328d7c 100644 --- a/regression-harness/src/extractor.jl +++ b/regression-harness/src/extractor.jl @@ -79,10 +79,19 @@ function apply_extraction(spec::QuantitySpec, raw)::ExtractedQuantity json_str = JSON.json(arr; allownan=true) return ExtractedQuantity(name, label, nothing, nothing, json_str, "json_array", threshold) + elseif startswith(spec.extract, "first_") && endswith(spec.extract, "_complex") + # "first_N_complex": [re, im] pairs of the leading N elements of a complex + # vector — the complex-preserving sibling of "first_N". + nkeep = parse(Int, spec.extract[(length("first_")+1):(end-length("_complex"))]) + pairs = [[real(x), imag(x)] for x in raw[1:min(nkeep, length(raw))]] + json_str = JSON.json(pairs; allownan=true) + return ExtractedQuantity(name, label, nothing, nothing, json_str, "json_array", threshold) + elseif startswith(spec.extract, "first_") # "first_N": real values of the leading N vector elements only. Used to # golden-pin the inner (trustworthy) rational surfaces while ignoring # edge surfaces where the Δ'/γ contour search is numerically unreliable. + # Real data only — use "first_N_complex" for complex datasets. nkeep = parse(Int, spec.extract[(length("first_")+1):end]) arr = Float64.(real.(raw))[1:min(nkeep, length(raw))] json_str = JSON.json(arr; allownan=true) diff --git a/regression-harness/src/runner.jl b/regression-harness/src/runner.jl index 8fc4f8a68..3485d7ede 100644 --- a/regression-harness/src/runner.jl +++ b/regression-harness/src/runner.jl @@ -134,11 +134,11 @@ t_start = time() pe = Equilibrium.setup_equilibrium(cfg) elapsed = time() - t_start h5open(ARGS[1], "w") do fid - fid["Equilibrium/psio"] = pe.psio - fid["Equilibrium/q0"] = pe.params.q0 - fid["Equilibrium/q95"] = pe.params.q95 - fid["Equilibrium/betat"] = pe.params.betat - fid["Equilibrium/betan"] = pe.params.betan + fid["Equilibrium/psi_total"] = pe.psio + fid["Equilibrium/q_axis"] = pe.params.q0 + fid["Equilibrium/q_95"] = pe.params.q95 + fid["Equilibrium/beta_t"] = pe.params.betat + fid["Equilibrium/beta_N"] = pe.params.betan end %RUNINFO% """ diff --git a/src/Analysis/Equilibrium.jl b/src/Analysis/Equilibrium.jl index 9777907f7..6dbfede40 100644 --- a/src/Analysis/Equilibrium.jl +++ b/src/Analysis/Equilibrium.jl @@ -30,8 +30,8 @@ A `Plots.jl` plot object. """ function plot_qprofile(h5path; show_singular=true, save_path=nothing) xs, q, q0, q95 = h5open(h5path, "r") do fid - read(fid["Equilibrium/Profiles/xs"]), read(fid["Equilibrium/Profiles/q"]), - read(fid["Equilibrium/q0"]), read(fid["Equilibrium/q95"]) + read(fid["Equilibrium/Profiles/psi"]), read(fid["Equilibrium/Profiles/q"]), + read(fid["Equilibrium/q_axis"]), read(fid["Equilibrium/q_95"]) end p = plot( @@ -51,7 +51,7 @@ function plot_qprofile(h5path; show_singular=true, save_path=nothing) if show_singular msing, psi_sing, q_sing = h5open(h5path, "r") do fid - read(fid["SingularSurfaces/msing"]), read(fid["SingularSurfaces/psi"]), read(fid["SingularSurfaces/q"]) + read(fid["SingularSurfaces/rational_count"]), read(fid["SingularSurfaces/rational_psi"]), read(fid["SingularSurfaces/rational_q"]) end for s in 1:msing vline!(p, [psi_sing[s]]; linestyle=:dash, color=:red, label=nothing) @@ -83,8 +83,8 @@ A `Plots.jl` plot object. """ function plot_pressure_profile(h5path; save_path=nothing) xs, mu0p, msing, psi_sing = h5open(h5path, "r") do fid - read(fid["Equilibrium/Profiles/xs"]), read(fid["Equilibrium/Profiles/mu0p"]), - read(fid["SingularSurfaces/msing"]), read(fid["SingularSurfaces/psi"]) + read(fid["Equilibrium/Profiles/psi"]), read(fid["Equilibrium/Profiles/mu0p"]), + read(fid["SingularSurfaces/rational_count"]), read(fid["SingularSurfaces/rational_psi"]) end p = plot( @@ -124,8 +124,8 @@ A `Plots.jl` plot object. """ function plot_f_profile(h5path; save_path=nothing) xs, twopif, msing, psi_sing = h5open(h5path, "r") do fid - read(fid["Equilibrium/Profiles/xs"]), read(fid["Equilibrium/Profiles/2piF"]), - read(fid["SingularSurfaces/msing"]), read(fid["SingularSurfaces/psi"]) + read(fid["Equilibrium/Profiles/psi"]), read(fid["Equilibrium/Profiles/2piF"]), + read(fid["SingularSurfaces/rational_count"]), read(fid["SingularSurfaces/rational_psi"]) end p = plot( @@ -172,9 +172,9 @@ A `Plots.jl` plot object. function plot_flux_surfaces(h5path; n_psi=11, n_theta=18, save_path=nothing) rcoords, offset_data, xs_rz, ys_rz, ro, zo, msing, psi_sing, q_sing = h5open(h5path, "r") do fid read(fid["Equilibrium/Geometry/rcoords"]), read(fid["Equilibrium/Geometry/offset"]), - read(fid["Equilibrium/Geometry/xs"]), read(fid["Equilibrium/Geometry/ys"]), - read(fid["Equilibrium/ro"]), read(fid["Equilibrium/zo"]), - read(fid["SingularSurfaces/msing"]), read(fid["SingularSurfaces/psi"]), read(fid["SingularSurfaces/q"]) + read(fid["Equilibrium/Geometry/psi"]), read(fid["Equilibrium/Geometry/theta"]), + read(fid["Equilibrium/R_axis"]), read(fid["Equilibrium/Z_axis"]), + read(fid["SingularSurfaces/rational_count"]), read(fid["SingularSurfaces/rational_psi"]), read(fid["SingularSurfaces/rational_q"]) end n_psi_grid = length(xs_rz) @@ -361,9 +361,9 @@ A `Plots.jl` plot object. """ function plot_equilibrium_summary(h5path; save_path=nothing) q0, q95, betat, betan, kappa, li1 = h5open(h5path, "r") do fid - read(fid["Equilibrium/q0"]), read(fid["Equilibrium/q95"]), - read(fid["Equilibrium/betat"]), read(fid["Equilibrium/betan"]), - read(fid["Equilibrium/kappa"]), read(fid["Equilibrium/li1"]) + read(fid["Equilibrium/q_axis"]), read(fid["Equilibrium/q_95"]), + read(fid["Equilibrium/beta_t"]), read(fid["Equilibrium/beta_N"]), + read(fid["Equilibrium/kappa"]), read(fid["Equilibrium/l_i_1"]) end title_str = "q0=$(round(q0,digits=2)) q95=$(round(q95,digits=2)) βₜ=$(round(betat,digits=3)) βₙ=$(round(betan,digits=3)) κ=$(round(kappa,digits=2)) li1=$(round(li1,digits=3))" diff --git a/src/Analysis/ForceFreeStates.jl b/src/Analysis/ForceFreeStates.jl index c92d6e698..f276c058a 100644 --- a/src/Analysis/ForceFreeStates.jl +++ b/src/Analysis/ForceFreeStates.jl @@ -129,7 +129,7 @@ function plot_energy_eigenvectors(h5path; matrix_type=:total, save_path=nothing) error("matrix_type=$matrix_type not supported; only :total has eigenvector matrix stored in HDF5 (ep/ev are eigenvalue vectors, not matrices)") wt, psio, mlow = h5open(h5path, "r") do fid - read(fid["ForceFreeStates/FreeBoundaryStability/W_freeboundary_eigenmodes"]), read(fid["Equilibrium/psio"]), read(fid["Info/mlow"]) + read(fid["ForceFreeStates/FreeBoundaryStability/W_freeboundary_eigenmodes"]), read(fid["Equilibrium/psi_total"]), read(fid["Info/mlow"]) end isempty(wt) && error("No vacuum data in $h5path; rerun with vac_flag = true") @@ -321,9 +321,9 @@ A `Plots.jl` plot object. """ function plot_delta_prime(h5path; save_path=nothing) msing, psi_sing, q_sing, ca_l, ca_r, psio, mn_index = h5open(h5path, "r") do fid - read(fid["SingularSurfaces/msing"]), read(fid["SingularSurfaces/psi"]), read(fid["SingularSurfaces/q"]), + read(fid["SingularSurfaces/rational_count"]), read(fid["SingularSurfaces/rational_psi"]), read(fid["SingularSurfaces/rational_q"]), read(fid["SingularSurfaces/ca_left"]), read(fid["SingularSurfaces/ca_right"]), - read(fid["Equilibrium/psio"]), read(fid["Info/mn_index"]) + read(fid["Equilibrium/psi_total"]), read(fid["Info/mn_index"]) end msing == 0 && return plot(; title="No singular surfaces found", legend=false) @@ -374,7 +374,7 @@ end Plot the BALOO-style infinite-n ballooning stability diagram: the experimental pressure gradient α (solid) and the first stability boundary α_crit (dashed) versus normalized poloidal flux ψ_N. Surfaces where the experimental α lies above the boundary -are ballooning-unstable. Reads `LocalStability/psi`, `LocalStability/alpha`, and +are ballooning-unstable. Reads `LocalStability/ballooning_psi`, `LocalStability/alpha`, and `LocalStability/alpha_critical` (populated when ForceFreeStates runs with `local_stability_flag = true`). @@ -395,7 +395,7 @@ A `Plots.jl` plot object. function plot_ballooning_alpha_boundary(h5path; save_path=nothing, psi_min=0.0) psi, alpha, alpha_crit = h5open(h5path, "r") do fid haskey(fid, "LocalStability/alpha") || return (Float64[], Float64[], Float64[]) - read(fid["LocalStability/psi"]), read(fid["LocalStability/alpha"]), read(fid["LocalStability/alpha_critical"]) + read(fid["LocalStability/ballooning_psi"]), read(fid["LocalStability/alpha"]), read(fid["LocalStability/alpha_critical"]) end isempty(alpha) && return plot(; title="No local stability data (set local_stability_flag)", legend=false) @@ -465,9 +465,9 @@ function plot_cond_fbar(h5path; save_path=nothing, zoom=false) read(kg["scan_cond"]), read(kg["scan_threshold"]), read(kg["psi"]), - read(fid["SingularSurfaces/psi"]), - read(fid["SingularSurfaces/q"]), - read(kg["kmsing"])) + read(fid["SingularSurfaces/rational_psi"]), + read(fid["SingularSurfaces/rational_q"]), + read(kg["rational_count"])) end if isempty(scan_psi) diff --git a/src/Analysis/PerturbedEquilibrium.jl b/src/Analysis/PerturbedEquilibrium.jl index 0e8675055..33f8bdc0f 100644 --- a/src/Analysis/PerturbedEquilibrium.jl +++ b/src/Analysis/PerturbedEquilibrium.jl @@ -184,15 +184,15 @@ end plot_driven_delta_prime(h5path; save_path=nothing) Scatter plot of `Re(Δ')` per rational surface vs ψ_N, computed by the perturbed -equilibrium module (from `PerturbedEquilibrium/SingularCoupling/delta_prime`). +equilibrium module (from `PerturbedEquilibrium/SingularCoupling/Delta_prime`). One marker series per toroidal mode n. Integer-valued q rational surfaces are annotated. This is the forcing-driven Δ' (response to the applied perturbation amplitudes in `intr.forcing_modes`); for the equilibrium-intrinsic Δ' from the STRIDE BVP, -read `SingularSurfaces/delta_prime_matrix` from the HDF5 directly. +read `SingularSurfaces/Delta_prime_matrix` from the HDF5 directly. -Requires `PerturbedEquilibrium/SingularCoupling/delta_prime` in the HDF5 file. +Requires `PerturbedEquilibrium/SingularCoupling/Delta_prime` in the HDF5 file. ### Arguments @@ -208,11 +208,11 @@ A `Plots.jl` plot object. """ function plot_driven_delta_prime(h5path; save_path=nothing) base = "PerturbedEquilibrium/SingularCoupling/" - _has_pe_data(h5path, base * "delta_prime") || + _has_pe_data(h5path, base * "Delta_prime") || return plot(; title="No PE Δ' data — run with perturbed equilibrium enabled", legend=false) delta_prime, rational_psi, rational_q, rational_n = h5open(h5path, "r") do fid - read(fid[base * "delta_prime"]), + read(fid[base * "Delta_prime"]), read(fid[base * "rational_psi"]), read(fid[base * "rational_q"]), read(fid[base * "rational_n"]) @@ -370,7 +370,7 @@ function plot_mode_spectrogram(h5path; component=:xi_psi, save_path=nothing) read(fid[dataset_path]), read(fid["ForceFreeStates/Solutions/ForwardIntegration/psi"]), read(fid["Info/mlow"]), read(fid["Info/mhigh"]), read(fid["Info/nhigh"]), - read(fid["Equilibrium/q95"]), + read(fid["Equilibrium/q_95"]), read(fid["PerturbedEquilibrium/SingularCoupling/rational_psi"]) end diff --git a/src/Analysis/PerturbedEquilibriumModes.jl b/src/Analysis/PerturbedEquilibriumModes.jl index b5e5a8e45..58aacdaf9 100644 --- a/src/Analysis/PerturbedEquilibriumModes.jl +++ b/src/Analysis/PerturbedEquilibriumModes.jl @@ -66,8 +66,8 @@ function modes_to_theta(h5_file::String, variable::String; if !keep_sfl_phi # Reconstruct ν spline from stored grid + nodal values (FastInterpolations v0.4 API) - rzphi_xs = read(f, "Equilibrium/Geometry/xs") - rzphi_ys = read(f, "Equilibrium/Geometry/ys") + rzphi_xs = read(f, "Equilibrium/Geometry/psi") + rzphi_ys = read(f, "Equilibrium/Geometry/theta") nu_vals = read(f, "Equilibrium/Geometry/nu") nu_spline = cubic_interp( (rzphi_xs, rzphi_ys), nu_vals; @@ -77,8 +77,8 @@ function modes_to_theta(h5_file::String, variable::String; psi_grid = read(f, "ForceFreeStates/Solutions/ForwardIntegration/psi") - bt_sign = haskey(f, "Equilibrium/bt_sign") ? read(f, "Equilibrium/bt_sign") : 1 - crnt = haskey(f, "Equilibrium/crnt") ? read(f, "Equilibrium/crnt") : 1.0 + bt_sign = haskey(f, "Equilibrium/B_T_sign") ? read(f, "Equilibrium/B_T_sign") : 1 + crnt = haskey(f, "Equilibrium/I_p") ? read(f, "Equilibrium/I_p") : 1.0 helicity = bt_sign * Int(sign(crnt)) hint = (Ref(1), Ref(1)) diff --git a/src/ForceFreeStates/ForceFreeStatesStructs.jl b/src/ForceFreeStates/ForceFreeStatesStructs.jl index 00eaf0ba7..5d7cbb920 100644 --- a/src/ForceFreeStates/ForceFreeStatesStructs.jl +++ b/src/ForceFreeStates/ForceFreeStatesStructs.jl @@ -252,7 +252,7 @@ gpec.toml. - `use_riccati::Bool` - Use the dual Riccati reformulation S = U₁·U₂⁻¹ instead of the standard U₁/U₂ ODE. Reduces stiffness for faster integration. See Glasser (2018) Phys. Plasmas 25, 032507. - `use_parallel::Bool` - Parallel fundamental matrix (propagator) integration using `Threads.@threads`. Each chunk is integrated independently from identity IC and assembled serially. Requires `singfac_min != 0`. Uses the same chunk bounds as the standard path but sub-divides chunks for load balancing. Crossings use the Riccati-style algorithm (no Gaussian reduction). - `parallel_threads::Int` - Cap on the number of threads the parallel BVP uses. **Default `2`** parallelises the FM chunks across two threads (the BVP has ~10 chunks; 2 threads is enough to amortize them — speedup saturates here, raising to 4 adds scheduling overhead). Set `parallel_threads = 1` to run the FM chunks SERIALLY (no `Threads.@threads`), which is bit-deterministic and immune to the thread-schedule sensitivity that can cause intermittent BVP divergence on numerically delicate equilibria. The parallel path produces bit-identical Δ′ across thread counts; `parallel_threads = 2` is about 20% faster than serial and saturates the speedup. If a parallel run diverges, drop to `parallel_threads = 1` rather than switching `use_parallel = false` — the latter is silently wrong. Capped at `Threads.nthreads()`. - - `populate_dense_xi::Bool` - When `use_parallel = true`, append a serial Euler-Lagrange pass after the propagator BVP so the returned `odet` carries dense axis-basis `u_store`/`du_store`/`xi_s_store` — the only convention PerturbedEquilibrium / FieldReconstruction consume correctly, and what fills HDF5 `ForceFreeStates/Solutions/ForwardIntegration/xi_*`. Δ' (`SingularSurfaces/delta_prime_matrix`) comes from the BVP and is identical either way; free-boundary energies (`ForceFreeStates/FreeBoundaryStability/eigenmode_*`) match a pure serial run when enabled (~0.12 % Riccati-vs-axis gap otherwise). **Default `false`** (skips the extra serial pass, ~1× BVP wall-clock); **PerturbedEquilibrium-using configs must set it `true`** when `use_parallel = true`, else PE reads Riccati-basis garbage. Auto-disabled when `force_termination = true`. + - `populate_dense_xi::Bool` - When `use_parallel = true`, append a serial Euler-Lagrange pass after the propagator BVP so the returned `odet` carries dense axis-basis `u_store`/`du_store`/`xi_s_store` — the only convention PerturbedEquilibrium / FieldReconstruction consume correctly, and what fills HDF5 `ForceFreeStates/Solutions/ForwardIntegration/xi_*`. Δ' (`SingularSurfaces/Delta_prime_matrix`) comes from the BVP and is identical either way; free-boundary energies (`ForceFreeStates/FreeBoundaryStability/eigenmode_*`) match a pure serial run when enabled (~0.12 % Riccati-vs-axis gap otherwise). **Default `false`** (skips the extra serial pass, ~1× BVP wall-clock); **PerturbedEquilibrium-using configs must set it `true`** when `use_parallel = true`, else PE reads Riccati-basis garbage. Auto-disabled when `force_termination = true`. - `extended_precision_bvp::Bool` - When `true` (default), promote the Δ' BVP linear system to `Complex{Double64}` (~31 digits) for the LU solve and PEST3 combination. Guards against catastrophic cancellation in the PEST3 four-term combination (dp_raw entries can be 10⁴–10⁵× larger than the result; the imaginary part of off-diagonal Δ' is particularly sensitive). Disabling (`false`) saves ~1.5–2× the BVP solve time but on DIIID-class equilibria the imaginary Δ' components can drift by factors of 2–5×; only disable for performance experiments on cases where Float64 has been validated against Double64. """ @kwdef struct ForceFreeStatesControl @@ -291,7 +291,7 @@ gpec.toml. save_interval::Int = 3 force_termination::Bool = false use_riccati::Bool = false - use_parallel::Bool = true # Default on: unlocks SingularSurfaces/delta_prime_matrix (STRIDE BVP Δ′ matrix) used by SLAYER/GGJ downstream. + use_parallel::Bool = true # Default on: unlocks SingularSurfaces/Delta_prime_matrix (STRIDE BVP Δ′ matrix) used by SLAYER/GGJ downstream. populate_dense_xi::Bool = false # When use_parallel=true, set to true ONLY if a PerturbedEquilibrium pipeline will consume dense ξ. Default false avoids the ~1× parallel-BVP serial-EL re-run for non-PE runs (Δ'/vacuum/ideal-stability only). See ForceFreeStatesControl docstring for the full trade-off (et[1] convention differs by ~0.12% on DIIID between populate=true vs false). extended_precision_bvp::Bool = true # Promote Δ' BVP to Complex{Double64}; default on (Float64 drifts the imaginary Δ' by 2–5× on DIIID-class cases). diff --git a/src/ForceFreeStates/Galerkin/GalerkinSolve.jl b/src/ForceFreeStates/Galerkin/GalerkinSolve.jl index 3802bb857..a67f1c797 100644 --- a/src/ForceFreeStates/Galerkin/GalerkinSolve.jl +++ b/src/ForceFreeStates/Galerkin/GalerkinSolve.jl @@ -238,30 +238,32 @@ ASCII/binary outputs. function write_galerkin!(out_h5, result::GalerkinResult) gal = "ForceFreeStates/Solutions/GalerkinIntegration" gdp = "SingularSurfaces/GalerkinDeltaPrime" - out_h5["$gal/msing"] = result.msing - result.msing == 0 && return nothing - out_h5["$gdp/delta"] = result.delta + out_h5["$gal/rational_count"] = result.msing + if result.msing == 0 + annotate_galerkin!(out_h5) + return nothing + end + out_h5["$gdp/Delta_prime_raw"] = result.delta out_h5["$gdp/pest3_A"] = result.Ap out_h5["$gdp/pest3_B"] = result.Bp out_h5["$gdp/pest3_Gamma"] = result.Gammap out_h5["$gdp/pest3_Delta"] = result.Deltap - out_h5["$gdp/sing_psi"] = result.sing_psi - out_h5["$gdp/sing_q"] = result.sing_q - out_h5["$gdp/sing_m"] = result.sing_m - out_h5["$gdp/sing_n"] = result.sing_n - out_h5["$gdp/di"] = result.di + out_h5["$gdp/rational_psi"] = result.sing_psi + out_h5["$gdp/rational_q"] = result.sing_q + out_h5["$gdp/rational_m"] = result.sing_m + out_h5["$gdp/rational_n"] = result.sing_n + out_h5["$gdp/D_I"] = result.di out_h5["$gdp/alpha"] = result.alpha if !isempty(result.delta_coil) - out_h5["$gdp/delta_coil"] = result.delta_coil + out_h5["$gdp/Delta_coil"] = result.delta_coil end if result.solution !== nothing sol = result.solution out_h5["$gal/Solution/psi"] = sol.psi - out_h5["$gal/Solution/q"] = sol.q - out_h5["$gal/Solution/issing"] = collect(sol.issing) - out_h5["$gal/Solution/xi"] = sol.xi - out_h5["$gal/Solution/xi_deriv"] = sol.xi_deriv - isempty(sol.xi_cut) || (out_h5["$gal/Solution/xi_cut"] = sol.xi_cut) + out_h5["$gal/Solution/is_rational"] = collect(sol.issing) + out_h5["$gal/Solution/xi_psi"] = sol.xi + out_h5["$gal/Solution/dxi_psidpsi"] = sol.xi_deriv + isempty(sol.xi_cut) || (out_h5["$gal/Solution/xi_psi_cut"] = sol.xi_cut) isempty(sol.cut_range) || (out_h5["$gal/Solution/cut_range"] = sol.cut_range) end if result.match !== nothing @@ -269,8 +271,8 @@ function write_galerkin!(out_h5, result::GalerkinResult) out_h5["$gal/Match/cout"] = m.cout out_h5["$gal/Match/cin"] = m.cin out_h5["$gal/Match/xi"] = m.xi - out_h5["$gal/Match/xi_deriv"] = m.xi_deriv - out_h5["$gal/Match/deltar"] = m.deltar + out_h5["$gal/Match/dxidpsi"] = m.xi_deriv + out_h5["$gal/Match/Delta_r"] = m.deltar out_h5["$gal/Match/bpen"] = m.bpen out_h5["$gal/Match/rpec_eig"] = m.rpec_eig # Per-surface inner-layer ξ_ψ(ψ) (match.f intotsol); ragged grids → one dataset pair per surface. @@ -281,10 +283,57 @@ function write_galerkin!(out_h5, result::GalerkinResult) end out_h5["$gal/Match/residual"] = m.residual if !isempty(m.inner_params) - for f in (:E, :F, :G, :H, :K, :M, :taua, :taur, :v1) + for f in (:E, :F, :G, :H, :K, :M) out_h5["$gal/Match/InnerParams/$(f)"] = [getfield(pp, f) for pp in m.inner_params] end + # Literature names, matching the Tearing PerSurface mapping for the same fields. + out_h5["$gal/Match/InnerParams/tau_A"] = [pp.taua for pp in m.inner_params] + out_h5["$gal/Match/InnerParams/tau_R"] = [pp.taur for pp in m.inner_params] + out_h5["$gal/Match/InnerParams/dVdpsi"] = [pp.v1 for pp in m.inner_params] end end + annotate_galerkin!(out_h5) + return nothing +end + +# Metadata tables for the Galerkin outputs (Match/** is debug-only and exempt from the +# metadata contract; see docs/development/hdf5-conventions.md). +const GALERKIN_H5_ANNOTATIONS = [ + "ForceFreeStates/Solutions/GalerkinIntegration/rational_count" => (; long_name="number of rational (singular) surfaces in the Galerkin solve"), + "ForceFreeStates/Solutions/GalerkinIntegration/Solution/psi" => (; long_name="normalized poloidal flux ψ_N grid of the Galerkin solution", scale="psi"), + "ForceFreeStates/Solutions/GalerkinIntegration/Solution/is_rational" => + (; long_name="flag: grid node lies on a rational surface", dims=("psi",), attach=(1 => "ForceFreeStates/Solutions/GalerkinIntegration/Solution/psi",)), + "ForceFreeStates/Solutions/GalerkinIntegration/Solution/xi_psi" => + (; long_name="Galerkin solution functions ξ^ψ (arbitrary amplitude; note the psi/solution axis order differs from ForwardIntegration)", dims=("mode", "psi", "solution")), + "ForceFreeStates/Solutions/GalerkinIntegration/Solution/dxi_psidpsi" => + (; long_name="ψ_N derivative of the Galerkin solution functions ξ^ψ (arbitrary amplitude)", dims=("mode", "psi", "solution")), + "ForceFreeStates/Solutions/GalerkinIntegration/Solution/xi_psi_cut" => + (; long_name="Galerkin solution functions ξ^ψ with the leading-order resonant response excised", dims=("mode", "psi", "solution")), + "ForceFreeStates/Solutions/GalerkinIntegration/Solution/cut_range" => + (; long_name="ψ_N bounds of the excised resonant + extension cells per surface", dims=("surface", "bound")), + "SingularSurfaces/GalerkinDeltaPrime/Delta_prime_raw" => + (; long_name="outer-region Δ' matrix (2msing×2msing, side-major [L_s1, R_s1, ...])", dims=("surface_side_row", "surface_side_col")), + "SingularSurfaces/GalerkinDeltaPrime/pest3_A" => (; long_name="PEST-3 matching block A' (Galerkin outer region)", dims=("surface_row", "surface_col")), + "SingularSurfaces/GalerkinDeltaPrime/pest3_B" => (; long_name="PEST-3 matching block B' (Galerkin outer region)", dims=("surface_row", "surface_col")), + "SingularSurfaces/GalerkinDeltaPrime/pest3_Gamma" => (; long_name="PEST-3 matching block Γ' (Galerkin outer region)", dims=("surface_row", "surface_col")), + "SingularSurfaces/GalerkinDeltaPrime/pest3_Delta" => (; long_name="PEST-3 matching block Δ' (Galerkin outer region)", dims=("surface_row", "surface_col")), + "SingularSurfaces/GalerkinDeltaPrime/rational_psi" => (; long_name="normalized poloidal flux ψ_N of each rational surface", scale="psi_rational"), + "SingularSurfaces/GalerkinDeltaPrime/rational_q" => + (; long_name="safety factor q = m/n at each rational surface", dims=("surface",), attach=(1 => "SingularSurfaces/GalerkinDeltaPrime/rational_psi",)), + "SingularSurfaces/GalerkinDeltaPrime/rational_m" => + (; long_name="resonant poloidal mode number m at each rational surface", dims=("surface",), attach=(1 => "SingularSurfaces/GalerkinDeltaPrime/rational_psi",)), + "SingularSurfaces/GalerkinDeltaPrime/rational_n" => + (; long_name="resonant toroidal mode number n at each rational surface", dims=("surface",), attach=(1 => "SingularSurfaces/GalerkinDeltaPrime/rational_psi",)), + "SingularSurfaces/GalerkinDeltaPrime/D_I" => + (; long_name="Mercier D_I at each rational surface", dims=("surface",), attach=(1 => "SingularSurfaces/GalerkinDeltaPrime/rational_psi",)), + "SingularSurfaces/GalerkinDeltaPrime/alpha" => + (; long_name="Frobenius small-solution exponent α at each rational surface", dims=("surface",), attach=(1 => "SingularSurfaces/GalerkinDeltaPrime/rational_psi",)), + "SingularSurfaces/GalerkinDeltaPrime/Delta_coil" => (; long_name="edge coil-response matrix (edge mode × surface-side; RPEC columns)", dims=("mode", "surface_side")) +] + +# Attach long_name/units/dims + dimension scales (declared in-table) to everything +# write_galerkin! wrote. +function annotate_galerkin!(out_h5) + Utilities.HDF5Annotations.annotate!(out_h5, GALERKIN_H5_ANNOTATIONS) return nothing end diff --git a/src/GeneralizedPerturbedEquilibrium.jl b/src/GeneralizedPerturbedEquilibrium.jl index 519dd952c..95f64277d 100755 --- a/src/GeneralizedPerturbedEquilibrium.jl +++ b/src/GeneralizedPerturbedEquilibrium.jl @@ -59,6 +59,7 @@ include("Analysis/Analysis.jl") import .Analysis as Analysis export Analysis +include("HDF5Schema.jl") include("Rerun.jl") # Import ForceFreeStates types and functions needed for main @@ -651,7 +652,7 @@ function main_from_inputs( if kf_ctrl.write_outputs_to_HDF5 h5open(joinpath(intr.dir_path, kf_ctrl.HDF5_filename), "cw") do h5file - KineticForces.write_to_hdf5!(h5file, kf_state) + KineticForces.write_to_hdf5!(h5file, kf_state; dVdpsi_spline=equil.profiles.dVdpsi_spline) end end end @@ -718,6 +719,9 @@ function write_outputs_to_HDF5( h5open(joinpath(intr.dir_path, ctrl.HDF5_filename), "w") do out_h5 + # File-level metadata contract (schema_version, Conventions, title, date). + Utilities.HDF5Annotations.write_root_attrs!(out_h5; title="GPEC output: $(basename(abspath(intr.dir_path)))") + # Store git version for reproducibility out_h5["Info/git_version"] = git_version @@ -756,27 +760,30 @@ function write_outputs_to_HDF5( out_h5["Info/mn_index"] = hcat(m, n) # (N, 2) matrix out_h5["Info/psilim"] = intr.psilim out_h5["Info/qlim"] = intr.qlim - out_h5["Info/q1lim"] = intr.q1lim - - # Write derived equilibrium parameters - for (key, val) in zip(fieldnames(Equilibrium.EquilibriumParameters), getfield.(Ref(equil.params), fieldnames(Equilibrium.EquilibriumParameters))) - if val !== nothing # TODO: looks like ro, zo, psio, and b_norm are not set, so skipping those for now but should fix eventually - out_h5["Equilibrium/$key"] = val - end + out_h5["Info/dqdpsi_lim"] = intr.q1lim + + # Write derived equilibrium parameters. The struct keeps its legacy field spellings; + # EQUIL_H5_NAMES maps them to literature dataset names and EQUIL_H5_SKIP drops + # duplicates and control-flag echoes. Fields left `nothing` are not written. + for f in fieldnames(Equilibrium.EquilibriumParameters) + f in EQUIL_H5_SKIP && continue + val = getfield(equil.params, f) + val === nothing && continue + out_h5["Equilibrium/$(get(EQUIL_H5_NAMES, f, String(f)))"] = val end - out_h5["Equilibrium/psio"] = equil.psio - out_h5["Equilibrium/ro"] = equil.ro - out_h5["Equilibrium/zo"] = equil.zo + out_h5["Equilibrium/psi_total"] = equil.psio + out_h5["Equilibrium/R_axis"] = equil.ro + out_h5["Equilibrium/Z_axis"] = equil.zo # Write equilibrium profile and geometry arrays (from the named splines) profiles = equil.profiles - out_h5["Equilibrium/Profiles/xs"] = profiles.xs + out_h5["Equilibrium/Profiles/psi"] = profiles.xs out_h5["Equilibrium/Profiles/2piF"] = profiles.F_spline.y out_h5["Equilibrium/Profiles/mu0p"] = profiles.P_spline.y out_h5["Equilibrium/Profiles/dVdpsi"] = profiles.dVdpsi_spline.y out_h5["Equilibrium/Profiles/q"] = profiles.q_spline.y - out_h5["Equilibrium/Geometry/xs"] = equil.rzphi_xs - out_h5["Equilibrium/Geometry/ys"] = equil.rzphi_ys + out_h5["Equilibrium/Geometry/psi"] = equil.rzphi_xs + out_h5["Equilibrium/Geometry/theta"] = equil.rzphi_ys # Extract grid point values from interpolants for HDF5 output out_h5["Equilibrium/Geometry/rcoords"] = equil.rzphi_rsquared.nodal_derivs.partials[1, :, :] out_h5["Equilibrium/Geometry/offset"] = equil.rzphi_offset.nodal_derivs.partials[1, :, :] @@ -784,23 +791,26 @@ function write_outputs_to_HDF5( out_h5["Equilibrium/Geometry/jac"] = equil.rzphi_jac.nodal_derivs.partials[1, :, :] # Write local stability data; always write all entries, using empty arrays when not computed. - # LocalStability/di = Mercier D_I (det(d0bar)); LocalStability/dr = resistive interchange D_R; + # LocalStability/D_I = Mercier D_I (det(d0bar)); LocalStability/D_R = resistive interchange D_R; # LocalStability/ballooning_Delta_prime = high-n ballooning Δ' (distinct from the Riccati - # tearing Δ' under PerturbedEquilibrium/SingularCoupling/delta_prime). + # tearing Δ' under PerturbedEquilibrium/SingularCoupling/Delta_prime). if locstab !== nothing locstab_xs = locstab.cache.x - out_h5["LocalStability/di"] = locstab.y[:, 1] ./ locstab_xs - out_h5["LocalStability/dr"] = locstab.y[:, 2] ./ locstab_xs + out_h5["LocalStability/psi"] = collect(locstab_xs) # cached vector → dense for HDF5 + out_h5["LocalStability/D_I"] = locstab.y[:, 1] ./ locstab_xs + out_h5["LocalStability/D_R"] = locstab.y[:, 2] ./ locstab_xs else - out_h5["LocalStability/di"] = Float64[] - out_h5["LocalStability/dr"] = Float64[] + out_h5["LocalStability/psi"] = Float64[] + out_h5["LocalStability/D_I"] = Float64[] + out_h5["LocalStability/D_R"] = Float64[] end - out_h5["SingularSurfaces/di0"] = (locstab !== nothing && !isempty(intr.sing)) ? + out_h5["SingularSurfaces/D_I"] = (locstab !== nothing && !isempty(intr.sing)) ? [locstab(sing.psifac)[1] / sing.psifac for sing in intr.sing] : Float64[] out_h5["LocalStability/ballooning_Delta_prime"] = locstab !== nothing ? locstab.y[:, 4] : Float64[] # First ballooning stability boundary: experimental α vs critical α (BALOO-style). - out_h5["LocalStability/psi"] = ballooning_boundary.psi + # Its own scan grid, distinct from the LocalStability/psi profile grid above. + out_h5["LocalStability/ballooning_psi"] = ballooning_boundary.psi out_h5["LocalStability/alpha"] = ballooning_boundary.alpha out_h5["LocalStability/alpha_critical"] = ballooning_boundary.alpha_critical @@ -812,7 +822,7 @@ function write_outputs_to_HDF5( out_h5["$fwd/q"] = odet.q_store out_h5["$fwd/xi_psi"] = odet.u_store[:, :, 1, :] out_h5["$fwd/u2"] = odet.u_store[:, :, 2, :] # TODO: what to name this? These are the "conjugate momenta" of u1 - out_h5["$fwd/dxi_psi"] = odet.du_store + out_h5["$fwd/dxi_psidpsi"] = odet.du_store out_h5["$fwd/xi_s"] = odet.xi_s_store out_h5["$fwd/crit"] = odet.crit_store @@ -830,10 +840,10 @@ function write_outputs_to_HDF5( end # Write singular surface data - out_h5["SingularSurfaces/msing"] = intr.msing - out_h5["SingularSurfaces/psi"] = [sing.psifac for sing in intr.sing] - out_h5["SingularSurfaces/q"] = [sing.q for sing in intr.sing] - out_h5["SingularSurfaces/q1"] = [sing.q1 for sing in intr.sing] + out_h5["SingularSurfaces/rational_count"] = intr.msing + out_h5["SingularSurfaces/rational_psi"] = [sing.psifac for sing in intr.sing] + out_h5["SingularSurfaces/rational_q"] = [sing.q for sing in intr.sing] + out_h5["SingularSurfaces/dqdpsi"] = [sing.q1 for sing in intr.sing] out_h5["SingularSurfaces/ca_left"] = odet.ca_l out_h5["SingularSurfaces/ca_right"] = odet.ca_r @@ -848,8 +858,8 @@ function write_outputs_to_HDF5( n_matrix[s, i] = sing.n[i] end end - out_h5["SingularSurfaces/m"] = m_matrix - out_h5["SingularSurfaces/n"] = n_matrix + out_h5["SingularSurfaces/rational_m"] = m_matrix + out_h5["SingularSurfaces/rational_n"] = n_matrix # Glasser-Greene-Johnson geometric coefficients + surface averages # (populated by ForceFreeStates.resist_eval_all! after sing_find!). @@ -866,9 +876,9 @@ function write_outputs_to_HDF5( out_h5["SingularSurfaces/M"] = [s.restype.M for s in intr.sing] out_h5["SingularSurfaces/avg_bsq_over_dpsisq"] = [s.restype.avg_bsq_over_dpsisq for s in intr.sing] out_h5["SingularSurfaces/avg_bsq"] = [s.restype.avg_bsq for s in intr.sing] - out_h5["SingularSurfaces/p_local"] = [s.restype.p_local for s in intr.sing] - out_h5["SingularSurfaces/p1_local"] = [s.restype.p1_local for s in intr.sing] - out_h5["SingularSurfaces/v1_local"] = [s.restype.v1_local for s in intr.sing] + out_h5["SingularSurfaces/mu0p"] = [s.restype.p_local for s in intr.sing] + out_h5["SingularSurfaces/dmu0pdpsi"] = [s.restype.p1_local for s in intr.sing] + out_h5["SingularSurfaces/dVdpsi"] = [s.restype.v1_local for s in intr.sing] end end @@ -877,16 +887,16 @@ function write_outputs_to_HDF5( # Write inter-surface Δ' matrix if computed (parallel FM path only). # Shape: [msing × msing] — PEST3-convention deltap (STRIDE BVP with vacuum coupling). if intr.msing > 0 && !isempty(intr.delta_prime_matrix) - out_h5["SingularSurfaces/delta_prime_matrix"] = intr.delta_prime_matrix + out_h5["SingularSurfaces/Delta_prime_matrix"] = intr.delta_prime_matrix end # Edge coil-response matrix, stored (numpert_total × 2msing) = (edge mode, surface-side) to match - # the SingularSurfaces/GalerkinDeltaPrime/delta_coil layout so H5Web heatmaps share axes + # the SingularSurfaces/GalerkinDeltaPrime/Delta_coil layout so H5Web heatmaps share axes # (x = edge mode, y = surface-side). # Internal intr.delta_coil stays (2msing × numpert_total); transpose only at write. if intr.msing > 0 && !isempty(intr.delta_coil) dc = permutedims(intr.delta_coil) - out_h5["SingularSurfaces/delta_coil"] = dc + out_h5["SingularSurfaces/Delta_coil"] = dc end # Write raw 2msing×2msing outer-region D' matrix in side-major ordering @@ -895,15 +905,15 @@ function write_outputs_to_HDF5( # Needed for the full det(D' − D(γ)) = 0 eigenvalue problem via # pest3_decompose to recover (A', B', Γ', Δ'). if intr.msing > 0 && !isempty(intr.delta_prime_raw) - out_h5["SingularSurfaces/delta_prime_raw"] = intr.delta_prime_raw + out_h5["SingularSurfaces/Delta_prime_raw"] = intr.delta_prime_raw end # Write kinetic singular surface data (det(F̄) near-zeros) and the cond(F̄) scan # used to find them. Populated only when kinetic crossings were searched for. - out_h5["SingularSurfaces/Kinetic/kmsing"] = intr.kmsing - out_h5["SingularSurfaces/Kinetic/psi"] = [s.psifac for s in intr.kinsing] - out_h5["SingularSurfaces/Kinetic/q"] = [s.q for s in intr.kinsing] - out_h5["SingularSurfaces/Kinetic/q1"] = [s.q1 for s in intr.kinsing] + out_h5["SingularSurfaces/Kinetic/rational_count"] = intr.kmsing + out_h5["SingularSurfaces/Kinetic/rational_psi"] = [s.psifac for s in intr.kinsing] + out_h5["SingularSurfaces/Kinetic/rational_q"] = [s.q for s in intr.kinsing] + out_h5["SingularSurfaces/Kinetic/dqdpsi"] = [s.q1 for s in intr.kinsing] out_h5["SingularSurfaces/Kinetic/scan_psi"] = intr.kinsing_scan_psi out_h5["SingularSurfaces/Kinetic/scan_cond"] = intr.kinsing_scan_cond out_h5["SingularSurfaces/Kinetic/scan_threshold"] = intr.kinsing_scan_threshold @@ -984,6 +994,9 @@ function write_outputs_to_HDF5( out_h5["$elm/Kinetic/G"] = _eval_mat_spline(ffit.gaats) end end + + # Self-describing metadata pass (long_name/units/dims + dimension scales). + apply_main_h5_metadata!(out_h5, intr) end end diff --git a/src/HDF5Schema.jl b/src/HDF5Schema.jl new file mode 100644 index 000000000..a8beacffa --- /dev/null +++ b/src/HDF5Schema.jl @@ -0,0 +1,286 @@ +# Metadata tables for the datasets written by write_outputs_to_HDF5 (the main gpec.h5 +# writer), applied post-write by Utilities.HDF5Annotations.annotate!. Sub-writers +# (Galerkin, PerturbedEquilibrium, KineticForces, Tearing) keep their tables next to +# their own writers. Paths absent from a given run are skipped automatically. +# +# Conventions (docs/development/hdf5-conventions.md): units are SI strings, "1" for +# dimensionless; ψ always means the normalized poloidal flux ψ_N ∈ [0, 1]; stability +# energies are power-normalized (per unit surface-averaged |ξ|², not joules); `dims` +# lists axis names in Julia (column-major) order, axis 1 first. + +# Equilibrium scalars are written by iterating EquilibriumParameters' fields. The struct keeps +# its legacy field spellings; these tables give the HDF5 datasets their literature names and +# drop fields that duplicate another dataset or echo a control flag (inputs live only under Input/). +const EQUIL_H5_NAMES = Dict( + :ro => "R_axis", + :zo => "Z_axis", + :b0 => "B_axis", + :bt0 => "B_T_axis", + :bwall => "B_T_wall", + :rmean => "R_mean", + :amean => "a_mean", + :aratio => "aspect_ratio", + :delta1 => "delta_upper", + :delta2 => "delta_lower", + :crnt => "I_p", + :rsep => "R_midplane", + :rext => "R_extremum", + :zext => "Z_extremum", + :q0 => "q_axis", + :qmin => "q_min", + :qmax => "q_max", + :qa => "q_edge", + :q95 => "q_95", + :qextrema_psi => "q_extrema_psi", + :qextrema_q => "q_extrema_q", + :mextrema => "q_extrema_count", + :betat => "beta_t", + :betan => "beta_N", + :betap1 => "beta_p_1", + :betap2 => "beta_p_2", + :betap3 => "beta_p_3", + :betaj => "beta_j", + :li1 => "l_i_1", + :li2 => "l_i_2", + :li3 => "l_i_3", + :bt_sign => "B_T_sign", + :psio => "psi_total" +) +const EQUIL_H5_SKIP = Set([:psi0, :psi_axis, :psi_axis_norm, :zsep, :verbose, :diagnose_src, :diagnose_maxima]) + +const MAIN_H5_ANNOTATIONS = [ + # --- Info/ --- + "Info/git_version" => (; long_name="GPEC git version that produced this file"), + "Info/mpert" => (; long_name="number of poloidal harmonics per toroidal mode"), + "Info/mlow" => (; long_name="lowest poloidal mode number m"), + "Info/mhigh" => (; long_name="highest poloidal mode number m"), + "Info/npert" => (; long_name="number of toroidal mode numbers"), + "Info/nlow" => (; long_name="lowest toroidal mode number n"), + "Info/nhigh" => (; long_name="highest toroidal mode number n"), + "Info/mn_index" => (; long_name="(m, n) mode numbers for each perturbation index", dims=("mode", "m_or_n")), + "Info/psilim" => (; long_name="normalized poloidal flux at the integration boundary"), + "Info/qlim" => (; long_name="safety factor q at the integration boundary"), + "Info/dqdpsi_lim" => (; long_name="dq/dψ_N at the integration boundary"), + # --- Equilibrium/ scalars (written per-field when set; superset listed) --- + "Equilibrium/R_axis" => (; long_name="R-coordinate of the magnetic axis", units="m"), + "Equilibrium/Z_axis" => (; long_name="Z-coordinate of the magnetic axis", units="m"), + "Equilibrium/psi_total" => (; long_name="total poloidal flux difference |ψ_axis − ψ_boundary|", units="Wb/rad"), + "Equilibrium/R_midplane" => (; long_name="R of the boundary at the inboard and outboard midplane crossings", units="m"), + "Equilibrium/R_extremum" => (; long_name="R of the boundary at its upper and lower Z extrema", units="m"), + "Equilibrium/Z_extremum" => (; long_name="Z of the boundary at its upper and lower Z extrema", units="m"), + "Equilibrium/B_axis" => (; long_name="total magnetic field strength at the axis", units="T"), + "Equilibrium/q_axis" => (; long_name="safety factor at the magnetic axis"), + "Equilibrium/q_min" => (; long_name="minimum safety factor in the plasma"), + "Equilibrium/q_max" => (; long_name="maximum safety factor in the plasma"), + "Equilibrium/q_edge" => (; long_name="safety factor at the plasma edge"), + "Equilibrium/q_95" => (; long_name="safety factor at the 95% flux surface"), + "Equilibrium/q_extrema_psi" => (; long_name="normalized poloidal flux at q-profile extrema"), + "Equilibrium/q_extrema_q" => (; long_name="safety factor at q-profile extrema"), + "Equilibrium/q_extrema_count" => (; long_name="number of extrema in the q-profile"), + "Equilibrium/R_mean" => (; long_name="mean major radius of the plasma", units="m"), + "Equilibrium/a_mean" => (; long_name="mean minor radius of the plasma", units="m"), + "Equilibrium/aspect_ratio" => (; long_name="aspect ratio R0/a"), + "Equilibrium/kappa" => (; long_name="plasma elongation"), + "Equilibrium/delta_upper" => (; long_name="upper triangularity"), + "Equilibrium/delta_lower" => (; long_name="lower triangularity"), + "Equilibrium/B_T_axis" => (; long_name="toroidal field at the axis", units="T"), + "Equilibrium/I_p" => (; long_name="plasma current", units="A"), + "Equilibrium/B_T_wall" => (; long_name="toroidal field at the wall", units="T"), + "Equilibrium/beta_t" => (; long_name="toroidal beta"), + "Equilibrium/beta_N" => (; long_name="normalized beta β_N"), + "Equilibrium/beta_p_1" => (; long_name="poloidal beta (definition 1)"), + "Equilibrium/beta_p_2" => (; long_name="poloidal beta (definition 2)"), + "Equilibrium/beta_p_3" => (; long_name="poloidal beta (definition 3)"), + "Equilibrium/beta_j" => (; long_name="current-weighted beta"), + "Equilibrium/l_i_1" => (; long_name="internal inductance (definition 1)"), + "Equilibrium/l_i_2" => (; long_name="internal inductance (definition 2)"), + "Equilibrium/l_i_3" => (; long_name="internal inductance (definition 3)"), + "Equilibrium/volume" => (; long_name="plasma volume", units="m^3"), + "Equilibrium/B_T_sign" => (; long_name="sign of the toroidal field"), + "Equilibrium/psi_norm" => (; long_name="normalized poloidal flux at the magnetic axis (0 by definition of ψ_N)"), + "Equilibrium/psi_boundary" => (; long_name="poloidal flux at the plasma boundary in the internal normalized convention (1 by construction, not a Wb/rad datum)"), + "Equilibrium/psi_boundary_norm" => (; long_name="normalized poloidal flux at the plasma boundary (1 by definition of ψ_N)"), + "Equilibrium/psi_axis_offset" => (; long_name="offset applied to the axis poloidal flux", units="Wb/rad"), + "Equilibrium/psi_boundary_offset" => (; long_name="offset applied to the boundary poloidal flux", units="Wb/rad"), + "Equilibrium/psi_axis_sign" => (; long_name="sign of the axis poloidal flux"), + "Equilibrium/psi_boundary_sign" => (; long_name="sign of the boundary poloidal flux"), + "Equilibrium/psi_boundary_zero" => (; long_name="flag: boundary poloidal flux is zero"), + # --- Equilibrium/Profiles/ (1-D profiles on the ψ_N grid) --- + "Equilibrium/Profiles/psi" => (; long_name="normalized poloidal flux ψ_N profile grid", scale="psi"), + "Equilibrium/Profiles/2piF" => (; long_name="2π F with F = R B_φ the toroidal field function", units="T*m", dims=("psi",), attach=(1 => "Equilibrium/Profiles/psi",)), + "Equilibrium/Profiles/mu0p" => (; long_name="μ0 × plasma pressure", units="T^2", dims=("psi",), attach=(1 => "Equilibrium/Profiles/psi",)), + "Equilibrium/Profiles/dVdpsi" => (; long_name="flux-surface volume derivative dV/dψ_N", units="m^3", dims=("psi",), attach=(1 => "Equilibrium/Profiles/psi",)), + "Equilibrium/Profiles/q" => (; long_name="safety factor profile", dims=("psi",), attach=(1 => "Equilibrium/Profiles/psi",)), + # --- Equilibrium/Geometry/ (2-D flux-coordinate maps on (ψ_N, θ/2π)) --- + "Equilibrium/Geometry/psi" => (; long_name="normalized poloidal flux ψ_N geometry grid", scale="psi"), + "Equilibrium/Geometry/theta" => (; long_name="normalized poloidal angle θ/2π geometry grid", scale="theta"), + "Equilibrium/Geometry/rcoords" => + (; long_name="squared minor-radius coordinate r² of the working coordinate map", units="m^2", dims=("psi", "theta"), + attach=(1 => "Equilibrium/Geometry/psi", 2 => "Equilibrium/Geometry/theta")), + "Equilibrium/Geometry/offset" => + (; long_name="poloidal-angle offset of the working coordinate map (fraction of 2π)", dims=("psi", "theta"), + attach=(1 => "Equilibrium/Geometry/psi", 2 => "Equilibrium/Geometry/theta")), + "Equilibrium/Geometry/nu" => + (; long_name="toroidal-angle offset ν = φ − 2πζ of the working coordinate map", units="rad", dims=("psi", "theta"), + attach=(1 => "Equilibrium/Geometry/psi", 2 => "Equilibrium/Geometry/theta")), + "Equilibrium/Geometry/jac" => + (; long_name="Jacobian of the (ψ_N, θ, ζ) working coordinates", units="m^3", dims=("psi", "theta"), + attach=(1 => "Equilibrium/Geometry/psi", 2 => "Equilibrium/Geometry/theta")), + # --- LocalStability/ --- + "LocalStability/psi" => (; long_name="normalized poloidal flux ψ_N of the local-stability profiles", scale="psi"), + "LocalStability/D_I" => (; long_name="Mercier ideal interchange criterion D_I", dims=("psi",), attach=(1 => "LocalStability/psi",)), + "LocalStability/D_R" => + (; long_name="Glasser-Greene-Johnson resistive interchange criterion D_R", dims=("psi",), attach=(1 => "LocalStability/psi",)), + "LocalStability/ballooning_Delta_prime" => + (; long_name="high-n ballooning Δ' (distinct from the tearing Δ')", dims=("psi",), attach=(1 => "LocalStability/psi",)), + "LocalStability/ballooning_psi" => + (; long_name="normalized poloidal flux ψ_N of the ballooning α boundary scan", scale="psi_ballooning"), + "LocalStability/alpha" => + (; long_name="experimental normalized pressure gradient α", dims=("psi_ballooning",), attach=(1 => "LocalStability/ballooning_psi",)), + "LocalStability/alpha_critical" => + (; long_name="critical normalized pressure gradient α for first ballooning stability", dims=("psi_ballooning",), + attach=(1 => "LocalStability/ballooning_psi",)), + # --- ForceFreeStates/Solutions/ForwardIntegration/ --- + "ForceFreeStates/Solutions/ForwardIntegration/nstep" => (; long_name="number of saved solution snapshots"), + "ForceFreeStates/Solutions/ForwardIntegration/nstep_total" => (; long_name="total ODE solver steps taken"), + "ForceFreeStates/Solutions/ForwardIntegration/psi" => (; long_name="normalized poloidal flux ψ_N at saved solution snapshots", scale="psi"), + "ForceFreeStates/Solutions/ForwardIntegration/q" => + (; long_name="safety factor at saved solution snapshots", dims=("psi",), attach=(1 => "ForceFreeStates/Solutions/ForwardIntegration/psi",)), + "ForceFreeStates/Solutions/ForwardIntegration/xi_psi" => + (; long_name="fundamental-matrix solutions ξ^ψ (arbitrary amplitude)", dims=("mode", "solution", "psi"), attach=(3 => "ForceFreeStates/Solutions/ForwardIntegration/psi",)), + "ForceFreeStates/Solutions/ForwardIntegration/u2" => + (; long_name="conjugate momenta of the fundamental-matrix solutions (arbitrary amplitude)", dims=("mode", "solution", "psi")), + "ForceFreeStates/Solutions/ForwardIntegration/dxi_psidpsi" => + (; long_name="ψ_N derivative of the fundamental-matrix solutions ξ^ψ (arbitrary amplitude)", dims=("mode", "solution", "psi")), + "ForceFreeStates/Solutions/ForwardIntegration/xi_s" => + (; long_name="Clebsch surface-displacement solutions Ξ_s (arbitrary amplitude)", dims=("mode", "solution", "psi"), + attach=(3 => "ForceFreeStates/Solutions/ForwardIntegration/psi",)), + "ForceFreeStates/Solutions/ForwardIntegration/crit" => + (; long_name="DCON zero-crossing criterion at saved snapshots", dims=("psi",), attach=(1 => "ForceFreeStates/Solutions/ForwardIntegration/psi",)), + # --- ForceFreeStates/EdgeScan/ (power-normalized (W, N) pencil energies) --- + "ForceFreeStates/EdgeScan/psi" => (; long_name="normalized poloidal flux ψ_N of the edge truncation scan", scale="psi"), + "ForceFreeStates/EdgeScan/q" => (; long_name="safety factor at scan points", dims=("psi",), attach=(1 => "ForceFreeStates/EdgeScan/psi",)), + "ForceFreeStates/EdgeScan/total_energy" => + (; long_name="power-normalized total energy of the least-stable free-boundary mode (per unit ⟨|ξ|²⟩)", dims=("psi",), attach=(1 => "ForceFreeStates/EdgeScan/psi",)), + "ForceFreeStates/EdgeScan/plasma_energy" => + (; long_name="power-normalized plasma energy of the least-stable mode (per unit ⟨|ξ|²⟩)", dims=("psi",), attach=(1 => "ForceFreeStates/EdgeScan/psi",)), + "ForceFreeStates/EdgeScan/vacuum_energy" => + (; long_name="power-normalized vacuum energy of the least-stable mode (per unit ⟨|ξ|²⟩)", dims=("psi",), attach=(1 => "ForceFreeStates/EdgeScan/psi",)), + "ForceFreeStates/EdgeScan/vacuum_eigenvalue" => + (; long_name="least vacuum eigenvalue of the (W, N) pencil at scan points", dims=("psi",), attach=(1 => "ForceFreeStates/EdgeScan/psi",)), + # --- SingularSurfaces/ --- + "SingularSurfaces/rational_count" => (; long_name="number of rational (singular) surfaces in the domain"), + "SingularSurfaces/rational_psi" => (; long_name="normalized poloidal flux ψ_N of each rational surface", scale="psi_rational"), + "SingularSurfaces/rational_q" => + (; long_name="safety factor q = m/n at each rational surface", dims=("surface",), scale="q_rational", attach=(1 => "SingularSurfaces/rational_psi",)), + "SingularSurfaces/dqdpsi" => + (; long_name="dq/dψ_N at each rational surface", dims=("surface",), attach=(1 => "SingularSurfaces/rational_psi", 1 => "SingularSurfaces/rational_q")), + "SingularSurfaces/rational_m" => (; long_name="resonant poloidal mode numbers per surface (0-padded)", dims=("surface", "mode")), + "SingularSurfaces/rational_n" => (; long_name="resonant toroidal mode numbers per surface (0-padded)", dims=("surface", "mode")), + "SingularSurfaces/D_I" => + (; long_name="Mercier D_I evaluated at each rational surface", dims=("surface",), attach=(1 => "SingularSurfaces/rational_psi", 1 => "SingularSurfaces/rational_q")), + "SingularSurfaces/ca_left" => + (; long_name="asymptotic large/small-solution coefficient matrices just left of each surface", dims=("mode", "solution", "large_small", "surface")), + "SingularSurfaces/ca_right" => + (; long_name="asymptotic large/small-solution coefficient matrices just right of each surface", dims=("mode", "solution", "large_small", "surface")), + "SingularSurfaces/E" => + (; long_name="Glasser-Greene-Johnson coefficient E per surface", dims=("surface",), attach=(1 => "SingularSurfaces/rational_psi", 1 => "SingularSurfaces/rational_q")), + "SingularSurfaces/F" => + (; long_name="Glasser-Greene-Johnson coefficient F per surface", dims=("surface",), attach=(1 => "SingularSurfaces/rational_psi", 1 => "SingularSurfaces/rational_q")), + "SingularSurfaces/G" => + (; long_name="Glasser-Greene-Johnson coefficient G per surface", dims=("surface",), attach=(1 => "SingularSurfaces/rational_psi", 1 => "SingularSurfaces/rational_q")), + "SingularSurfaces/H" => + (; long_name="Glasser-Greene-Johnson coefficient H per surface", dims=("surface",), attach=(1 => "SingularSurfaces/rational_psi", 1 => "SingularSurfaces/rational_q")), + "SingularSurfaces/K" => + (; long_name="Glasser-Greene-Johnson coefficient K per surface", dims=("surface",), attach=(1 => "SingularSurfaces/rational_psi", 1 => "SingularSurfaces/rational_q")), + "SingularSurfaces/M" => + (; long_name="Glasser-Greene-Johnson coefficient M per surface", dims=("surface",), attach=(1 => "SingularSurfaces/rational_psi", 1 => "SingularSurfaces/rational_q")), + "SingularSurfaces/avg_bsq_over_dpsisq" => + (; long_name="flux-surface average ⟨B²/|∇ψ_N|²⟩ per surface", units="T^2*m^2", dims=("surface",), + attach=(1 => "SingularSurfaces/rational_psi", 1 => "SingularSurfaces/rational_q")), + "SingularSurfaces/avg_bsq" => + (; long_name="flux-surface average ⟨B²⟩ per surface", units="T^2", dims=("surface",), attach=(1 => "SingularSurfaces/rational_psi", 1 => "SingularSurfaces/rational_q")), + "SingularSurfaces/mu0p" => + (; long_name="μ0 × local pressure at each surface", units="T^2", dims=("surface",), attach=(1 => "SingularSurfaces/rational_psi", 1 => "SingularSurfaces/rational_q")), + "SingularSurfaces/dmu0pdpsi" => + (; long_name="μ0 × dp/dψ_N at each surface", units="T^2", dims=("surface",), attach=(1 => "SingularSurfaces/rational_psi", 1 => "SingularSurfaces/rational_q")), + "SingularSurfaces/dVdpsi" => + (; long_name="dV/dψ_N at each surface", units="m^3", dims=("surface",), attach=(1 => "SingularSurfaces/rational_psi", 1 => "SingularSurfaces/rational_q")), + "SingularSurfaces/Delta_prime_matrix" => (; long_name="inter-surface Δ' matrix (PEST3 convention, STRIDE BVP with vacuum coupling)", dims=("surface_row", "surface_col")), + "SingularSurfaces/Delta_prime_raw" => + (; long_name="raw 2msing×2msing outer-region D' matrix, side-major ordering [L_s1, R_s1, ...]", dims=("surface_side_row", "surface_side_col")), + "SingularSurfaces/Delta_coil" => (; long_name="edge coil-response matrix (edge mode × surface-side)", dims=("mode", "surface_side")), + # --- SingularSurfaces/Kinetic/ --- + "SingularSurfaces/Kinetic/rational_count" => (; long_name="number of kinetic singular surfaces (det(F̄) near-zeros)"), + "SingularSurfaces/Kinetic/rational_psi" => (; long_name="normalized poloidal flux ψ_N of kinetic singular surfaces", scale="psi_kinetic_rational"), + "SingularSurfaces/Kinetic/rational_q" => + (; long_name="safety factor at kinetic singular surfaces", dims=("surface",), attach=(1 => "SingularSurfaces/Kinetic/rational_psi",)), + "SingularSurfaces/Kinetic/dqdpsi" => + (; long_name="dq/dψ_N at kinetic singular surfaces", dims=("surface",), attach=(1 => "SingularSurfaces/Kinetic/rational_psi",)), + "SingularSurfaces/Kinetic/scan_psi" => (; long_name="ψ_N grid of the cond(F̄) scan", scale="psi_scan"), + "SingularSurfaces/Kinetic/scan_cond" => + (; long_name="condition number of F̄ along the scan", dims=("psi_scan",), attach=(1 => "SingularSurfaces/Kinetic/scan_psi",)), + "SingularSurfaces/Kinetic/scan_threshold" => (; long_name="cond(F̄) threshold used to flag kinetic singular surfaces"), + # --- ForceFreeStates/FreeBoundaryStability/ (power-normalized (W, N) pencil) --- + "ForceFreeStates/FreeBoundaryStability/W_freeboundary" => (; long_name="power-normalized free-boundary energy matrix W (per unit ⟨|ξ|²⟩)", dims=("mode_row", "mode_col")), + "ForceFreeStates/FreeBoundaryStability/W_plasma" => (; long_name="power-normalized plasma energy matrix (per unit ⟨|ξ|²⟩)", dims=("mode_row", "mode_col")), + "ForceFreeStates/FreeBoundaryStability/W_vacuum" => (; long_name="power-normalized vacuum energy matrix (per unit ⟨|ξ|²⟩)", dims=("mode_row", "mode_col")), + "ForceFreeStates/FreeBoundaryStability/W_freeboundary_eigenmodes" => + (; long_name="generalized eigenvectors of the (W, N) pencil, columns sorted most-unstable first, unit power norm", dims=("mode", "eigenmode")), + "ForceFreeStates/FreeBoundaryStability/eigenmode_energies" => + (; long_name="generalized eigenvalues of the (W, N) pencil: total energy per unit ⟨|ξ|²⟩, coordinate-invariant", dims=("eigenmode",)), + "ForceFreeStates/FreeBoundaryStability/eigenmode_plasma_energies" => (; long_name="plasma contribution to the power-normalized eigenmode energies", dims=("eigenmode",)), + "ForceFreeStates/FreeBoundaryStability/eigenmode_vacuum_energies" => (; long_name="vacuum contribution to the power-normalized eigenmode energies", dims=("eigenmode",)), + "ForceFreeStates/FreeBoundaryStability/vacuum_eigenvalue" => (; long_name="least eigenvalue of the vacuum energy matrix"), + # --- SurfaceGeometries/ --- + "SurfaceGeometries/Plasma/x" => (; long_name="Cartesian x of plasma-surface point cloud", units="m"), + "SurfaceGeometries/Plasma/y" => (; long_name="Cartesian y of plasma-surface point cloud", units="m"), + "SurfaceGeometries/Plasma/z" => (; long_name="Cartesian z of plasma-surface point cloud", units="m"), + "SurfaceGeometries/Wall/x" => (; long_name="Cartesian x of wall point cloud", units="m"), + "SurfaceGeometries/Wall/y" => (; long_name="Cartesian y of wall point cloud", units="m"), + "SurfaceGeometries/Wall/z" => (; long_name="Cartesian z of wall point cloud", units="m"), +] + +# Euler-Lagrange operator matrices: same wording per letter, Ideal/ and Kinetic/ variants. +const _ELM_IDEAL_LETTERS = [ + ("A", "Euler-Lagrange primitive coefficient matrix A"), + ("B", "Euler-Lagrange primitive coefficient matrix B"), + ("C", "Euler-Lagrange primitive coefficient matrix C"), + ("D", "Euler-Lagrange primitive coefficient matrix D"), + ("E", "Euler-Lagrange primitive coefficient matrix E"), + ("H", "Euler-Lagrange primitive coefficient matrix H"), + ("F", "Euler-Lagrange derived coefficient matrix F"), + ("K", "Euler-Lagrange derived coefficient matrix K"), + ("G", "Euler-Lagrange derived coefficient matrix G"), +] +# The kinetic branch overwrites only A, B, C, K, G and adds f0; D, E, H, F are shared +# unchanged from the ideal set and are not re-emitted. +const _ELM_KINETIC_LETTERS = [ + ("A", "Euler-Lagrange primitive coefficient matrix A"), + ("B", "Euler-Lagrange primitive coefficient matrix B"), + ("C", "Euler-Lagrange primitive coefficient matrix C"), + ("K", "Euler-Lagrange derived coefficient matrix K"), + ("G", "Euler-Lagrange derived coefficient matrix G"), + ("f0", "raw kinetic component matrix f0"), +] +const ELM_H5_ANNOTATIONS = vcat( + ["ForceFreeStates/EulerLagrangeMatrices/psi" => (; long_name="normalized poloidal flux ψ_N grid of the operator matrices", scale="psi")], + ["ForceFreeStates/EulerLagrangeMatrices/Ideal/$l" => + (; long_name="ideal " * d, dims=("psi", "mode_row", "mode_col"), attach=(1 => "ForceFreeStates/EulerLagrangeMatrices/psi",)) for (l, d) in _ELM_IDEAL_LETTERS], + ["ForceFreeStates/EulerLagrangeMatrices/Kinetic/$l" => + (; long_name="kinetic-modified " * d, dims=("psi", "mode_row", "mode_col"), attach=(1 => "ForceFreeStates/EulerLagrangeMatrices/psi",)) for (l, d) in _ELM_KINETIC_LETTERS] +) + +""" + apply_main_h5_metadata!(out_h5, intr) + +Apply the self-describing metadata contract to the datasets written by +`write_outputs_to_HDF5`: `long_name`/`units`/`dims` attributes plus HDF5 Dimension +Scales for the shared coordinate datasets (ψ_N grids, rational-surface ψ). +""" +function apply_main_h5_metadata!(out_h5, intr) + # Scales and attachments are declared in-table via the scale/attach entry fields. + Utilities.HDF5Annotations.annotate!(out_h5, MAIN_H5_ANNOTATIONS) + Utilities.HDF5Annotations.annotate!(out_h5, ELM_H5_ANNOTATIONS) + return out_h5 +end diff --git a/src/InnerLayer/SLAYER/LayerParameters.jl b/src/InnerLayer/SLAYER/LayerParameters.jl index 57ee30763..5dbf9163f 100644 --- a/src/InnerLayer/SLAYER/LayerParameters.jl +++ b/src/InnerLayer/SLAYER/LayerParameters.jl @@ -40,8 +40,8 @@ de-normalization. The parametrization uses `P_perp`, `P_tor`, and | `R0` | Major radius [m] | | `bt` | Toroidal field [T] | | `sval_r` | r-based magnetic shear r_s · (dq/dr) / q (Fitzpatrick convention) | -| `dr_val` | Radial width parameter at surface (input to dc_tmp) | -| `dgeo_val` | Geometric Δ (Shafranov shift factor) | +| `dr_val` | Resistive interchange D_R = E + F + H² (critical-Δ input; auto-derived from GGJ coefficients unless overridden) | +| `dgeo_val` | Connor-Hastie-Helander 2015 Eq. 59 geometric factor (0 unless supplied) | | `eta` | Parallel resistivity entering τ_R = μ₀r_s²/η [Ω·m] | | `d_beta` | Beta-weighted ion length scale c_β · d_i [m] | | `dc_tmp` | Critical-Δ offset from chi_parallel matching | diff --git a/src/KineticForces/Output.jl b/src/KineticForces/Output.jl index 502be6b77..dc613c04f 100644 --- a/src/KineticForces/Output.jl +++ b/src/KineticForces/Output.jl @@ -7,22 +7,25 @@ then write to gpec.h5 in a single pass. """ """ - write_to_hdf5!(h5file::HDF5.File, state::KineticForcesState) + write_to_hdf5!(h5file::HDF5.File, state::KineticForcesState; dVdpsi_spline=nothing) Write all KineticForces results to the "KineticForces" group in gpec.h5. # Arguments - `h5file::HDF5.File`: Open HDF5 file handle - `state::KineticForcesState`: Accumulated computation results +- `dVdpsi_spline`: Optional dV/dψ_N profile interpolant; when given, dV/dψ_N is + written at the quadrature points so the torque density dT/dV = (dT/dψ)/(dV/dψ) + is directly available """ -function write_to_hdf5!(h5file::HDF5.File, state::KineticForcesState) +function write_to_hdf5!(h5file::HDF5.File, state::KineticForcesState; dVdpsi_spline=nothing) g = create_group(h5file, "KineticForces") for (method_name, result) in state.method_results mg = create_group(g, method_name) - mg["nn"] = result.nn + mg["n"] = result.nn # Torque and kinetic energy are the two real physical quantities packed into the - # complex T (Re = T_φ, Im = 2n·δW_k); store each once as a real scalar. + # complex T (Re = T_φ, Im = 2n·δW_k); total_energy stores δW_k = Im(T)/(2n). mg["total_torque"] = real(result.total_torque) mg["total_energy"] = result.total_energy mg["psi_nsteps"] = result.psi_nsteps @@ -35,14 +38,15 @@ function write_to_hdf5!(h5file::HDF5.File, state::KineticForcesState) mg["resonance_psi"] = result.resonance_psis end - # Per-ψ torque profiles from quadrature evaluation points. + # Per-ψ complex torque profiles from quadrature evaluation points: # dT/dψ integrand values and cumulative T(ψ) via trapezoidal integration. if !isempty(result.psi_grid) mg["psi"] = result.psi_grid - mg["dTdpsi_real"] = real.(result.dtdpsi) - mg["dTdpsi_imag"] = imag.(result.dtdpsi) - mg["T_real"] = real.(result.t_cumulative) - mg["T_imag"] = imag.(result.t_cumulative) + mg["dTdpsi"] = result.dtdpsi + mg["T"] = result.t_cumulative + if dVdpsi_spline !== nothing + mg["dVdpsi"] = [dVdpsi_spline(p) for p in result.psi_grid] + end end if !isempty(result.records) @@ -54,12 +58,54 @@ function write_to_hdf5!(h5file::HDF5.File, state::KineticForcesState) for (method_name, mat) in state.kinetic_matrices method_g = haskey(g, method_name) ? g[method_name] : create_group(g, method_name) mat_g = create_group(method_g, "KineticMatrices") - for k in 1:6 - mat_g["matrix_$k"] = mat[:, :, k] + for (letter, k) in _KINETIC_MATRIX_LETTERS + mat_g[letter] = mat[:, :, k] end end + + # Metadata pass: method tokens are data-driven, so annotate each method group. + for method_name in keys(g) + Utilities.HDF5Annotations.annotate!(g[method_name], KF_METHOD_H5_ANNOTATIONS) + end end +# Storage slice → Logan 2015 matrix letter (Eqs 7.30-7.35); see BounceAveraging.jl +# for the Hermitian/full packing details. +const _KINETIC_MATRIX_LETTERS = (("A", 1), ("B", 2), ("C", 3), ("D", 4), ("E", 5), ("H", 6)) + +# Metadata table per KineticForces// group (paths relative to the method group). +# The NTV torque and kinetic energy follow Logan et al. (2013); the six drift-kinetic +# coefficient matrices are Logan 2015 Eqs 7.30-7.35, stored in the energy (δW) +# normalization (torque integrand divided by 2in). +const KF_METHOD_H5_ANNOTATIONS = [ + "n" => (; long_name="toroidal mode number n of this torque calculation"), + "total_torque" => (; long_name="total NTV toroidal torque T_φ", units="N*m"), + "total_energy" => (; long_name="total perturbed kinetic energy δW_k = Im(T)/(2n)", units="J"), + "psi_nsteps" => (; long_name="number of ψ_N quadrature evaluations"), + "panel_psi" => (; long_name="ψ_N panel boundaries of the radial quadrature"), + "resonance_psi" => (; long_name="ψ_N of located kinetic-resonance surfaces"), + "psi" => (; long_name="normalized poloidal flux ψ_N at quadrature evaluation points", scale="psi"), + "dTdpsi" => (; long_name="complex torque density dT/dψ_N = dT_φ/dψ_N + 2i·n·dδW_k/dψ_N; divide by dVdpsi for dT/dV", units="N*m", dims=("psi",), attach=(1 => "psi",)), + "T" => (; long_name="cumulative complex torque T(ψ_N) = T_φ + 2i·n·δW_k (trapezoidal)", units="N*m", dims=("psi",), attach=(1 => "psi",)), + "dVdpsi" => (; long_name="flux-surface volume derivative dV/dψ_N at quadrature points", units="m^3", dims=("psi",), attach=(1 => "psi",)), + "EnergyIntegrals/psi" => (; long_name="ψ_N of each energy-integration record"), + "EnergyIntegrals/lambda" => (; long_name="pitch λ = μB0/E of each record"), + "EnergyIntegrals/ell" => (; long_name="bounce harmonic ℓ of each record"), + "EnergyIntegrals/leff" => (; long_name="effective bounce harmonic ℓ_eff = ℓ + n·q (circulating) or ℓ (trapped)"), + "EnergyIntegrals/torque" => (; long_name="complex torque contribution of the record", units="N*m"), + "EnergyIntegrals/kinetic_energy" => (; long_name="complex kinetic energy contribution of the record", units="J"), + "EnergyIntegrals/trajectory_offsets" => (; long_name="ragged-array offsets: record k spans offsets[k]+1:offsets[k+1] of the *_all arrays"), + "EnergyIntegrals/x_all" => (; long_name="normalized energy x = E/T abscissae of all integration trajectories (concatenated)"), + "EnergyIntegrals/integrand_all" => (; long_name="complex energy-space torque integrand along all trajectories (concatenated)"), + "EnergyIntegrals/integral_all" => (; long_name="complex cumulative energy-space integral along all trajectories (concatenated)"), + "KineticMatrices/A" => (; long_name="Logan 2015 drift-kinetic coefficient matrix A = W_Z†W_Z (energy normalization)", dims=("mode_row", "mode_col")), + "KineticMatrices/B" => (; long_name="Logan 2015 drift-kinetic coefficient matrix B = W_Z†W_X (energy normalization)", dims=("mode_row", "mode_col")), + "KineticMatrices/C" => (; long_name="Logan 2015 drift-kinetic coefficient matrix C = W_Z†W_Y (energy normalization)", dims=("mode_row", "mode_col")), + "KineticMatrices/D" => (; long_name="Logan 2015 drift-kinetic coefficient matrix D = W_X†W_X (energy normalization)", dims=("mode_row", "mode_col")), + "KineticMatrices/E" => (; long_name="Logan 2015 drift-kinetic coefficient matrix E = W_X†W_Y (energy normalization)", dims=("mode_row", "mode_col")), + "KineticMatrices/H" => (; long_name="Logan 2015 drift-kinetic coefficient matrix H = W_Y†W_Y (energy normalization)", dims=("mode_row", "mode_col")) +] + """ write_integration_records!(mg::HDF5.Group, records::Vector{EnergyIntegrationResult}) @@ -78,10 +124,8 @@ function write_integration_records!(mg::HDF5.Group, records::Vector{EnergyIntegr rg["lambda"] = [r.lambda for r in records] rg["ell"] = [r.ell for r in records] rg["leff"] = [r.leff for r in records] - rg["torque_real"] = [real(r.torque) for r in records] - rg["torque_imag"] = [imag(r.torque) for r in records] - rg["kinetic_energy_real"] = [real(r.kinetic_energy) for r in records] - rg["kinetic_energy_imag"] = [imag(r.kinetic_energy) for r in records] + rg["torque"] = [r.torque for r in records] + rg["kinetic_energy"] = [r.kinetic_energy for r in records] # Variable-length trajectories: concatenate all, store offsets for indexing lengths = [length(r.x_trajectory) for r in records] @@ -90,10 +134,8 @@ function write_integration_records!(mg::HDF5.Group, records::Vector{EnergyIntegr if sum(lengths) > 0 rg["x_all"] = vcat([r.x_trajectory for r in records]...) - rg["integrand_real_all"] = vcat([real.(r.integrand_trajectory) for r in records]...) - rg["integrand_imag_all"] = vcat([imag.(r.integrand_trajectory) for r in records]...) - rg["integral_real_all"] = vcat([real.(r.integral_trajectory) for r in records]...) - rg["integral_imag_all"] = vcat([imag.(r.integral_trajectory) for r in records]...) + rg["integrand_all"] = vcat([r.integrand_trajectory for r in records]...) + rg["integral_all"] = vcat([r.integral_trajectory for r in records]...) end end diff --git a/src/PerturbedEquilibrium/Utils.jl b/src/PerturbedEquilibrium/Utils.jl index 9f14bd457..562b0b0f6 100644 --- a/src/PerturbedEquilibrium/Utils.jl +++ b/src/PerturbedEquilibrium/Utils.jl @@ -61,14 +61,14 @@ PerturbedEquilibrium/ ├── forcing_b / forcing_b_root_area / forcing_b_area # control-surface forcing spectrum (b, b̃, b̄) [numpert_total], tesla ├── response_b / response_b_root_area / response_b_area # control-surface response spectrum (b, b̃, b̄) [numpert_total], tesla ├── Response/ -│ ├── psi_n # Radial abscissa ψ_N [npsi] shared by every response profile below +│ ├── psi # Radial abscissa ψ_N [npsi] shared by every response profile below │ ├── xi_psi # Radial displacement ξ^ψ = ξ·∇ψ (ComplexF64 [npsi, mpert]) -│ ├── xi_psi_J # J·ξ^ψ Jacobian-weighted (from gpeq_contra) +│ ├── Jxi_psi # J·ξ^ψ Jacobian-weighted (from gpeq_contra) │ ├── b_psi_area_weighted # b^ψ / ⟨J·|∇ψ|⟩_θ area-normalized (ComplexF64 [npsi, mpert]) │ ├── b_n # Physical normal field b_n (ComplexF64 [npsi, mpert]) │ ├── xi_n # Physical normal displacement xi_n (ComplexF64 [npsi, mpert]) -│ ├── b_theta -│ └── b_zeta +│ ├── Jb_theta +│ └── Jb_zeta ├── ResponseMatrices/ # [numpert_total × numpert_total], root-area-weighted field (b̃) space; R = S·A │ ├── plasma_inductance # Λ̃ = R⁻¹·Λ·R⁻† │ ├── surface_inductance # L̃ = R⁻¹·L·R⁻† @@ -81,17 +81,17 @@ PerturbedEquilibrium/ │ ├── C_resonant_current │ ├── C_island_width_sq │ ├── C_penetrated_area_weighted_field -│ ├── C_delta_prime +│ ├── C_Delta_prime │ ├── resonant_area_weighted_field # [n_rational] applied vector = C̃ · b̃_x (resonant area-weighted field b^r [T]) │ ├── resonant_current │ ├── island_width_sq │ ├── penetrated_area_weighted_field -│ ├── delta_prime +│ ├── Delta_prime │ ├── island_half_width # [n_rational] Float64 │ ├── chirikov_parameter │ ├── rational_psi # [n_rational] surface metadata │ ├── rational_q -│ ├── rational_m_res +│ ├── rational_m │ └── rational_n └── Energies/ ├── vacuum_energy @@ -137,46 +137,46 @@ function write_outputs_to_HDF5( # Response fields (ComplexF64 directly) response_group = haskey(pe_group, "Response") ? pe_group["Response"] : create_group(pe_group, "Response") - !isempty(state.psi_grid) && (response_group["psi_n"] = state.psi_grid) + !isempty(state.psi_grid) && (response_group["psi"] = state.psi_grid) have_xi = !isnothing(state.xi_modes) have_b = have_xi && !isnothing(state.b_modes) response_group["xi_psi"] = have_xi ? state.xi_modes.psi : ComplexF64[] response_group["b_psi_area_weighted"] = have_b ? state.b_modes.b_psi_area_weighted : ComplexF64[] - response_group["b_theta"] = have_b ? state.b_modes.theta : ComplexF64[] - response_group["b_zeta"] = have_b ? state.b_modes.zeta : ComplexF64[] + response_group["Jb_theta"] = have_b ? state.b_modes.theta : ComplexF64[] + response_group["Jb_zeta"] = have_b ? state.b_modes.zeta : ComplexF64[] response_group["b_n"] = !isnothing(state.b_n_modes) ? state.b_n_modes : ComplexF64[] response_group["xi_n"] = !isnothing(state.xi_n_modes) ? state.xi_n_modes : ComplexF64[] # Clebsch displacements for PENTRC (matches Fortran gpout_xclebsch) if have_xi - response_group["clebsch_psi"] = state.xi_modes.clebsch_psi - response_group["clebsch_psi1"] = state.xi_modes.clebsch_psi1 - response_group["clebsch_alpha"] = state.xi_modes.clebsch_alpha + response_group["xi_clebsch_psi"] = state.xi_modes.clebsch_psi + response_group["dxi_clebsch_psidpsi"] = state.xi_modes.clebsch_psi1 + response_group["xi_clebsch_alpha"] = state.xi_modes.clebsch_alpha end # Contravariant displacement (from gpeq_contra, all J-weighted) if have_xi - response_group["xi_psi_J"] = state.xi_modes.psi_J - response_group["xi_theta"] = state.xi_modes.theta - response_group["xi_zeta"] = state.xi_modes.zeta + response_group["Jxi_psi"] = state.xi_modes.psi_J + response_group["Jxi_theta"] = state.xi_modes.theta + response_group["Jxi_zeta"] = state.xi_modes.zeta end # Covariant components (from gpeq_cova) if have_xi - response_group["xi_cova_psi"] = state.xi_modes.cova_psi - response_group["xi_cova_theta"] = state.xi_modes.cova_theta - response_group["xi_cova_zeta"] = state.xi_modes.cova_zeta + response_group["xi_cov_psi"] = state.xi_modes.cova_psi + response_group["xi_cov_theta"] = state.xi_modes.cova_theta + response_group["xi_cov_zeta"] = state.xi_modes.cova_zeta end if have_xi - response_group["xi_theta_reg"] = state.xi_modes.theta_reg - response_group["xi_zeta_reg"] = state.xi_modes.zeta_reg + response_group["Jxi_theta_reg"] = state.xi_modes.theta_reg + response_group["Jxi_zeta_reg"] = state.xi_modes.zeta_reg end if have_b - response_group["b_theta_reg"] = state.b_modes.theta_reg - response_group["b_zeta_reg"] = state.b_modes.zeta_reg - response_group["b_cova_psi"] = state.b_modes.cova_psi - response_group["b_cova_theta"] = state.b_modes.cova_theta - response_group["b_cova_zeta"] = state.b_modes.cova_zeta + response_group["Jb_theta_reg"] = state.b_modes.theta_reg + response_group["Jb_zeta_reg"] = state.b_modes.zeta_reg + response_group["b_cov_psi"] = state.b_modes.cova_psi + response_group["b_cov_theta"] = state.b_modes.cova_theta + response_group["b_cov_zeta"] = state.b_modes.cova_zeta end # R,Z,φ cylindrical components in mode-space (from gpeq_rzphi) @@ -199,14 +199,14 @@ function write_outputs_to_HDF5( !isempty(state.C_resonant_current) && (coupling_group["C_resonant_current"] = state.C_resonant_current) !isempty(state.C_island_width_sq) && (coupling_group["C_island_width_sq"] = state.C_island_width_sq) !isempty(state.C_penetrated_area_weighted_field) && (coupling_group["C_penetrated_area_weighted_field"] = state.C_penetrated_area_weighted_field) - !isempty(state.C_delta_prime) && (coupling_group["C_delta_prime"] = state.C_delta_prime) + !isempty(state.C_delta_prime) && (coupling_group["C_Delta_prime"] = state.C_delta_prime) # Applied resonant vectors [n_rational] !isempty(state.resonant_area_weighted_field) && (coupling_group["resonant_area_weighted_field"] = state.resonant_area_weighted_field) !isempty(state.resonant_current) && (coupling_group["resonant_current"] = state.resonant_current) !isempty(state.island_width_sq) && (coupling_group["island_width_sq"] = state.island_width_sq) !isempty(state.penetrated_area_weighted_field) && (coupling_group["penetrated_area_weighted_field"] = state.penetrated_area_weighted_field) - !isempty(state.delta_prime) && (coupling_group["delta_prime"] = state.delta_prime) + !isempty(state.delta_prime) && (coupling_group["Delta_prime"] = state.delta_prime) !isempty(state.forcing_solution_weights) && (coupling_group["forcing_solution_weights"] = state.forcing_solution_weights) !isempty(state.rational_area) && (coupling_group["rational_area"] = state.rational_area) !isempty(state.island_half_width) && (coupling_group["island_half_width"] = state.island_half_width) @@ -215,7 +215,7 @@ function write_outputs_to_HDF5( # Metadata [n_rational] !isempty(state.rational_psi) && (coupling_group["rational_psi"] = state.rational_psi) !isempty(state.rational_q) && (coupling_group["rational_q"] = state.rational_q) - !isempty(state.rational_m_res) && (coupling_group["rational_m_res"] = state.rational_m_res) + !isempty(state.rational_m_res) && (coupling_group["rational_m"] = state.rational_m_res) !isempty(state.rational_n) && (coupling_group["rational_n"] = state.rational_n) # Energies @@ -224,5 +224,133 @@ function write_outputs_to_HDF5( energy_group["surface_energy"] = state.surface_energy energy_group["plasma_energy"] = state.plasma_energy energy_group["toroidal_torque"] = state.toroidal_torque + + annotate_pe!(pe_group) end end + +# Metadata tables for the PerturbedEquilibrium group, applied post-write (paths are +# relative to the PerturbedEquilibrium group). Field-representation naming follows +# docs/src/conventions.md (Pharr 2026): b = bare, b̄ = area-weighted, b̃ = +# root-area-weighted; all in tesla. +const PE_H5_ANNOTATIONS = [ + "ForcingModes/n" => (; long_name="toroidal mode number of each forcing mode"), + "ForcingModes/m" => (; long_name="poloidal mode number of each forcing mode"), + "ForcingModes/amplitude" => (; long_name="complex forcing amplitude of each mode", units="T"), + "forcing_b" => (; long_name="control-surface forcing spectrum, bare normal field b", units="T", dims=("mode",)), + "forcing_b_root_area" => (; long_name="control-surface forcing spectrum, root-area-weighted field b̃ (coordinate-invariant)", units="T", dims=("mode",)), + "forcing_b_area" => (; long_name="control-surface forcing spectrum, area-weighted field b̄ (Φ = A·b̄)", units="T", dims=("mode",)), + "response_b" => (; long_name="control-surface response spectrum, bare normal field b", units="T", dims=("mode",)), + "response_b_root_area" => (; long_name="control-surface response spectrum, root-area-weighted field b̃ (coordinate-invariant)", units="T", dims=("mode",)), + "response_b_area" => (; long_name="control-surface response spectrum, area-weighted field b̄ (Φ = A·b̄)", units="T", dims=("mode",)), + "ResponseMatrices/plasma_inductance" => (; long_name="plasma inductance Λ̃ in root-area-weighted field space", dims=("mode_row", "mode_col")), + "ResponseMatrices/surface_inductance" => (; long_name="surface inductance L̃ in root-area-weighted field space", dims=("mode_row", "mode_col")), + "ResponseMatrices/permeability" => (; long_name="permeability P̃ = Λ̃·L̃⁻¹ in root-area-weighted field space", dims=("mode_row", "mode_col")), + "ResponseMatrices/reluctance" => (; long_name="reluctance ϱ̃ in root-area-weighted field space", dims=("mode_row", "mode_col")), + "ResponseMatrices/rootarea_to_area_weight_operator" => (; long_name="operator S = Σ/√A at ψ_lim; b̄ = S·b̃", dims=("mode_row", "mode_col")), + "ResponseMatrices/surface_area" => (; long_name="control-surface scalar area A = ∮J|∇ψ|dθ; Φ = A·b̄", units="m^2"), + "Response/psi" => (; long_name="normalized poloidal flux ψ_N grid shared by the response profiles", scale="psi"), + "Response/xi_psi" => (; long_name="contravariant radial displacement ξ^ψ = ξ·∇ψ_N", dims=("psi", "mode"), attach=(1 => "Response/psi",)), + "Response/Jxi_psi" => (; long_name="Jacobian-weighted contravariant radial displacement J·ξ^ψ", units="m^3", dims=("psi", "mode"), attach=(1 => "Response/psi",)), + "Response/Jxi_theta" => (; long_name="Jacobian-weighted contravariant poloidal displacement J·ξ^θ", units="m^3", dims=("psi", "mode"), attach=(1 => "Response/psi",)), + "Response/Jxi_zeta" => (; long_name="Jacobian-weighted contravariant toroidal displacement J·ξ^ζ", units="m^3", dims=("psi", "mode"), attach=(1 => "Response/psi",)), + "Response/Jxi_theta_reg" => + (; long_name="regularized Jacobian-weighted contravariant poloidal displacement J·ξ^θ", units="m^3", dims=("psi", "mode"), attach=(1 => "Response/psi",)), + "Response/Jxi_zeta_reg" => + (; long_name="regularized Jacobian-weighted contravariant toroidal displacement J·ξ^ζ", units="m^3", dims=("psi", "mode"), attach=(1 => "Response/psi",)), + "Response/xi_cov_psi" => (; long_name="covariant radial displacement ξ_ψ = ξ·e_ψ", units="m^2", dims=("psi", "mode"), attach=(1 => "Response/psi",)), + "Response/xi_cov_theta" => (; long_name="covariant poloidal displacement ξ_θ = ξ·e_θ", units="m^2", dims=("psi", "mode"), attach=(1 => "Response/psi",)), + "Response/xi_cov_zeta" => (; long_name="covariant toroidal displacement ξ_ζ = ξ·e_ζ", units="m^2", dims=("psi", "mode"), attach=(1 => "Response/psi",)), + "Response/xi_clebsch_psi" => (; long_name="Clebsch displacement component ξ^ψ (PENTRC input, gpout_xclebsch convention)", dims=("psi", "mode"), attach=(1 => "Response/psi",)), + "Response/dxi_clebsch_psidpsi" => + (; long_name="regularized ψ_N derivative of ξ^ψ (× singfac²/(singfac²+reg_spot²))", dims=("psi", "mode"), attach=(1 => "Response/psi",)), + "Response/xi_clebsch_alpha" => + (; long_name="Clebsch displacement component ξ^α/χ₁ (PENTRC input, gpout_xclebsch convention)", dims=("psi", "mode"), attach=(1 => "Response/psi",)), + "Response/xi_n" => (; long_name="physical normal displacement ξ_n", units="m", dims=("psi", "mode"), attach=(1 => "Response/psi",)), + "Response/xi_R" => (; long_name="cylindrical displacement component ξ_R (mode space)", units="m", dims=("psi", "mode"), attach=(1 => "Response/psi",)), + "Response/xi_Z" => (; long_name="cylindrical displacement component ξ_Z (mode space)", units="m", dims=("psi", "mode"), attach=(1 => "Response/psi",)), + "Response/xi_phi" => (; long_name="cylindrical displacement component ξ_φ (mode space)", units="m", dims=("psi", "mode"), attach=(1 => "Response/psi",)), + "Response/b_psi_area_weighted" => + (; long_name="area-weighted normal field b̄ = (J·b^ψ)_m/A (Φ = A·b̄; same representation as response_b_area)", units="T", + dims=("psi", "mode"), attach=(1 => "Response/psi",)), + "Response/b_n" => (; long_name="physical normal field b_n", units="T", dims=("psi", "mode"), attach=(1 => "Response/psi",)), + "Response/Jb_theta" => + (; long_name="Jacobian-weighted contravariant poloidal field J·b^θ (flux-like density, not a tesla field amplitude)", units="T*m^2", + dims=("psi", "mode"), attach=(1 => "Response/psi",)), + "Response/Jb_zeta" => + (; long_name="Jacobian-weighted contravariant toroidal field J·b^ζ (flux-like density, not a tesla field amplitude)", units="T*m^2", + dims=("psi", "mode"), attach=(1 => "Response/psi",)), + "Response/Jb_theta_reg" => + (; long_name="regularized Jacobian-weighted contravariant poloidal field J·b^θ (flux-like density, not a tesla field amplitude)", units="T*m^2", dims=("psi", "mode"), + attach=(1 => "Response/psi",)), + "Response/Jb_zeta_reg" => + (; long_name="regularized Jacobian-weighted contravariant toroidal field J·b^ζ (flux-like density, not a tesla field amplitude)", units="T*m^2", + dims=("psi", "mode"), attach=(1 => "Response/psi",)), + "Response/b_cov_psi" => (; long_name="covariant radial field b_ψ = b·e_ψ", units="T*m", dims=("psi", "mode"), attach=(1 => "Response/psi",)), + "Response/b_cov_theta" => (; long_name="covariant poloidal field b_θ = b·e_θ", units="T*m", dims=("psi", "mode"), attach=(1 => "Response/psi",)), + "Response/b_cov_zeta" => (; long_name="covariant toroidal field b_ζ = b·e_ζ", units="T*m", dims=("psi", "mode"), attach=(1 => "Response/psi",)), + "Response/b_R" => (; long_name="cylindrical field component b_R (mode space)", units="T", dims=("psi", "mode"), attach=(1 => "Response/psi",)), + "Response/b_Z" => (; long_name="cylindrical field component b_Z (mode space)", units="T", dims=("psi", "mode"), attach=(1 => "Response/psi",)), + "Response/b_phi" => (; long_name="cylindrical field component b_φ (mode space)", units="T", dims=("psi", "mode"), attach=(1 => "Response/psi",)), + "SingularCoupling/C_resonant_area_weighted_field" => + (; long_name="coupling matrix: applied b̃ → resonant area-weighted field b̄^r = Φ^r/A^r", dims=("surface", "mode"), + attach=(1 => "SingularCoupling/rational_psi", 1 => "SingularCoupling/rational_q")), + "SingularCoupling/C_resonant_current" => + (; long_name="coupling matrix: applied b̃ → pitch-resonant current", units="A/T", dims=("surface", "mode"), + attach=(1 => "SingularCoupling/rational_psi", 1 => "SingularCoupling/rational_q")), + "SingularCoupling/C_island_width_sq" => + (; long_name="coupling matrix: applied b̃ → squared island half-width", units="1/T", dims=("surface", "mode"), + attach=(1 => "SingularCoupling/rational_psi", 1 => "SingularCoupling/rational_q")), + "SingularCoupling/C_penetrated_area_weighted_field" => + (; long_name="coupling matrix: applied b̃ → penetrated area-weighted field", dims=("surface", "mode"), + attach=(1 => "SingularCoupling/rational_psi", 1 => "SingularCoupling/rational_q")), + "SingularCoupling/C_Delta_prime" => + (; long_name="coupling matrix: applied b̃ → forcing-driven Δ'", units="1/T", dims=("surface", "mode"), + attach=(1 => "SingularCoupling/rational_psi", 1 => "SingularCoupling/rational_q")), + "SingularCoupling/resonant_area_weighted_field" => + (; long_name="resonant area-weighted field b̄^r = Φ^r/A^r per rational surface (coordinate-invariant)", units="T", dims=("surface",), + attach=(1 => "SingularCoupling/rational_psi", 1 => "SingularCoupling/rational_q")), + "SingularCoupling/resonant_current" => + (; long_name="pitch-resonant current per rational surface", units="A", dims=("surface",), + attach=(1 => "SingularCoupling/rational_psi", 1 => "SingularCoupling/rational_q")), + "SingularCoupling/island_width_sq" => + (; long_name="squared island half-width per rational surface (in ψ_N²)", dims=("surface",), + attach=(1 => "SingularCoupling/rational_psi", 1 => "SingularCoupling/rational_q")), + "SingularCoupling/penetrated_area_weighted_field" => + (; long_name="penetrated area-weighted field per rational surface", units="T", dims=("surface",), + attach=(1 => "SingularCoupling/rational_psi", 1 => "SingularCoupling/rational_q")), + "SingularCoupling/Delta_prime" => + (; long_name="forcing-driven tearing Δ' per rational surface (Riccati; response to applied forcing)", dims=("surface",), + attach=(1 => "SingularCoupling/rational_psi", 1 => "SingularCoupling/rational_q")), + "SingularCoupling/forcing_solution_weights" => + (; long_name="weights of the forcing solutions in the singular-coupling decomposition", dims=("surface",), + attach=(1 => "SingularCoupling/rational_psi", 1 => "SingularCoupling/rational_q")), + "SingularCoupling/rational_area" => + (; long_name="scalar surface area A^r of each rational surface", units="m^2", dims=("surface",), + attach=(1 => "SingularCoupling/rational_psi", 1 => "SingularCoupling/rational_q")), + "SingularCoupling/island_half_width" => + (; long_name="island half-width per rational surface (in ψ_N)", dims=("surface",), attach=(1 => "SingularCoupling/rational_psi", 1 => "SingularCoupling/rational_q")), + "SingularCoupling/chirikov_parameter" => + (; long_name="Chirikov overlap parameter: island half-width / half-distance to the neighbouring rational surface", dims=("surface",), + attach=(1 => "SingularCoupling/rational_psi", 1 => "SingularCoupling/rational_q")), + "SingularCoupling/rational_psi" => (; long_name="normalized poloidal flux ψ_N of each rational surface", scale="psi_rational"), + "SingularCoupling/rational_q" => + (; long_name="safety factor q = m/n at each rational surface", dims=("surface",), scale="q_rational", attach=(1 => "SingularCoupling/rational_psi",)), + "SingularCoupling/rational_m" => + (; long_name="resonant poloidal mode number m at each rational surface", dims=("surface",), + attach=(1 => "SingularCoupling/rational_psi", 1 => "SingularCoupling/rational_q")), + "SingularCoupling/rational_n" => + (; long_name="resonant toroidal mode number n at each rational surface", + dims=("surface",), attach=(1 => "SingularCoupling/rational_psi", 1 => "SingularCoupling/rational_q")), + "Energies/vacuum_energy" => (; long_name="perturbed vacuum energy", units="J"), + "Energies/surface_energy" => (; long_name="perturbed surface energy", units="J"), + "Energies/plasma_energy" => (; long_name="perturbed plasma energy", units="J"), + "Energies/toroidal_torque" => (; long_name="net toroidal torque on the plasma", units="N*m") +] + +# Attach long_name/units/dims + dimension scales (declared in-table) to the +# PerturbedEquilibrium group. +function annotate_pe!(pe_group) + Utilities.HDF5Annotations.annotate!(pe_group, PE_H5_ANNOTATIONS) + return nothing +end diff --git a/src/Tearing/Runner/HDF5Output.jl b/src/Tearing/Runner/HDF5Output.jl index 2c1ee6611..9bbbf3292 100644 --- a/src/Tearing/Runner/HDF5Output.jl +++ b/src/Tearing/Runner/HDF5Output.jl @@ -12,7 +12,7 @@ # ├── PerSurface/ -- struct-of-arrays for SLAYERParameters fields # │ ├── psi, q, q1, ... # │ └── ... -# ├── Roots/ -- Q_root (real, imag), omega_Hz, gamma_Hz +# ├── Roots/ -- complex Q_root, omega_Hz, gamma_Hz # ├── Diagnostics/ -- ValidRoots, Poles, FilteredRoots # │ (flat-plus-offsets ragged encoding) # └── Scan/ -- optional: full Q/Δ scan data @@ -34,8 +34,14 @@ function write_slayer_hdf5!(parent::Union{HDF5.File,HDF5.Group}, end g = create_group(parent, "Tearing") g["enabled"] = Int(result.enabled) + # Which inner-layer model produced PerSurface/: the SLAYER and GGJ branches write + # disjoint field sets, so readers must not have to infer it from the schema. + attrs(g)["layer_model"] = result.enabled ? _layer_model_token(eltype(result.params)) : "none" - result.enabled || return g # nothing else to write + if !result.enabled # nothing else to write + _annotate_tearing!(g) + return g + end _write_per_surface!(g, result.params, result.dp_matrix) _write_roots!(g, result) @@ -44,32 +50,137 @@ function write_slayer_hdf5!(parent::Union{HDF5.File,HDF5.Group}, if result.control.store_scan && !isempty(result.scan_data) _write_scan_data!(g, result) end + _annotate_tearing!(g) return g end +# Token recorded in the Tearing group's layer_model attribute; keyed by the +# per-surface parameter type since the SLAYER and GGJ branches write disjoint fields. +_layer_model_token(::Type{SLAYERParameters}) = "slayer" +_layer_model_token(::Type{GGJParameters}) = "ggj" + +# Metadata table for the Tearing group (paths relative to it); ragged Diagnostics +# subgroups and Scan/Surface_ groups are annotated by iteration below. +const TEARING_H5_ANNOTATIONS = [ + "enabled" => (; long_name="flag: SLAYER/tearing stage ran (1) or was disabled (0)"), + "PerSurface/rational_index" => (; long_name="rational-surface index of each row", dims=("surface",)), + "PerSurface/m" => (; long_name="resonant poloidal mode number m per surface", dims=("surface",)), + "PerSurface/n" => (; long_name="resonant toroidal mode number n per surface", dims=("surface",)), + "PerSurface/tau" => (; long_name="temperature ratio τ = T_i/T_e per surface", dims=("surface",)), + "PerSurface/lu" => (; long_name="Lundquist number S per surface", dims=("surface",)), + "PerSurface/c_beta" => (; long_name="compressibility factor c_β = √(β_local/(1+β_local)) per surface", dims=("surface",)), + "PerSurface/D_norm" => (; long_name="Fitzpatrick normalized ion-sound/drift scale D = (d_β/r_s)·S^(1/3)·√(τ/(1+τ)) per surface", dims=("surface",)), + "PerSurface/P_perp" => (; long_name="perpendicular magnetic Prandtl number per surface", dims=("surface",)), + "PerSurface/P_tor" => (; long_name="toroidal (momentum) magnetic Prandtl number per surface", dims=("surface",)), + "PerSurface/Q_e" => (; long_name="normalized electron diamagnetic frequency Q_e per surface", dims=("surface",)), + "PerSurface/Q_i" => (; long_name="normalized ion diamagnetic frequency Q_i per surface", dims=("surface",)), + "PerSurface/iota_e" => (; long_name="electron fraction ι_e = Q_e/(Q_e − Q_i) per surface", dims=("surface",)), + "PerSurface/tau_k" => + (; long_name="Q-normalization time S^(1/3)·τ_H per surface (Q = τ_k·ω; diamagnetic inputs Q_e, Q_i carry the opposite sign by convention)", units="s", dims=("surface",)), + "PerSurface/tau_R" => (; long_name="resistive diffusion time τ_R = μ₀r_s²/η per surface", units="s", dims=("surface",)), + "PerSurface/Delta_prime_norm" => (; long_name="Δ'-normalization factor S^(1/3)/r_s per surface", units="1/m", dims=("surface",)), + "PerSurface/rs" => (; long_name="minor radius of each rational surface", units="m", dims=("surface",)), + "PerSurface/R0" => (; long_name="major radius", units="m", dims=("surface",)), + "PerSurface/bt" => (; long_name="toroidal field", units="T", dims=("surface",)), + "PerSurface/sval_r" => (; long_name="r-based magnetic shear r_s·(dq/dr)/q (Fitzpatrick convention)", dims=("surface",)), + "PerSurface/D_R" => + (; long_name="resistive interchange D_R = E + F + H² for the critical-Δ formula (auto-derived from GGJ coefficients unless overridden)", dims=("surface",)), + "PerSurface/D_geo" => (; long_name="Connor-Hastie-Helander 2015 Eq. 59 geometric factor (0 unless supplied)", dims=("surface",)), + "PerSurface/eta" => (; long_name="parallel resistivity at each surface", units="Ohm*m", dims=("surface",)), + "PerSurface/d_beta" => (; long_name="β-weighted ion drift scale d_β", units="m", dims=("surface",)), + "PerSurface/D_c_offset" => (; long_name="critical-Δ offset from χ_∥/χ_⊥ matching (Connor-Hastie-Helander 2015 Eq. 59)", dims=("surface",)), + "PerSurface/D_c_type" => (; long_name="per-surface D_c prescription label", dims=("surface",)), + "PerSurface/E" => (; long_name="Glasser-Greene-Johnson coefficient E per surface", dims=("surface",)), + "PerSurface/F" => (; long_name="Glasser-Greene-Johnson coefficient F per surface", dims=("surface",)), + "PerSurface/G" => (; long_name="Glasser-Greene-Johnson coefficient G per surface", dims=("surface",)), + "PerSurface/H" => (; long_name="Glasser-Greene-Johnson coefficient H per surface", dims=("surface",)), + "PerSurface/K" => (; long_name="Glasser-Greene-Johnson coefficient K per surface", dims=("surface",)), + "PerSurface/M" => (; long_name="Glasser-Greene-Johnson coefficient M per surface", dims=("surface",)), + "PerSurface/tau_A" => (; long_name="Alfvén time τ_A per surface (GGJ layer parameters)", units="s", dims=("surface",)), + "PerSurface/dVdpsi" => (; long_name="dV/dψ_N at each surface", units="m^3", dims=("surface",)), + "PerSurface/Delta_prime_matrix" => (; long_name="full complex Δ' matrix coupling the rational surfaces", dims=("surface_row", "surface_col")), + "Roots/Q_root" => (; long_name="complex dispersion-root normalized frequency Q (NaN = no root)", dims=("surface",)), + "Roots/omega" => + (; long_name="mode rotation angular frequency ω = Re(Q)/τ_k of each root", units="rad/s", dims=("surface",)), + "Roots/gamma" => + (; long_name="growth rate γ = Im(Q)/τ_k of each root, positive = unstable (an e-folding rate, so no 2π distinction applies)", units="1/s", dims=("surface",)), + "Roots/no_root" => (; long_name="flag: no usable dispersion root found (Q_root is NaN, ω/γ are placeholders)", dims=("surface",)), + "LayerWidths/rational_index" => (; long_name="rational-surface index of each row", dims=("surface",)), + "LayerWidths/m" => (; long_name="resonant poloidal mode number m per surface", dims=("surface",)), + "LayerWidths/n" => (; long_name="resonant toroidal mode number n per surface", dims=("surface",)), + "LayerWidths/delta_s_over_d_beta" => (; long_name="complex dimensionless layer thickness δ_s/d_β", dims=("surface",)), + "LayerWidths/delta_s" => (; long_name="complex resistive layer thickness δ_s (Riccati)", dims=("surface",)), + "LayerWidths/delta_s_abs" => (; long_name="physical resistive layer thickness |δ_s|", units="m", dims=("surface",)), + "LayerWidths/d_beta" => (; long_name="β-weighted ion drift scale d_β", units="m", dims=("surface",)) +] + +const TEARING_RAGGED_H5_ANNOTATIONS = [ + "flat" => (; long_name="concatenated complex entries (rows delimited by offsets)"), + "offsets" => (; long_name="ragged-array offsets: row k spans offsets[k]+1:offsets[k+1]") +] + +const TEARING_SCAN_H5_ANNOTATIONS = [ + "kind" => (; long_name="scan kind: brute_force or amr"), + "Q" => (; long_name="sampled complex normalized frequency Q"), + "Delta" => (; long_name="complex inner-layer matching Δ(Q)"), + "re_axis" => (; long_name="Re(Q) axis of the brute-force scan grid"), + "im_axis" => (; long_name="Im(Q) axis of the brute-force scan grid"), + "n_cells" => (; long_name="number of AMR cells sampled"), + "truncated" => (; long_name="flag: AMR refinement stopped at the cell cap") +] + +# Attach long_name/units/dims to everything write_slayer_hdf5! wrote. +function _annotate_tearing!(g) + ann = Utilities.HDF5Annotations + ann.annotate!(g, TEARING_H5_ANNOTATIONS) + if haskey(g, "Diagnostics") + for sub in keys(g["Diagnostics"]) + ann.annotate!(g["Diagnostics"][sub], TEARING_RAGGED_H5_ANNOTATIONS) + end + end + if haskey(g, "Scan") + for sub in keys(g["Scan"]) + sg = g["Scan"][sub] + ann.annotate!(sg, TEARING_SCAN_H5_ANNOTATIONS) + # Brute-force scans store 2-D (re, im) grids; AMR stores flat samples. + if haskey(sg, "Q") && ndims(sg["Q"]) == 2 + for a in ("Q", "Delta") + haskey(sg, a) && (attrs(sg[a])["dims"] = "(re_axis, im_axis)") + end + end + end + end + return nothing +end + # ---------- per-surface layer parameters ---------- function _write_per_surface!(g, params::AbstractVector{SLAYERParameters}, dp_matrix::Matrix{ComplexF64}) ps = create_group(g, "PerSurface") # Scalar struct-of-arrays for all Float64 / Int fields - for fname in (:ising, :m, :n) + ps["rational_index"] = Int[p.ising for p in params] + for fname in (:m, :n) ps[String(fname)] = Int[getfield(p, fname) for p in params] end for fname in (:tau, :lu, :c_beta, :D_norm, :P_perp, :P_tor, :Q_e, :Q_i, :iota_e, - :tauk, :tau_r, :delta_n, - :rs, :R0, :bt, :sval_r, :dr_val, :dgeo_val, - :eta, :d_beta, :dc_tmp) + :rs, :R0, :bt, :sval_r, + :eta, :d_beta) ps[String(fname)] = Float64[getfield(p, fname) for p in params] end + # Literature/physics dataset names where the struct fields kept legacy spellings. + ps["tau_k"] = Float64[p.tauk for p in params] + ps["tau_R"] = Float64[p.tau_r for p in params] + ps["Delta_prime_norm"] = Float64[p.delta_n for p in params] + ps["D_c_offset"] = Float64[p.dc_tmp for p in params] + # Physics names for the interchange/geometric inputs (struct fields keep *_val). + ps["D_R"] = Float64[p.dr_val for p in params] + ps["D_geo"] = Float64[p.dgeo_val for p in params] # Store dc_type per-surface as string array - ps["dc_type"] = String[String(p.dc_type) for p in params] + ps["D_c_type"] = String[String(p.dc_type) for p in params] - # Full Δ' matrix, split real/imag - dp = create_group(ps, "DpMatrix") - dp["real"] = real.(dp_matrix) - dp["imag"] = imag.(dp_matrix) + ps["Delta_prime_matrix"] = dp_matrix return nothing end @@ -78,23 +189,23 @@ end function _write_per_surface!(g, params::AbstractVector{GGJParameters}, dp_matrix::Matrix{ComplexF64}) ps = create_group(g, "PerSurface") - ps["ising"] = Int[p.ising for p in params] - for fname in (:E, :F, :G, :H, :K, :M, :taua, :taur, :v1) + ps["rational_index"] = Int[p.ising for p in params] + for fname in (:E, :F, :G, :H, :K, :M) ps[String(fname)] = Float64[getfield(p, fname) for p in params] end - dp = create_group(ps, "DpMatrix") - dp["real"] = real.(dp_matrix) - dp["imag"] = imag.(dp_matrix) + ps["tau_A"] = Float64[p.taua for p in params] + ps["tau_R"] = Float64[p.taur for p in params] + ps["dVdpsi"] = Float64[p.v1 for p in params] + ps["Delta_prime_matrix"] = dp_matrix return nothing end # ---------- eigenvalue roots ---------- function _write_roots!(g, r::SLAYERResult) roots = create_group(g, "Roots") - roots["Q_root_real"] = real.(r.Q_root) - roots["Q_root_imag"] = imag.(r.Q_root) - roots["omega_Hz"] = r.omega_Hz - roots["gamma_Hz"] = r.gamma_Hz + roots["Q_root"] = r.Q_root + roots["omega"] = r.omega_Hz + roots["gamma"] = r.gamma_Hz # `no_root[k] == 1` flags entries where the extraction found NO usable # root (Q_root is NaN; omega_Hz/gamma_Hz are 0 placeholders, not a true # γ≈0 result). Aligned element-wise with Q_root/omega_Hz/gamma_Hz. @@ -108,16 +219,15 @@ end # ---------- resistive layer thickness (del_s Riccati) ---------- function _write_layer_widths!(g, widths::Vector{LayerWidths}) lw = create_group(g, "LayerWidths") - for fname in (:ising, :m, :n) + lw["rational_index"] = Int[w.ising for w in widths] + for fname in (:m, :n) lw[String(fname)] = Int[getfield(w, fname) for w in widths] end - # Dimensionless del_s/d_beta and the complex layer thickness, split re/im. - lw["dels_db_real"] = Float64[real(w.dels_db) for w in widths] - lw["dels_db_imag"] = Float64[imag(w.dels_db) for w in widths] - lw["delta_s_real"] = Float64[real(w.delta_s) for w in widths] - lw["delta_s_imag"] = Float64[imag(w.delta_s) for w in widths] + # Dimensionless del_s/d_beta and the complex layer thickness. + lw["delta_s_over_d_beta"] = ComplexF64[w.dels_db for w in widths] + lw["delta_s"] = ComplexF64[w.delta_s for w in widths] # Physical thickness [m] and the β-weighted ion drift scale [m]. - lw["delta_s_m"] = Float64[w.delta_s_m for w in widths] + lw["delta_s_abs"] = Float64[w.delta_s_m for w in widths] lw["d_beta"] = Float64[w.d_beta for w in widths] return nothing end @@ -141,22 +251,19 @@ function _write_diagnostics!(g, r::SLAYERResult) return nothing end -# Write a ragged vector-of-vectors of ComplexF64 as (flat_re, flat_im, -# offsets) — `offsets[k+1] - offsets[k]` is the length of row `k`. This -# avoids HDF5 VLEN types, which have patchy cross-language support. +# Write a ragged vector-of-vectors of ComplexF64 as (flat, offsets) — +# `offsets[k+1] - offsets[k]` is the length of row `k`. This avoids HDF5 VLEN +# types, which have patchy cross-language support. function _write_ragged_complex!(parent, name::String, data::Vector{Vector{ComplexF64}}) g = create_group(parent, name) - flat_re = Float64[] - flat_im = Float64[] + flat = ComplexF64[] offsets = Int[0] for v in data - append!(flat_re, real.(v)) - append!(flat_im, imag.(v)) + append!(flat, v) push!(offsets, offsets[end] + length(v)) end - g["flat_real"] = flat_re - g["flat_imag"] = flat_im + g["flat"] = flat g["offsets"] = offsets return nothing end @@ -173,10 +280,8 @@ end function _write_single_scan!(g, data::ScanResult) g["kind"] = "brute_force" - g["Q_real"] = real.(data.Q) - g["Q_imag"] = imag.(data.Q) - g["Delta_real"] = real.(data.Δ) - g["Delta_imag"] = imag.(data.Δ) + g["Q"] = data.Q + g["Delta"] = data.Δ g["re_axis"] = data.re_axis g["im_axis"] = data.im_axis return nothing @@ -184,10 +289,8 @@ end function _write_single_scan!(g, data::AMRResult) g["kind"] = "amr" - g["Q_real"] = real.(data.Q) - g["Q_imag"] = imag.(data.Q) - g["Delta_real"] = real.(data.Δ) - g["Delta_imag"] = imag.(data.Δ) + g["Q"] = data.Q + g["Delta"] = data.Δ g["n_cells"] = length(data.cells) g["truncated"] = Int(data.truncated) return nothing diff --git a/src/Utilities/HDF5Annotations.jl b/src/Utilities/HDF5Annotations.jl new file mode 100644 index 000000000..1fff44e62 --- /dev/null +++ b/src/Utilities/HDF5Annotations.jl @@ -0,0 +1,122 @@ +""" + HDF5Annotations + +Self-describing metadata for `gpec.h5` (the contract in +`docs/development/hdf5-conventions.md`): every dataset carries a `long_name` and +`units` attribute, array datasets carry a `dims` axis-name attribute, and coordinate +datasets are marked as HDF5 Dimension Scales (netCDF-4 coordinate variables) attached +to the arrays that share their axis, so h5py/xarray/HDFView read the file unaided. + +Writers stay table-driven: each writer keeps a table of `path => (; long_name, units, +dims)` entries next to it and calls [`annotate!`](@ref) once after its datasets are +written. Paths absent from the file are skipped silently (many writes are conditional). +""" +module HDF5Annotations + +using HDF5 +using Dates + +export annotate!, make_scale!, attach_scale!, write_root_attrs! + +""" + annotate!(parent, table) + +Apply a metadata table to datasets under `parent` (an open `HDF5.File` or group). +`table` iterates `path => meta` pairs where `meta` is a NamedTuple with fields +`long_name` (required), `units` (default `"1"` = dimensionless), and optionally: + + - `dims` — tuple of axis names in Julia (column-major) order, axis 1 first, stored + as the greppable string attribute `dims = "(psi, m)"`. + - `scale` — mark this dataset as an HDF5 Dimension Scale with the given name + (netCDF-4 coordinate variable). + - `attach` — tuple of `axis => scale_path` pairs attaching declared scales to this + dataset's Julia axes; `scale_path` is relative to `parent` and its dimension + label comes from that entry's `scale` name (falling back to the path basename). + +Scales are marked in a first pass so attachments within the same table resolve. +Missing paths are skipped throughout. +""" +function annotate!(parent::Union{HDF5.File,HDF5.Group}, table) + for (path, meta) in table + sc = get(meta, :scale, nothing) + sc === nothing || make_scale!(parent, path, String(sc)) + end + for (path, meta) in table + haskey(parent, path) || continue + a = attrs(parent[path]) + a["long_name"] = String(meta.long_name) + a["units"] = String(get(meta, :units, "1")) + d = get(meta, :dims, nothing) + d === nothing || (a["dims"] = "(" * join(d, ", ") * ")") + att = get(meta, :attach, nothing) + att === nothing && continue + for (axis, scale_path) in att + attach_scale!(parent, path, axis, scale_path, _scale_label(table, scale_path)) + end + end + return parent +end + +# Dimension label for an attachment: the target entry's `scale` name when the table +# declares it, else the scale path's basename. +function _scale_label(table, scale_path) + for (path, meta) in table + path == scale_path || continue + sc = get(meta, :scale, nothing) + sc === nothing || return String(sc) + end + return String(last(split(scale_path, '/'))) +end + +""" + make_scale!(parent, path, name) + +Mark the dataset at `path` as an HDF5 Dimension Scale named `name`. No-op when the +path is absent. +""" +function make_scale!(parent::Union{HDF5.File,HDF5.Group}, path::AbstractString, name::AbstractString) + haskey(parent, path) || return nothing + HDF5.API.h5ds_set_scale(parent[path], String(name)) + return nothing +end + +""" + attach_scale!(parent, path, julia_axis, scale_path, label) + +Attach the Dimension Scale at `scale_path` to Julia axis `julia_axis` (axis 1 first) +of the dataset at `path`, and label that dimension. The H5DS C API indexes file +(row-major) dimensions, so Julia axis `k` of an `N`-d dataset is C index `N - k`. +No-op when either path is absent or the axis lengths disagree. +""" +function attach_scale!(parent::Union{HDF5.File,HDF5.Group}, path::AbstractString, julia_axis::Int, + scale_path::AbstractString, label::AbstractString) + (haskey(parent, path) && haskey(parent, scale_path)) || return nothing + dset = parent[path] + sc = parent[scale_path] + size(dset, julia_axis) == length(sc) || return nothing + cdim = ndims(dset) - julia_axis + HDF5.API.h5ds_attach_scale(dset, sc, cdim) + # h5ds_set_label's wrapper types the C `const char*` as Ref{UInt8}; pass a + # NUL-terminated byte buffer instead of a String. + HDF5.API.h5ds_set_label(dset, cdim, Vector{UInt8}(codeunits(String(label) * "\0"))) + return nothing +end + +""" + write_root_attrs!(file; title) + +Stamp the file-level contract: `schema_version`, `Conventions`, `references`, +`title` (run description), `date_created` (ISO 8601 UTC). The code version lives in +`Info/git_version`. +""" +function write_root_attrs!(file::HDF5.File; title::AbstractString) + a = attrs(file) + a["schema_version"] = "2.0" + a["Conventions"] = "GPEC-HDF5-2.0" + a["references"] = "docs/development/hdf5-conventions.md; https://openfusiontoolkit.github.io/GPEC/dev/" + a["title"] = String(title) + a["date_created"] = Dates.format(Dates.now(UTC), dateformat"yyyy-mm-dd\THH:MM:SS\Z") + return file +end + +end # module HDF5Annotations diff --git a/src/Utilities/Utilities.jl b/src/Utilities/Utilities.jl index 5dab6bdb2..1dfaf8bd1 100644 --- a/src/Utilities/Utilities.jl +++ b/src/Utilities/Utilities.jl @@ -13,6 +13,8 @@ mathematical utilities. - `PhysicalConstants`: SI physical constants matching Fortran GPEC/SLAYER values - `NeoclassicalResistivity`: Spitzer/Sauter/Redl resistivity closures shared by the GGJ and SLAYER inner-layer models + - `HDF5Annotations`: self-describing metadata (long_name/units/dims attributes and + HDF5 Dimension Scales) for the gpec.h5 output """ module Utilities @@ -22,6 +24,7 @@ include("PhysicalConstants.jl") include("KineticProfiles.jl") include("NeoclassicalResistivity.jl") include("GridUtilities.jl") +include("HDF5Annotations.jl") using .FourierTransforms export FourierTransform, inverse, compute_fourier_coefficients @@ -34,6 +37,9 @@ export MU_0, M_E, M_P, E_CHG, K_B, EPS_0 export KineticProfiles +using .HDF5Annotations +export HDF5Annotations + using .NeoclassicalResistivity export NeoclassicalResistivity export NeoResistivityModel, SpitzerModel, SpitzerHarmModel, SauterNeoModel, RedlNeoModel diff --git a/test/h5_metadata_check.jl b/test/h5_metadata_check.jl new file mode 100644 index 000000000..c2e1fb7de --- /dev/null +++ b/test/h5_metadata_check.jl @@ -0,0 +1,27 @@ +# Shared metadata-contract walker (docs/development/hdf5-conventions.md): every +# dataset carries long_name + units, and rank ≥ 2 datasets carry a dims axis-name +# attribute. Exempt: the Input/ raw snapshot and the debug-only GalerkinIntegration +# Match/ group. Included by runtests_h5_schema.jl (full-run walk) and +# runtests_slayer_runner.jl (Tearing/ writer walk, which no full-run test deck covers). + +_metadata_exempt(path) = startswith(path, "Input/") || occursin("/Match/", path) + +function _collect_metadata_violations(h5) + bad = String[] + function walk(node, prefix) + for k in keys(node) + child = node[k] + full = isempty(prefix) ? k : prefix * "/" * k + if child isa HDF5.Group + walk(child, full) + elseif !_metadata_exempt(full) + a = attrs(child) + haskey(a, "long_name") || push!(bad, "$full: missing long_name") + haskey(a, "units") || push!(bad, "$full: missing units") + ndims(child) >= 2 && !haskey(a, "dims") && push!(bad, "$full: missing dims") + end + end + end + walk(h5, "") + return bad +end diff --git a/test/runtests_h5_schema.jl b/test/runtests_h5_schema.jl index b8e4662a3..aba9d32c5 100644 --- a/test/runtests_h5_schema.jl +++ b/test/runtests_h5_schema.jl @@ -32,6 +32,8 @@ function _collect_bad_groups(h5) return bad end +include("h5_metadata_check.jl") + # The full-run walk below only exercises an ideal deck, which never produces the # data-driven group names — pin the whitelist rules directly. @testset "gpec.h5 schema naming: group-name rule" begin @@ -70,6 +72,24 @@ end # Inputs live only under Input/; spot-check the rerun-critical paths. @test haskey(h5, "Input/gpec_toml_raw") @test haskey(h5, "Info/git_version") + + # Metadata contract: long_name/units everywhere, dims on rank ≥ 2 arrays. + viol = _collect_metadata_violations(h5) + isempty(viol) || @error "metadata contract violations in gpec.h5" viol + @test isempty(viol) + + # File-level attributes. + ra = attrs(h5) + for k in ("schema_version", "Conventions", "title", "date_created") + @test haskey(ra, k) + end + @test ra["schema_version"] == "2.0" + + # Dimension scales: the ψ_N coordinate of the forward integration is a + # scale and is attached to its q profile (netCDF-4 pattern). + fwd = "ForceFreeStates/Solutions/ForwardIntegration" + @test HDF5.API.h5ds_is_scale(h5["$fwd/psi"]) + @test HDF5.API.h5ds_is_attached(h5["$fwd/q"], h5["$fwd/psi"], 0) end end end diff --git a/test/runtests_parallel_integration.jl b/test/runtests_parallel_integration.jl index d23a00790..34ccc7256 100644 --- a/test/runtests_parallel_integration.jl +++ b/test/runtests_parallel_integration.jl @@ -265,7 +265,7 @@ using TOML @test isapprox(et_par, et_std; rtol=0.02) # Per-surface Δ' assertions were removed: per-surface Δ' is a stub calculation # left in the code for future work but no longer reported, output, or tested. - # The STRIDE BVP Δ' matrix (`SingularSurfaces/delta_prime_matrix`) is the canonical + # The STRIDE BVP Δ' matrix (`SingularSurfaces/Delta_prime_matrix`) is the canonical # Δ', regression-tested via the DIIID-like fixture which has well-conditioned # values; Solovev is near marginal stability and BVP Δ' is pathological there. end diff --git a/test/runtests_slayer_runner.jl b/test/runtests_slayer_runner.jl index f74d774c9..25a28f4fc 100644 --- a/test/runtests_slayer_runner.jl +++ b/test/runtests_slayer_runner.jl @@ -5,6 +5,8 @@ using GeneralizedPerturbedEquilibrium.Runner using HDF5 + include("h5_metadata_check.jl") + # ------- Helper: build a synthetic SLAYERParameters with full control function _mk_params(; rs=0.5, lu=1e7, tauk=1e-4, Q_e=-1.0, Q_i=0.5, m=2, n=1, ising=1, @@ -173,6 +175,16 @@ h5open(path, "w") do f write_slayer_hdf5!(f, r) end + + # Metadata contract must hold for Tearing/ too — the full-run schema test + # only exercises an ideal deck, which writes no Tearing/ group. + h5open(path, "r") do f + viol = _collect_metadata_violations(f) + isempty(viol) || @error "metadata contract violations in Tearing/" viol + @test isempty(viol) + @test attrs(f["Tearing"])["layer_model"] == "slayer" + end + h5open(path, "r") do f g = f["Tearing"] @test haskey(g, "enabled") && read(g["enabled"]) == 1 @@ -184,25 +196,24 @@ @test haskey(g, "Scan") # Per-surface arrays have the right length - @test length(read(g["PerSurface/ising"])) == 2 - @test read(g["PerSurface/ising"]) == [1, 2] + @test length(read(g["PerSurface/rational_index"])) == 2 + @test read(g["PerSurface/rational_index"]) == [1, 2] @test read(g["PerSurface/lu"])[1] ≈ 1.0e7 @test read(g["PerSurface/lu"])[2] ≈ 2.0e7 # Roots arrays - @test length(read(g["Roots/Q_root_real"])) == 1 # coupled - @test length(read(g["Roots/omega_Hz"])) == 1 + @test length(read(g["Roots/Q_root"])) == 1 # coupled + @test length(read(g["Roots/omega"])) == 1 # Layer-thickness diagnostic: one entry per surface, with # the physical thickness [m] and the drift scale. - @test length(read(g["LayerWidths/delta_s_m"])) == 2 - @test all(read(g["LayerWidths/delta_s_m"]) .>= 0) - @test haskey(g["LayerWidths"], "dels_db_real") + @test length(read(g["LayerWidths/delta_s_abs"])) == 2 + @test all(read(g["LayerWidths/delta_s_abs"]) .>= 0) + @test haskey(g["LayerWidths"], "delta_s_over_d_beta") @test haskey(g["LayerWidths"], "d_beta") # Ragged diagnostics use flat+offsets encoding - @test haskey(g["Diagnostics/ValidRoots"], "flat_real") - @test haskey(g["Diagnostics/ValidRoots"], "flat_imag") + @test haskey(g["Diagnostics/ValidRoots"], "flat") @test haskey(g["Diagnostics/ValidRoots"], "offsets") # Scan group present (store_scan=true) @@ -212,6 +223,26 @@ end end + @testset "GGJ per-surface writer: metadata contract" begin + # The GGJ branch of _write_per_surface! is unreachable from the SLAYER + # round-trip above; enforce its annotation-table rows on a synthetic write. + params = [GGJParameters(; E=0.1, F=0.2, G=0.3, H=0.4, K=0.5, taua=1e-6, taur=1.0, ising=i) for i in 1:2] + dp = ComplexF64[1.0 0.0; 0.0 2.0] + mktemp() do path, io + close(io) + h5open(path, "w") do f + g = HDF5.create_group(f, "Tearing") + Runner._write_per_surface!(g, params, dp) + Runner._annotate_tearing!(g) + end + h5open(path, "r") do f + viol = _collect_metadata_violations(f) + isempty(viol) || @error "metadata contract violations in GGJ PerSurface/" viol + @test isempty(viol) + end + end + end + @testset "write_slayer_hdf5!: disabled result still emits enabled=0" begin c = SLAYERControl(; enabled=false) r = empty_slayer_result(c) diff --git a/test/test_data/regression_solovev_ideal_example/gpec.toml b/test/test_data/regression_solovev_ideal_example/gpec.toml index a16ab6898..f2fc6d1f0 100644 --- a/test/test_data/regression_solovev_ideal_example/gpec.toml +++ b/test/test_data/regression_solovev_ideal_example/gpec.toml @@ -46,7 +46,7 @@ singfac_min = 1e-4 # Fractional distance from rational q at which ide ucrit = 1e3 # Column-norm threshold that triggers solution renormalization # Δ' BVP + parallel integration (see ForceFreeStatesControl docstring for details) -use_parallel = true # Run parallel FM-propagator BVP path (unlocks SingularSurfaces/delta_prime_matrix) +use_parallel = true # Run parallel FM-propagator BVP path (unlocks SingularSurfaces/Delta_prime_matrix) parallel_threads = 2 # BVP thread cap (1 = serial/bit-deterministic; 2 ≈ +20% speedup; ≥3 saturates) populate_dense_xi = true # Append serial-EL pass so dense ξ is stored — REQUIRED with a [PerturbedEquilibrium] section set_psilim_via_dmlim = false # FALSE for limited/analytical equilibria — rationals sparse, dmlim would chop too much edge diff --git a/test/test_data/regression_solovev_kinetic_example/gpec.toml b/test/test_data/regression_solovev_kinetic_example/gpec.toml index 1c020a086..cef6ad546 100644 --- a/test/test_data/regression_solovev_kinetic_example/gpec.toml +++ b/test/test_data/regression_solovev_kinetic_example/gpec.toml @@ -46,7 +46,7 @@ singfac_min = 1e-4 # Fractional distance from rational q at which ide ucrit = 1e3 # Column-norm threshold that triggers solution renormalization # Δ' BVP + parallel integration (see ForceFreeStatesControl docstring for details) -use_parallel = true # Run parallel FM-propagator BVP path (unlocks SingularSurfaces/delta_prime_matrix) +use_parallel = true # Run parallel FM-propagator BVP path (unlocks SingularSurfaces/Delta_prime_matrix) parallel_threads = 2 # BVP thread cap (1 = serial/bit-deterministic; 2 ≈ +20% speedup; ≥3 saturates) populate_dense_xi = true # Append serial-EL pass so dense ξ is stored — REQUIRED with a [PerturbedEquilibrium] section set_psilim_via_dmlim = false # FALSE for limited/analytical equilibria — rationals sparse, dmlim would chop too much edge